1SUBROUTINE LinearFormsAssembly( Model,Solver,dt,TransientSimulation )
2!------------------------------------------------------------------------------
3!******************************************************************************
4!
5!  Unit test for linear form / vectorized basis function computation
6!  routines in Elmer
7!
8!  ARGUMENTS:
9!
10!  TYPE(Model_t) :: Model,
11!     INPUT: All model information (mesh, materials, BCs, etc...)
12!
13!  TYPE(Solver_t) :: Solver
14!     INPUT: Linear & nonlinear equation solver options
15!
16!  REAL(KIND=dp) :: dt,
17!     INPUT: Timestep size for time dependent simulations
18!
19!  LOGICAL :: TransientSimulation
20!     INPUT: Steady state or transient simulation
21!
22!******************************************************************************
23    USE DefUtils
24    USE LinearForms
25    USE ISO_C_BINDING
26#ifdef _OPENMP
27    USE omp_lib
28#endif
29!------------------------------------------------------------------------------
30    IMPLICIT NONE
31!------------------------------------------------------------------------------
32    TYPE(Model_t) :: Model
33    TYPE(Solver_t) :: Solver
34
35    REAL(KIND=dp) :: dt
36    LOGICAL :: TransientSimulation
37!------------------------------------------------------------------------------
38! Local variables
39!------------------------------------------------------------------------------
40    REAL(KIND=dp), PARAMETER :: tol1d = 1D-12, tol2d=1D-12, tol3d=1D-12
41    REAL(KIND=dp) :: float_P
42    INTEGER :: nerror, netest, P
43    LOGICAL :: Found
44
45    nerror = 0
46    float_P = ListGetCReal(GetSolverParams(), 'P', Found)
47    IF (.NOT. Found) float_P = 6.0_dp
48    P = NINT(float_P)
49
50    ! 1D tests
51    netest = TestLineElement(Solver, P, tol1d)
52    IF (netest /= 0) THEN
53       CALL Warn('LinearFormsAssembly','Line element contained errors')
54    END IF
55    nerror = nerror + netest
56
57    ! ! 2D tests
58    netest = TestTriangleElement(Solver, P, tol2d)
59    IF (netest /= 0) THEN
60      CALL Warn('LinearFormsAssembly','Triangle element contained errors')
61    END IF
62    nerror = nerror + netest
63
64    netest = TestQuadElement(Solver, P, tol2d)
65    IF (netest /= 0) THEN
66      CALL Warn('LinearFormsAssembly','Quad element contained errors')
67    END IF
68    nerror = nerror + netest
69
70    ! ! 3D tests
71    netest = TestTetraElement(Solver, P, tol3d)
72    IF (netest /= 0) THEN
73      CALL Warn('LinearFormsAssembly','Tetra element contained errors')
74    END IF
75    nerror = nerror + netest
76
77    netest = TestWedgeElement(Solver, P, tol3d)
78    IF (netest /= 0) THEN
79      CALL Warn('LinearFormsAssembly','Wedge element contained errors')
80    END IF
81    nerror = nerror + netest
82
83    netest = TestBrickElement(Solver, P, tol3d)
84    IF (netest /= 0) THEN
85      CALL Warn('LinearFormsAssembly','Brick element contained errors')
86    END IF
87    nerror = nerror + netest
88
89    ! Build solution norm for error checking
90    Solver % Variable % Norm = REAL(1+nerror,dp)
91    Solver % Variable % Values = REAL(1+nerror,dp)
92
93CONTAINS
94
95  FUNCTION TestLineElement(Solver, P, tol) RESULT(nerror)
96    IMPLICIT NONE
97
98    TYPE(Solver_t) :: Solver
99    INTEGER, INTENT(IN) :: P
100    REAL(kind=dp), INTENT(IN) :: tol
101    INTEGER :: nerror
102
103    nerror = TestElement(Solver, 202, P, tol)
104  END FUNCTION TestLineElement
105
106  FUNCTION TestTriangleElement(Solver, P, tol) RESULT(nerror)
107    IMPLICIT NONE
108
109    TYPE(Solver_t) :: Solver
110    INTEGER, INTENT(IN) :: P
111    REAL(kind=dp), INTENT(IN) :: tol
112    INTEGER :: nerror
113
114    nerror = TestElement(Solver, 303, P, tol)
115  END FUNCTION TestTriangleElement
116
117  FUNCTION TestQuadElement(Solver, P, tol) RESULT(nerror)
118    IMPLICIT NONE
119
120    TYPE(Solver_t) :: Solver
121    INTEGER, INTENT(IN) :: P
122    REAL(kind=dp), INTENT(IN) :: tol
123    INTEGER :: nerror
124
125    nerror = TestElement(Solver, 404, P, tol)
126  END FUNCTION TestQuadElement
127
128  FUNCTION TestTetraElement(Solver, P, tol) RESULT(nerror)
129    IMPLICIT NONE
130
131    TYPE(Solver_t) :: Solver
132    INTEGER, INTENT(IN) :: P
133    REAL(kind=dp), INTENT(IN) :: tol
134    INTEGER :: nerror
135
136    nerror = TestElement(Solver, 504, P, tol)
137  END FUNCTION TestTetraElement
138
139  FUNCTION TestWedgeElement(Solver, P, tol) RESULT(nerror)
140    IMPLICIT NONE
141
142    TYPE(Solver_t) :: Solver
143    INTEGER, INTENT(IN) :: P
144    REAL(kind=dp), INTENT(IN) :: tol
145    INTEGER :: nerror
146
147    nerror = TestElement(Solver, 706, P, tol)
148  END FUNCTION TestWedgeElement
149
150  FUNCTION TestBrickElement(Solver, P, tol) RESULT(nerror)
151    IMPLICIT NONE
152
153    TYPE(Solver_t) :: Solver
154    INTEGER, INTENT(IN) :: P
155    REAL(kind=dp), INTENT(IN) :: tol
156    INTEGER :: nerror
157
158    nerror = TestElement(Solver, 808, P, tol)
159  END FUNCTION TestBrickElement
160
161  FUNCTION TestElement(Solver, ecode, P, tol) RESULT(nerror)
162    IMPLICIT NONE
163
164    TYPE(Solver_t) :: Solver
165    INTEGER, INTENT(IN) :: ecode
166    INTEGER, INTENT(IN) :: P
167    REAL(kind=dp), INTENT(IN) :: tol
168
169    TYPE(Element_t), POINTER :: Element, SingleElement
170    TYPE(Mesh_t), POINTER :: NewMesh, OldMesh
171    REAL(KIND=dp), ALLOCATABLE :: STIFF(:,:), LOAD(:), FORCE(:), &
172            STIFFvec(:,:), FORCEvec(:)
173
174    INTEGER :: i, j, k, l, q, nerror, nbasis, nndof, allocstat, tag, nthr, &
175            nbasisvec, ndbasisdxvec, rep, dim, lm_eval, lm_eval_vec, NumGP
176
177    INTEGER, PARAMETER :: NREP = 100
178    REAL(kind=dp) :: t_start, t_end, t_tot, t_startvec, t_endvec, t_tot_vec
179    TYPE(GaussIntegrationPoints_t) :: Quadrature
180
181    nerror = 0
182    lm_eval = 0
183    lm_eval_vec = 0
184    t_tot = REAL(0,dp)
185    t_tot_vec = REAL(0,dp)
186
187    ! Create a mesh with a single element
188    OldMesh => Solver % Mesh
189    NewMesh => NULL()
190    CALL AllocateMeshAndPElement(NewMesh, ecode, P, SingleElement)
191    Solver % Mesh => NewMesh
192
193    ! Insert P element definitions to Solver mapping (sets P elements as "active")
194    IF (ALLOCATED(Solver % Def_Dofs)) THEN
195      tag = ecode / 100
196      Solver % Def_Dofs(tag,1,6) = P
197    END IF
198    IF (ecode==404) THEN
199      Quadrature = GaussPointsQuad((P+1)**2)
200    ELSE
201      Quadrature = GaussPoints(SingleElement, PReferenceElement=.TRUE.)
202    END IF
203    NumGP = Quadrature % N
204
205    !$OMP PARALLEL SHARED(Solver, SingleElement, ecode, tol) &
206    !$OMP PRIVATE(STIFF, FORCE, STIFFvec, FORCEvec, &
207    !$OMP         LOAD, Element, nndof, nbasis, &
208    !$OMP         allocstat, rep, t_start, t_end, &
209    !$OMP         t_startvec, t_endvec) &
210    !$OMP REDUCTION(+:nerror,t_tot,t_tot_vec,lm_eval,lm_eval_vec) &
211    !$OMP DEFAULT(NONE)
212
213    ! Construct a temporary element based on the one created previously
214    Element => ClonePElement(SingleElement)
215
216    nndof = Element % Type % NumberOfNodes
217    nbasis = nndof + GetElementNOFDOFs( Element )
218
219    ! Reserve workspace
220    ALLOCATE(STIFF(nbasis, nbasis), FORCE(nbasis), &
221            STIFFvec(nbasis, nbasis), FORCEvec(nbasis), &
222            LOAD(nndof), &
223            STAT=allocstat)
224    IF (allocstat /= 0) THEN
225      CALL Fatal('LinearForms',&
226              'Storage allocation for local matrices failed')
227    END IF
228
229    ! Initialize artificial load vector
230    LOAD = REAL(1,dp)
231
232    ! Warmup
233    CALL LocalMatrix( STIFF, FORCE, LOAD, Element, nndof, nbasis)
234    !$OMP BARRIER
235    t_start = ftimer()
236    DO rep=1,NREP
237      ! Construct local matrix
238!DIR$ NOINLINE
239      CALL LocalMatrix( STIFF, FORCE, LOAD, Element, nndof, nbasis)
240    END DO
241    t_end = ftimer()
242    lm_eval = NREP
243
244    ! Warmup
245    CALL LocalMatrixVec( STIFFvec, FORCEvec, LOAD, Element, nndof, nbasis)
246    !$OMP BARRIER
247    t_startvec = ftimer()
248    DO rep=1,NREP
249      ! Construct local matrix
250!DIR$ NOINLINE
251      CALL LocalMatrixVec( STIFFvec, FORCEvec, LOAD, Element, nndof, nbasis)
252    END DO
253    t_endvec = ftimer()
254    lm_eval_vec = NREP
255
256    nerror = TestLocalMatrix(nbasis, STIFF, STIFFvec, tol)
257    nerror = nerror + TestLocalForce(nbasis, FORCE, FORCEvec, tol)
258
259    t_tot = t_end - t_start
260    t_tot_vec = t_endvec - t_startvec
261
262    CALL DeallocatePElement(Element)
263    DEALLOCATE(STIFF, FORCE, STIFFvec, FORCEvec, LOAD)
264
265    !$OMP END PARALLEL
266
267    ! Normalize the times / thread
268    nthr = 1
269    !$ nthr = omp_get_max_threads()
270    t_tot = t_tot / nthr
271    t_tot_vec = t_tot_vec / nthr
272    CALL PrintTestData(SingleElement, NumGP, t_tot, lm_eval, &
273            t_tot_vec, lm_eval_vec)
274
275    CALL DeallocateTemporaryMesh(NewMesh)
276    Solver % Mesh => OldMesh
277    CALL DeallocatePElement(SingleElement)
278
279    IF (ALLOCATED(Solver % Def_Dofs)) THEN
280      tag = ecode / 100
281      Solver % Def_Dofs(tag,1,6) = 0
282    END IF
283  END FUNCTION TestElement
284
285  SUBROUTINE LocalMatrix( STIFF, FORCE, LOAD, Element, n, nd )
286    IMPLICIT NONE
287    REAL(KIND=dp) CONTIG :: STIFF(:,:), FORCE(:)
288    REAL(KIND=dp) CONTIG, INTENT(IN) :: LOAD(:)
289    INTEGER :: n, nd
290    TYPE(Element_t), POINTER :: Element
291!------------------------------------------------------------------------------
292    REAL(KIND=dp) :: Basis(nd),dBasisdx(nd,3),DetJ,LoadAtIP,Weight
293    LOGICAL :: Stat
294    INTEGER :: i,j,t,p,q,dim
295    TYPE(GaussIntegrationPoints_t) :: IP
296
297    TYPE(Nodes_t), SAVE :: Nodes
298    !$OMP THREADPRIVATE(Nodes)
299!------------------------------------------------------------------------------
300    CALL GetReferenceElementNodes( Nodes, Element )
301    STIFF = 0.0d0
302    FORCE = 0.0d0
303
304    dim = Element % Type % Dimension
305    !Numerical integration:
306    !----------------------
307    IF (Element % TYPE % ElementCode / 100 == 4) THEN
308      IP = GaussPointsQuad((Element % PDefs % P+1)**2)
309    ELSE
310      IP = GaussPoints( Element, PReferenceElement=.TRUE. )
311    END IF
312
313    DO t=1,IP % n
314
315      ! Basis function values & derivatives at the integration point:
316      !--------------------------------------------------------------
317      stat = ElementInfo( Element, Nodes, IP % U(t), IP % V(t), &
318                           IP % W(t), detJ, Basis, dBasisdx)
319
320      ! The source term at the integration point:
321      !------------------------------------------
322      LoadAtIP = SUM( Basis(1:n) * LOAD(1:n) )
323      Weight = IP % s(t) * DetJ
324      ! STIFF=STIFF+(grad u, grad v)
325      STIFF(1:nd,1:nd) = STIFF(1:nd,1:nd) + Weight * &
326            MATMUL( dBasisdx(1:nd,1:dim), TRANSPOSE( dBasisdx(1:nd,1:dim) ) )
327
328      DO p=1,nd
329        DO q=1,nd
330          ! STIFF=STIFF+(grad u,v)
331          ! -----------------------------------
332          STIFF (p,q) = STIFF(p,q) + Weight * SUM(dBasisdx(q,1:dim)) * Basis(p)
333
334          ! STIFF=STIFF+(u,v)
335          STIFF(p,q) = STIFF(p,q) + Weight * Basis(q) * Basis(p)
336        END DO
337      END DO
338
339      FORCE(1:nd) = FORCE(1:nd) + IP % s(t) * DetJ * LoadAtIP * Basis(1:nd)
340    END DO
341  END SUBROUTINE LocalMatrix
342
343  SUBROUTINE LocalMatrixVec( STIFF, FORCE, LOAD, Element, n, nd )
344!------------------------------------------------------------------------------
345    IMPLICIT NONE
346    REAL(KIND=dp) CONTIG :: STIFF(:,:), FORCE(:)
347    REAL(KIND=dp) CONTIG, INTENT(IN) :: LOAD(:)
348    INTEGER, INTENT(IN) :: n, nd
349    TYPE(Element_t), POINTER :: Element
350!------------------------------------------------------------------------------
351    LOGICAL :: Stat
352    INTEGER :: i, j, ngp, allocstat, gp
353    TYPE(GaussIntegrationPoints_t) :: IP
354
355    TYPE(Nodes_t), SAVE :: Nodes
356    REAL(KIND=dp), ALLOCATABLE, SAVE :: Basis(:,:), dBasisdx(:,:,:), &
357            DetJ(:), LoadAtIPs(:)
358    !$OMP THREADPRIVATE(Nodes, Basis, dBasisdx, DetJ, LoadAtIPs)
359!------------------------------------------------------------------------------
360    CALL GetReferenceElementNodes( Nodes, Element )
361
362    STIFF = REAL(0, dp)
363    FORCE = REAL(0, dp)
364
365    ! Get integration points
366    IF (Element % TYPE % ElementCode / 100 == 4) THEN
367      IP = GaussPointsQuad((Element % PDefs % P+1)**2)
368    ELSE
369      IP = GaussPoints( Element, PReferenceElement=.TRUE. )
370    END IF
371    ngp = IP % n
372
373    ! Reserve workspace
374    IF (.NOT. ALLOCATED(Basis)) THEN
375      ALLOCATE(Basis(ngp,nd), dBasisdx(ngp,nd,3), &
376              DetJ(ngp), LoadAtIPs(ngp), STAT=allocstat)
377      IF (allocstat /= 0) THEN
378        CALL Fatal('LocalMatrixVec',&
379                   'Storage allocation for local element basis failed')
380      END IF
381    ELSE IF (SIZE(Basis,1) /= ngp .OR. SIZE(Basis,2) /= nd) THEN
382      DEALLOCATE(Basis, dBasisdx, DetJ, LoadAtIPs)
383      ALLOCATE(Basis(ngp,nd), dBasisdx(ngp,nd,3), &
384              DetJ(ngp), LoadAtIPs(ngp), STAT=allocstat)
385      IF (allocstat /= 0) THEN
386        CALL Fatal('LocalMatrixVec',&
387                'Storage allocation for local element basis failed')
388      END IF
389    END IF
390
391    ! Compute values of all basis functions at all integration points
392    stat = ElementInfoVec( Element, Nodes, ngp, &
393            IP % U, IP % V, IP % W, DetJ, SIZE(Basis,2), Basis, dBasisdx )
394
395
396    ! Compute actual integration weights (recycle memory space of DetJ)
397    DO i=1,ngp
398       DetJ(i) = Ip % s(i)*Detj(i)
399    END DO
400
401    ! STIFF=STIFF+(grad u, grad u)
402    CALL LinearForms_GradUdotGradU(ngp, nd, Element % TYPE % DIMENSION, &
403             dBasisdx, DetJ, STIFF)
404    ! STIFF=STIFF+(u,u)
405    CALL LinearForms_UdotU(ngp, nd, Element % TYPE % DIMENSION, Basis, DetJ, STIFF)
406    ! STIFF=STIFF+(grad u,v)
407    CALL LinearForms_GradUdotU(ngp, nd, Element % TYPE % DIMENSION, dBasisdx, Basis, DetJ, STIFF)
408
409    ! Source terms at IPs
410    !------------------------------------------
411    ! LoadAtIPs(1:ngp) = MATMUL( Basis(1:ngp,1:n), LOAD(1:n) )
412    CALL LinearForms_ProjectToU(ngp, n, Basis, LOAD, LoadAtIPs)
413
414    ! FORCE=FORCE+(u,f)
415    CALL LinearForms_UdotF(ngp, nd, Basis, DetJ, LoadAtIPs, FORCE)
416  END SUBROUTINE LocalMatrixVec
417
418  FUNCTION TestLocalMatrix(nbasis, STIFF1, STIFF2, tol) RESULT(nerror)
419    IMPLICIT NONE
420
421    INTEGER, INTENT(IN) :: nbasis
422    REAL(KIND=dp) CONTIG, INTENT(IN) :: STIFF1(:,:), STIFF2(:,:)
423    REAL(kind=dp), INTENT(IN) :: tol
424    INTEGER :: nerror
425
426    INTEGER :: i, j, dim
427
428    nerror = 0
429    ! Test element of local matrix
430    DO j=1,nbasis
431      DO i=1,nbasis
432        IF (ABS(STIFF1(i,j)-STIFF2(i,j)) >= tol) THEN
433          nerror = nerror + 1
434          WRITE (*,*) 'STIFF:', i,j,STIFF1(i,j), STIFF2(i,j)
435        END IF
436      END DO
437    END DO
438  END FUNCTION TestLocalMatrix
439
440  FUNCTION TestLocalForce(nbasis, FORCE1, FORCE2, tol) RESULT(nerror)
441    IMPLICIT NONE
442
443    INTEGER, INTENT(IN) :: nbasis
444    REAL(KIND=dp) CONTIG, INTENT(IN) :: FORCE1(:), FORCE2(:)
445    REAL(kind=dp), INTENT(IN) :: tol
446    INTEGER :: nerror
447
448    INTEGER :: i, dim
449
450    nerror = 0
451    ! Test element of local force vector
452    DO i=1,nbasis
453      IF (ABS(FORCE1(i)-FORCE2(i)) >= tol) THEN
454        nerror = nerror + 1
455        WRITE (*,*) 'FORCE:', i, FORCE1(i), FORCE2(i)
456      END IF
457    END DO
458  END FUNCTION TestLocalForce
459
460  SUBROUTINE GetReferenceElementNodes( ElementNodes, Element )
461     TYPE(Nodes_t), TARGET :: ElementNodes
462     TYPE(Element_t) :: Element
463
464     INTEGER :: i, n, padn, sz, astat
465
466     n = Element % TYPE % NumberOfNodes
467     padn = n
468
469     ! TODO: Implement padding
470     IF (.NOT. ALLOCATED( ElementNodes % xyz)) THEN
471       ! Deallocate old storage
472       IF (ASSOCIATED(ElementNodes % x)) DEALLOCATE(ElementNodes % x)
473       IF (ASSOCIATED(ElementNodes % y)) DEALLOCATE(ElementNodes % y)
474       IF (ASSOCIATED(ElementNodes % z)) DEALLOCATE(ElementNodes % z)
475
476       ! Allocate new storage
477       ALLOCATE(ElementNodes % xyz(padn,3))
478       ElementNodes % xyz = REAL(0,dp)
479       ElementNodes % x => ElementNodes % xyz(1:n,1)
480       ElementNodes % y => ElementNodes % xyz(1:n,2)
481       ElementNodes % z => ElementNodes % xyz(1:n,3)
482     ELSE IF (SIZE(ElementNodes % xyz, 1)<padn) THEN
483       DEALLOCATE(ElementNodes % xyz)
484       ALLOCATE(ElementNodes % xyz(padn,3))
485       ElementNodes % xyz = REAL(0,dp)
486       ElementNodes % x => ElementNodes % xyz(1:n,1)
487       ElementNodes % y => ElementNodes % xyz(1:n,2)
488       ElementNodes % z => ElementNodes % xyz(1:n,3)
489     ELSE
490       ElementNodes % x => ElementNodes % xyz(1:n,1)
491       ElementNodes % y => ElementNodes % xyz(1:n,2)
492       ElementNodes % z => ElementNodes % xyz(1:n,3)
493       sz = SIZE(ElementNodes % xyz,1)
494       SELECT CASE(Element % TYPE % DIMENSION)
495       CASE(1)
496         DO i=n+1,sz
497           ElementNodes % xyz(i,1) = REAL(0,dp)
498         END DO
499         DO i=1,sz
500           ElementNodes % xyz(i,2) = REAL(0,dp)
501           ElementNodes % xyz(i,3) = REAL(0,dp)
502         END DO
503       CASE(2)
504         DO i=n+1,sz
505           ElementNodes % xyz(i,1) = REAL(0,dp)
506           ElementNodes % xyz(i,2) = REAL(0,dp)
507         END DO
508         DO i=1,sz
509           ElementNodes % xyz(i,3) = REAL(0,dp)
510         END DO
511       CASE(3)
512         DO i=n+1,sz
513           ElementNodes % xyz(i,1) = REAL(0,dp)
514           ElementNodes % xyz(i,2) = REAL(0,dp)
515           ElementNodes % xyz(i,3) = REAL(0,dp)
516         END DO
517       CASE DEFAULT
518         CALL Fatal('GetReferenceElementNodes','Unsupported element dimension')
519       END SELECT
520     END IF
521
522     IF (isPElement(Element)) THEN
523       CALL GetRefPElementNodes(Element % Type, &
524               ElementNodes % x, &
525               ElementNodes % y, &
526               ElementNodes % z)
527     ELSE
528       IF (ALLOCATED(Element % Type % NodeU)) THEN
529         !DIR$ IVDEP
530         DO i=1,n
531           ElementNodes % x(i) = Element % Type % NodeU(i)
532         END DO
533       END IF
534       IF (ALLOCATED(Element % Type % NodeV)) THEN
535         !DIR$ IVDEP
536         DO i=1,n
537           ElementNodes % y(i) = Element % Type % NodeV(i)
538         END DO
539       END IF
540       IF (ALLOCATED(Element % Type % NodeW)) THEN
541         !DIR$ IVDEP
542         DO i=1,n
543           ElementNodes % z(i) = Element % Type % NodeW(i)
544         END DO
545       END IF
546     END IF
547  END SUBROUTINE GetReferenceElementNodes
548
549  SUBROUTINE PrintTestData(Element, ngp, t_n1, evals1, t_n2, evals2)
550
551    IMPLICIT NONE
552    TYPE(Element_t) :: Element
553    INTEGER, INTENT(IN) :: ngp, evals1, evals2
554    REAL(kind=dp), INTENT(IN) :: t_n1, t_n2
555
556    WRITE (*,'(A,I0)') 'Element type=', Element % TYPE % ElementCode
557    IF (ASSOCIATED(Element % PDefs)) THEN
558      WRITE (*,'(A,I0)') 'Element polynomial degree=', Element % PDefs % P
559    END IF
560    WRITE (*,'(A,L1)') 'Active P element=', isActivePElement(Element)
561    WRITE (*,'(A,I0)') 'Element number of nodes=', Element % TYPE % NumberOfNodes
562    WRITE (*,'(A,I0)') 'Element number of nonnodal dofs=', GetElementNOFDOFs(Element)
563    WRITE (*,'(A,I0)') 'Element number of bubble dofs=', GetElementNOFBDOFs()
564    WRITE (*,'(A,I0)') 'Number of Gauss points=', ngp
565    WRITE (*,'(A,I0)') 'Nodal basis, number of local matrix evaluations=', evals1
566    WRITE (*,'(A,F12.9)') 'Nodal  basis, local matrix assembly t(s):', t_n1
567    WRITE (*,'(A,F12.2)') 'Nodal  basis, local matrix evaluations/sec:', evals1/t_n1
568    WRITE (*,'(A,I0)') 'Vector basis, number of local matrix evaluations=', evals2
569    WRITE (*,'(A,F12.9)') 'Vector basis, local matrix assembly t(s):', t_n2
570    WRITE (*,'(A,F12.2)') 'Vector basis, local matrix evaluations/sec:', evals2/t_n2
571
572  END SUBROUTINE PrintTestData
573
574  SUBROUTINE AllocateMeshAndPElement(Mesh, ElementCode, P, PElement)
575    IMPLICIT NONE
576
577    TYPE(Mesh_t), POINTER :: Mesh
578    INTEGER, INTENT(IN) :: ElementCode, P
579    TYPE(Element_t), POINTER :: PElement
580    INTEGER :: node, edge, face, astat
581
582    ! Allocate a mesh with a single element
583    Mesh => AllocateMesh()
584    ALLOCATE(Mesh % Elements(1), STAT=astat)
585    IF (astat /= 0) THEN
586      CALL Fatal('AllocateMeshAndElement','Allocation of mesh element failed')
587    END IF
588    Mesh % NumberOfBulkElements = 1
589
590    ! Construct P element
591    PElement => AllocateElement()
592    PElement % ElementIndex = 1
593    PElement % BodyId = 1
594    PElement % Type => GetElementType( ElementCode )
595    CALL AllocatePDefinitions(PElement)
596
597    ! Add element node indexes to element
598    ALLOCATE(PElement % NodeIndexes(PElement % Type % NumberOfNodes), STAT=astat)
599    IF (astat /= 0) THEN
600      CALL Fatal('AllocateMeshAndElement','Allocation of node indices failed')
601    END IF
602    PElement % NodeIndexes(:) = [(node, node=1,PElement % Type % NumberOfNodes)]
603
604    IF (PElement % Type % Dimension > 1) THEN
605      ! Add element edge indexes to element
606      ALLOCATE(PElement % EdgeIndexes(PElement % Type % NumberOfEdges), STAT=astat)
607      IF (astat /= 0) THEN
608        CALL Fatal('AllocateMeshAndElement','Allocation of edge indices failed')
609      END IF
610      PElement % EdgeIndexes(:) = [(edge, edge=1,PElement % Type % NumberOfEdges)]
611
612      ! Add all element edges to mesh
613      ALLOCATE(Mesh % Edges(PElement % Type % NumberOfEdges), STAT=astat)
614      IF (astat /= 0) THEN
615        CALL Fatal('AllocateMeshAndElement','Allocation of mesh edges failed')
616      END IF
617      Mesh % NumberOfEdges = PElement % Type % NumberOfEdges
618      ! Mesh % MinEdgeDofs = HUGE(Mesh % MinEdgeDofs)
619      Mesh % MinEdgeDofs = 0
620      Mesh % MaxEdgeDofs = 0
621      DO edge=1,PElement % Type % NumberOfEdges
622        CALL InitializePElement(Mesh % Edges(edge))
623        ! Allocate edge element (always 202, even in 3D)
624        Mesh % Edges(edge) % Type => GetElementType(202, .FALSE.)
625        CALL AllocatePDefinitions(Mesh % Edges(edge))
626        Mesh % Edges(edge) % PDefs % P = P
627        Mesh % Edges(edge) % PDefs % isEdge = .TRUE.
628        Mesh % Edges(edge) % PDefs % GaussPoints = (P+1) ** Mesh % Edges(edge) % Type % DIMENSION
629        Mesh % Edges(edge) % PDefs % LocalNumber = edge
630        Mesh % Edges(edge) % BDOFs = GetBubbleDofs(Mesh % Edges(edge), P)
631        Mesh % MinEdgeDofs = MIN(Mesh % MinEdgeDofs, Mesh % Edges(edge) % BDOFs)
632        Mesh % MaxEdgeDofs = MAX(Mesh % MaxEdgeDofs, Mesh % Edges(edge) % BDOFs)
633      END DO
634    END IF
635
636    IF (PElement % Type % Dimension > 2) THEN
637      ! Add element face indexes to element
638      ALLOCATE(PElement % FaceIndexes(PElement % Type % NumberOfFaces))
639      IF (astat /= 0) THEN
640        CALL Fatal('AllocateMeshAndElement','Allocation of face indices failed')
641      END IF
642      PElement % FaceIndexes(:) = [(face, face=1,PElement % Type % NumberOfFaces)]
643
644      ! Add element faces to mesh
645      ALLOCATE(Mesh % Faces(PElement % Type % NumberOfFaces), STAT=astat)
646      IF (astat /= 0) THEN
647        CALL Fatal('AllocateMeshAndElement','Allocation of mesh faces failed')
648      END IF
649      Mesh % NumberOfFaces = PElement % Type % NumberOfFaces
650      Mesh % MinFaceDofs = HUGE(Mesh % MinFaceDofs)
651      Mesh % MaxFaceDofs = 0
652      DO face=1,PElement % Type % NumberOfFaces
653        CALL InitializePElement(Mesh % Faces(face))
654        ! Allocate edge element (always 303 or 404, depending on element)
655        SELECT CASE (ElementCode/100)
656        CASE(5)
657          Mesh % Faces(face) % Type => GetElementType(303, .FALSE.)
658        CASE(6)
659          IF ( face == 1 ) THEN
660            Mesh % Faces(Face) % Type => GetElementType( 404, .FALSE. )
661          ELSE
662            Mesh % Faces(Face) % Type => GetElementType( 303, .FALSE. )
663          END IF
664        CASE(7)
665          IF ( face <= 2 ) THEN
666            Mesh % Faces(Face) % Type => GetElementType( 303, .FALSE. )
667          ELSE
668            Mesh % Faces(Face) % Type => GetElementType( 404, .FALSE. )
669          END IF
670        CASE(8)
671          Mesh % Faces(Face) % Type => GetElementType( 404, .FALSE.)
672        CASE DEFAULT
673          CALL Fatal('AllocateMeshAndElement','Unknown element type')
674        END SELECT
675        CALL AllocatePDefinitions(Mesh % Faces(face))
676        Mesh % Faces(face) % PDefs % P = P
677        Mesh % Faces(face) % PDefs % isEdge = .TRUE.
678        Mesh % Faces(face) % PDefs % GaussPoints = (P+1) ** Mesh % Faces(face) % Type % DIMENSION
679        Mesh % Faces(face) % PDefs % LocalNumber = face
680        Mesh % Faces(face) % BDOFs = GetBubbleDofs(Mesh % Faces(face), P)
681        Mesh % MinFaceDofs = MAX(Mesh % MinFaceDofs, Mesh % Faces(face) % BDOFs)
682        Mesh % MaxFaceDofs = MAX(Mesh % MaxFaceDofs, Mesh % Faces(face) % BDOFs)
683      END DO
684    END IF
685
686    PElement % BDofs = GetBubbleDofs( PElement, P )
687    PElement % PDefs % P = P
688    IF (ElementCode == 504) THEN
689      PElement % PDefs % TetraType = 1
690    ELSE
691      PElement % PDefs % TetraType = 0
692    END IF
693    PElement % PDefs % isEdge = .FALSE.
694    ! Set max element dofs to mesh
695    Mesh % MaxElementDOFs = PElement % Type % NumberOfNodes + &
696           PElement % Type % NumberOfEdges * Mesh % MaxEdgeDOFs + &
697           PElement % Type % NumberOfFaces * Mesh % MaxFaceDOFs + &
698           PElement % BDOFs
699
700    PElement % PDefs % GaussPoints = (P+1) ** PElement % Type % DIMENSION
701  END SUBROUTINE AllocateMeshAndPElement
702
703  SUBROUTINE InitializePElement(Element)
704    IMPLICIT NONE
705    TYPE(Element_t) :: Element
706
707    Element % BDOFs    =  0
708    Element % NDOFs    =  0
709    Element % BodyId   = 1
710    Element % Splitted =  0
711    Element % hK = 0
712    Element % ElementIndex = 0
713    Element % StabilizationMk = 0
714    NULLIFY( Element % TYPE )
715    NULLIFY( Element % PDefs )
716    NULLIFY( Element % BubbleIndexes )
717    NULLIFY( Element % DGIndexes )
718    NULLIFY( Element % NodeIndexes )
719    NULLIFY( Element % EdgeIndexes )
720    NULLIFY( Element % FaceIndexes )
721    NULLIFY( Element % BoundaryInfo )
722  END SUBROUTINE InitializePElement
723
724  FUNCTION ClonePElement(Element) RESULT(ClonedElement)
725    IMPLICIT NONE
726    TYPE(Element_t), POINTER :: Element, ClonedElement
727
728    ALLOCATE(ClonedElement)
729    CALL InitializePElement(ClonedElement)
730    ClonedElement % BDOFs = Element % BDOFs
731    ClonedElement % NDOFs =  Element % NDOFs
732    ClonedElement % BodyId = Element % BodyId
733    ClonedElement % Splitted = Element % Splitted
734    ClonedElement % hK = Element % hK
735    ClonedElement % ElementIndex = Element % ElementIndex
736    ClonedElement % StabilizationMk = Element % StabilizationMk
737    ClonedElement % Type => Element % Type
738    IF (ASSOCIATED(Element % PDefs)) THEN
739      CALL AllocatePDefinitions(ClonedElement)
740      ClonedElement % PDefs % P = Element % PDefs % P
741      ClonedElement % PDefs % TetraType = Element % PDefs % TetraType
742      ClonedElement % PDefs % isEdge = Element % PDefs % isEdge
743      ClonedElement % PDefs % GaussPoints = Element % PDefs % GaussPoints
744      ClonedElement % PDefs % pyramidQuadEdge = Element % PDefs % pyramidQuadEdge
745      ClonedElement % PDefs % localNumber = Element % PDefs % localNumber
746    END IF
747    IF (ASSOCIATED( Element % NodeIndexes )) THEN
748      ALLOCATE(ClonedElement % NodeIndexes( SIZE(Element % NodeIndexes)))
749      ClonedElement % NodeIndexes(:) = Element % NodeIndexes(:)
750    END IF
751    IF (ASSOCIATED( Element % EdgeIndexes )) THEN
752      ALLOCATE(ClonedElement % EdgeIndexes( SIZE(Element % EdgeIndexes)))
753      ClonedElement % EdgeIndexes(:) = Element % EdgeIndexes(:)
754    END IF
755    IF (ASSOCIATED( Element % FaceIndexes )) THEN
756      ALLOCATE(ClonedElement % FaceIndexes( SIZE(Element % FaceIndexes)))
757      ClonedElement % FaceIndexes(:) = Element % FaceIndexes(:)
758    END IF
759  END FUNCTION ClonePElement
760
761  SUBROUTINE DeallocatePElement(PElement)
762    IMPLICIT NONE
763
764    TYPE(Element_t), POINTER :: PElement
765
766    IF (ASSOCIATED(PElement % PDefs)) DEALLOCATE(PElement % PDefs)
767    IF (ASSOCIATED(PElement % NodeIndexes)) DEALLOCATE(PElement % NodeIndexes)
768    IF (ASSOCIATED(PElement % EdgeIndexes)) DEALLOCATE(PElement % EdgeIndexes)
769    IF (ASSOCIATED(PElement % FaceIndexes)) DEALLOCATE(PElement % FaceIndexes)
770    DEALLOCATE(PElement)
771  END SUBROUTINE DeallocatePElement
772
773  SUBROUTINE DeallocateTemporaryMesh(Mesh)
774    IMPLICIT NONE
775
776     TYPE(Mesh_t), POINTER :: Mesh
777     INTEGER :: edge, face
778
779     IF (ASSOCIATED(Mesh % Elements)) DEALLOCATE(Mesh % Elements)
780     DO edge=1,Mesh % NumberOfEdges
781       IF (ASSOCIATED(Mesh % Edges(edge) % PDefs)) DEALLOCATE(Mesh % Edges(edge) % PDefs)
782     END DO
783     IF (ASSOCIATED(Mesh % Edges)) DEALLOCATE(Mesh % Edges)
784     DO face=1,Mesh % NumberOfFaces
785       IF (ASSOCIATED(Mesh % Faces(face) % PDefs)) DEALLOCATE(Mesh % Faces(face) % PDefs)
786     END DO
787     IF (ASSOCIATED(Mesh % Faces)) DEALLOCATE(Mesh % Faces)
788
789     DEALLOCATE( Mesh % Nodes )
790     DEALLOCATE( Mesh )
791  END SUBROUTINE DeallocateTemporaryMesh
792
793  ! Portable wall-clock timer
794  FUNCTION ftimer() RESULT(timerval)
795    IMPLICIT NONE
796
797    REAL(KIND=dp) :: timerval
798    INTEGER(KIND=8) :: t, rate
799
800#ifdef _OPENMP
801    timerval = OMP_GET_WTIME()
802#else
803    CALL SYSTEM_CLOCK(t,count_rate=rate)
804    timerval = REAL(t,dp)/rate
805#endif
806  END FUNCTION ftimer
807
808!------------------------------------------------------------------------------
809END SUBROUTINE LinearFormsAssembly
810!------------------------------------------------------------------------------
811