CoCalc -- ElementUtils.F90

GitHub Repository: ElmerCSC/elmerfem
Path: blob/devel/fem/src/ElementUtils.F90
⁵²¹³ views
1
!/*****************************************************************************/
2
! *
3
! *  Elmer, A Finite Element Software for Multiphysical Problems
4
! *
5
! *  Copyright 1st April 1995 - , CSC - IT Center for Science Ltd., Finland
6
! * 
7
! *  This library is free software; you can redistribute it and/or
8
! *  modify it under the terms of the GNU Lesser General Public
9
! *  License as published by the Free Software Foundation; either
10
! *  version 2.1 of the License, or (at your option) any later version.
11
! *
12
! *  This library is distributed in the hope that it will be useful,
13
! *  but WITHOUT ANY WARRANTY; without even the implied warranty of
14
! *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
15
! *  Lesser General Public License for more details.
16
! * 
17
! *  You should have received a copy of the GNU Lesser General Public
18
! *  License along with this library (in file ../LGPL-2.1); if not, write 
19
! *  to the Free Software Foundation, Inc., 51 Franklin Street, 
20
! *  Fifth Floor, Boston, MA  02110-1301  USA
21
! *
22
! *****************************************************************************/
23
!
24
!/******************************************************************************
25
! *
26
! *  Authors: Juha Ruokolainen
27
! *  Email:   [email protected]
28
! *  Web:     http://www.csc.fi/elmer
29
! *  Address: CSC - IT Center for Science Ltd.
30
! *           Keilaranta 14
31
! *           02101 Espoo, Finland 
32
! *
33
! *  Original Date: 01 Oct 1996
34
! *
35
! ******************************************************************************/
36

37
!> \ingroup ElmerLib
38
!> \}
39

40
!--------------------------------------------------------------------------------
41
!>  Some basic finite element utilities.
42
!--------------------------------------------------------------------------------
43
#include "../config.h"
44

45
MODULE ElementUtils
46

47

48
    USE DirectSolve
49
    USE ListMatrixArray
50
    USE Integration
51
    USE Lists
52
    USE Interpolation
53
    USE BandwidthOptimize
54
    USE ElementDescription, ONLY : getEdgeDOFs,GetBubbleDOFs,getFaceDOFs, &
55
      CrossProduct, NormalVector, InterpolateInElement, mGetElementDOFs
56
            
57
    IMPLICIT NONE
58

59
CONTAINS
60

61
!------------------------------------------------------------------------------
62
!> Frees structures of the matrix.
63
!------------------------------------------------------------------------------
64
   RECURSIVE SUBROUTINE FreeMatrix(Matrix)
65
!------------------------------------------------------------------------------
66
     TYPE(Matrix_t), POINTER :: Matrix
67
!------------------------------------------------------------------------------
68
     TYPE(Solver_t) :: Solver
69
     REAL(KIND=dp) :: x(1), b(1)
70
     INTEGER :: i
71

72
     TYPE(SplittedMatrixT), POINTER :: s
73
     TYPE(BasicMatrix_t), POINTER :: m
74
     TYPE(SParIterSolverGlobalD_t), POINTER :: p
75

76
     LOGICAL :: Found
77

78

79
#ifdef HAVE_HYPRE
80
    INTERFACE
81
      !! destroy the data structures (should be called when the matrix has
82
      !! to be updated and SolveHYPRE1 has to be called again).
83
      SUBROUTINE SolveHYPRE4(hypreContainer, verbosity ) BIND(C,name="solvehypre4")
84
        USE, INTRINSIC :: iso_c_binding
85
        INTEGER(KIND=C_INTPTR_T) :: hypreContainer
86
        INTEGER(KIND=c_int) :: verbosity
87
      END SUBROUTINE SolveHYPRE4
88
    END INTERFACE
89
#endif
90

91
!------------------------------------------------------------------------------
92

93
     IF ( .NOT. ASSOCIATED( Matrix ) ) RETURN
94

95
     CALL DirectSolver( Matrix,x,b,Solver,Free_Fact=.TRUE.)
96

97
     IF ( ASSOCIATED( Matrix % Perm ) )        DEALLOCATE( Matrix % Perm )
98
     IF ( ASSOCIATED( Matrix % InvPerm ) )     DEALLOCATE( Matrix % InvPerm )
99

100
     IF ( ASSOCIATED( Matrix % Rows ) ) THEN
101
        IF ( ASSOCIATED(Matrix % Rows, Matrix % ILURows) .OR. SIZE(Matrix % Rows)==0 )  &
102
           Matrix % ILURows => Null()
103
        DEALLOCATE( Matrix % Rows )
104
     END IF
105

106
     IF ( ASSOCIATED( Matrix % Cols ) ) THEN
107
        IF ( ASSOCIATED(Matrix % Cols, Matrix % ILUCols) .OR. SIZE(Matrix % Cols)==0 ) &
108
           Matrix % ILUCols => Null()
109
        DEALLOCATE( Matrix % Cols )
110
     END IF
111

112
     IF ( ASSOCIATED( Matrix % Diag ) ) THEN
113
        IF ( ASSOCIATED(Matrix % Diag, Matrix % ILUDiag) .OR. SIZE(Matrix % Diag)==0 ) &
114
          Matrix % ILUDiag => Null()
115
        DEALLOCATE( Matrix % Diag )
116
     END IF
117

118
     IF ( ASSOCIATED( Matrix % RHS   ) )       DEALLOCATE( Matrix % RHS )
119
     IF ( ASSOCIATED( Matrix % Force ) )       DEALLOCATE( Matrix % Force )
120
     IF ( ASSOCIATED( Matrix % RHS_im ) )      DEALLOCATE( Matrix % RHS_im )
121

122
     IF ( ALLOCATED( Matrix % ExtraVals ) )    DEALLOCATE( Matrix % ExtraVals )
123
     IF ( ASSOCIATED( Matrix % Values ) )      DEALLOCATE( Matrix % Values )
124
     IF ( ASSOCIATED( Matrix % Values_im ) )   DEALLOCATE( Matrix % Values_im )
125
     IF ( ASSOCIATED( Matrix % MassValues ) )  DEALLOCATE( Matrix % MassValues )
126
     IF ( ASSOCIATED( Matrix % DampValues ) )  DEALLOCATE( Matrix % DampValues )
127
     IF ( ASSOCIATED( Matrix % PrecValues ) )  DEALLOCATE( Matrix % PrecValues )
128
     IF ( ASSOCIATED( Matrix % BulkValues ) )  DEALLOCATE( Matrix % BulkValues )
129
     IF ( ASSOCIATED( Matrix % BulkRHS   ) )   DEALLOCATE( Matrix % BulkRHS )
130

131
     IF ( ASSOCIATED( Matrix % ILUValues ) )   DEALLOCATE( Matrix % ILUValues )
132
     IF ( ASSOCIATED( Matrix % ILURows ) )     DEALLOCATE( Matrix % ILURows )
133
     IF ( ASSOCIATED( Matrix % ILUCols ) )     DEALLOCATE( Matrix % ILUCols )
134
     IF ( ASSOCIATED( Matrix % ILUDiag ) )     DEALLOCATE( Matrix % ILUDiag )
135

136
     IF ( ASSOCIATED( Matrix % CRHS   ) )      DEALLOCATE( Matrix % CRHS )
137
     IF ( ASSOCIATED( Matrix % CForce ) )      DEALLOCATE( Matrix % CForce )
138

139
     IF ( ASSOCIATED( Matrix % CValues ) )     DEALLOCATE( Matrix % CValues )
140
     IF ( ASSOCIATED( Matrix % CILUValues ) )  DEALLOCATE( Matrix % CILUValues )
141

142
     IF ( ASSOCIATED(Matrix % CMassValues) )  DEALLOCATE( Matrix % CMassValues )
143
     IF ( ASSOCIATED(Matrix % CDampValues) )  DEALLOCATE( Matrix % CDampValues )
144

145
     IF ( ALLOCATED( Matrix % GRows ) )      DEALLOCATE( Matrix % GRows )
146
     IF ( ALLOCATED( Matrix % RowOwner ) )   DEALLOCATE( Matrix % RowOwner )
147
     IF ( ASSOCIATED( Matrix % GOrder) )      DEALLOCATE( Matrix % GOrder )
148

149
     CALL FreeMatrix( Matrix % CircuitMatrix )
150
     CALL FreeMatrix( Matrix % AddMatrix )
151
     CALL FreeMatrix( Matrix % EMatrix )
152
     CALL FreeMatrix( Matrix % ConstraintMatrix )
153
     CALL FreeMatrix( Matrix % CollectionMatrix )
154

155
     IF  (ALLOCATED(Matrix % Dvalues) ) DEALLOCATE(Matrix % DValues)
156
     IF  (ALLOCATED(Matrix % ConstrainedDOF) ) DEALLOCATE(Matrix % ConstrainedDOF)
157
     IF  (ASSOCIATED(Matrix % EPerm) ) DEALLOCATE(Matrix % EPerm)
158
     IF  (ASSOCIATED(Matrix % BulkResidual) ) DEALLOCATE(Matrix % BulKResidual)
159
     IF  (ASSOCIATED(Matrix % DiagScaling) ) DEALLOCATE(Matrix % DiagScaling)
160
     IF  (ASSOCIATED(Matrix % Tvalues) ) DEALLOCATE(Matrix % Tvalues)
161
     IF  (ASSOCIATED(Matrix % BulkMassValues) ) DEALLOCATE(Matrix % BulkMassValues)
162
     IF  (ASSOCIATED(Matrix % BulkDampValues) ) DEALLOCATE(Matrix % BulkDampValues)
163
     IF  (ASSOCIATED(Matrix % HaloValues) ) DEALLOCATE(Matrix % HaloValues)
164
     IF  (ASSOCIATED(Matrix % HaloMassValues) ) DEALLOCATE(Matrix % HaloMassValues)
165

166
     IF(ASSOCIATED(Matrix % ParallelInfo)) THEN
167
       IF(ASSOCIATED(Matrix % ParallelInfo % GlobalDOFs)) DEALLOCATE(Matrix % ParallelInfo % GlobalDOFs)
168
       IF(ASSOCIATED(Matrix % ParallelInfo % GInterface)) DEALLOCATE(Matrix % ParallelInfo % GInterface)
169
       IF(ASSOCIATED(Matrix % ParallelInfo % Gorder)) DEALLOCATE(Matrix % ParallelInfo % GOrder)
170
 
171
       IF(ASSOCIATED(Matrix % ParallelInfo % NeighbourList)) THEN
172
         DO i=1,SIZE(Matrix % ParallelInfo % NeighbourList)
173
           IF (ASSOCIATED(Matrix % ParallelInfo % NeighbourList(i) % Neighbours)) &
174
             DEALLOCATE(Matrix % ParallelInfo % NeighbourList(i) % Neighbours)
175
         END DO
176
         DEALLOCATE(Matrix % ParallelInfo % NeighbourList)
177
       END IF
178

179
       IF(ASSOCIATED(Matrix % ParallelInfo % FaceNeighbourList)) THEN
180
         DO i=1,SIZE(Matrix % ParallelInfo % FaceNeighbourList)
181
           IF (ASSOCIATED(Matrix % ParallelInfo % FaceNeighbourList(i) % Neighbours)) &
182
             DEALLOCATE(Matrix % ParallelInfo % FaceNeighbourList(i) % Neighbours)
183
         END DO
184
         DEALLOCATE(Matrix % ParallelInfo % FaceNeighbourList)
185
       END IF
186
       IF(ASSOCIATED(Matrix % ParallelInfo % FaceInterface)) DEALLOCATE(Matrix % ParallelInfo % FaceInterface)
187

188
       IF(ASSOCIATED(Matrix % ParallelInfo % EdgeNeighbourList)) THEN
189
         DO i=1,SIZE(Matrix % ParallelInfo % EdgeNeighbourList)
190
           IF (ASSOCIATED(Matrix % ParallelInfo % EdgeNeighbourList(i) % Neighbours)) &
191
             DEALLOCATE(Matrix % ParallelInfo % EdgeNeighbourList(i) % Neighbours)
192
         END DO
193
         DEALLOCATE(Matrix % ParallelInfo % EdgeNeighbourList)
194
       END IF
195
       IF(ASSOCIATED(Matrix % ParallelInfo % EdgeInterface)) DEALLOCATE(Matrix % ParallelInfo % EdgeInterface)
196
       DEALLOCATE(Matrix % ParallelInfo)
197
     END IF
198
         
199

200
     p=>Matrix % ParMatrix
201
     IF(ASSOCIATED(p)) THEN
202
       s => Matrix % ParMatrix % SplittedMatrix
203

204
       s % InsideMatrix % EMatrix => NULL()
205
       CALL FreeMatrix(s % InsideMatrix)
206

207
       DO i=1,SIZE(s % IfMatrix)
208
         m => s % IfMatrix(i)
209
         IF(ASSOCIATED(m)) THEN
210
           IF(ALLOCATED(m % Rows)) DEALLOCATE(m % Rows)
211
           IF(ALLOCATED(m % Cols)) DEALLOCATE(m % Cols)
212
           IF(ALLOCATED(m % Diag)) DEALLOCATE(m % Diag)
213
           IF(ALLOCATED(m % Grows)) DEALLOCATE(m % Grows)
214
           IF(ALLOCATED(m % Values)) DEALLOCATE(m % Values)
215
           IF(ALLOCATED(m % RowOwner)) DEALLOCATE(m % RowOwner)
216
           IF(ALLOCATED(m % ILUValues)) DEALLOCATE(m % ILUValues)
217
           IF(ALLOCATED(m % MassValues)) DEALLOCATE(m % MassValues)
218
           IF(ALLOCATED(m % DampValues)) DEALLOCATE(m % DampValues)
219
           IF(ALLOCATED(m % PrecValues)) DEALLOCATE(m % PrecValues)
220
         END IF
221
       END DO
222
       DEALLOCATE(s % IfMatrix)
223

224
       DO i=1,SIZE(s % NbsIfMatrix)
225
         m => s % NbsIfMatrix(i)
226
         IF(ASSOCIATED(m)) THEN
227
           IF(ALLOCATED(m % Rows)) DEALLOCATE(m % Rows)
228
           IF(ALLOCATED(m % Cols)) DEALLOCATE(m % Cols)
229
           IF(ALLOCATED(m % Diag)) DEALLOCATE(m % Diag)
230
           IF(ALLOCATED(m % Grows)) DEALLOCATE(m % Grows)
231
           IF(ALLOCATED(m % Values)) DEALLOCATE(m % Values)
232
           IF(ALLOCATED(m % RowOwner)) DEALLOCATE(m % RowOwner)
233
           IF(ALLOCATED(m % ILUValues)) DEALLOCATE(m % ILUValues)
234
           IF(ALLOCATED(m % MassValues)) DEALLOCATE(m % MassValues)
235
           IF(ALLOCATED(m % DampValues)) DEALLOCATE(m % DampValues)
236
           IF(ALLOCATED(m % PrecValues)) DEALLOCATE(m % PrecValues)
237
         END IF
238
       END DO
239
       DEALLOCATE(s % NbsIfMatrix)
240

241
       IF(ASSOCIATED(s % VecIndices)) THEN
242
         DO i=1,SIZE(s % VecIndices)
243
           IF(ASSOCIATED(s % Vecindices(i) % RevInd)) DEALLOCATE(s % VecIndices(i) % RevInd)
244
         END DO
245
         DEALLOCATE(s % VecIndices)
246
       END IF
247

248
       IF(ASSOCIATED(s % IfVecs)) THEN
249
         DO i=1,SIZE(s % IfVecs)
250
           IF(ASSOCIATED(s % IfVecs(i) % IfVec)) DEALLOCATE(s % IfVecs(i) % IfVec)
251
         END DO
252
         DEALLOCATE(s % ifVecs)
253
       END IF
254

255
       IF(ASSOCIATED(s % IfORows)) THEN
256
         DO i=1,SIZE(s % IfORows)
257
           IF(ASSOCIATED(s % IfORows(i) % IfVec)) DEALLOCATE(s % IfORows(i) % IfVec)
258
         END DO
259
         DEALLOCATE(s % ifORows)
260
       END IF
261

262
       IF(ASSOCIATED(s % IfLCols)) THEN
263
         DO i=1,SIZE(s % IfLCols)
264
           IF(ASSOCIATED(s % IfLCols(i) % IfVec)) DEALLOCATE(s % IfLCols(i) % IfVec)
265
         END DO
266
         DEALLOCATE(s % ifLCols)
267
       END IF
268

269
       IF(ASSOCIATED(s % ResBuf)) THEN
270
         DO i=1,SIZE(s % ResBuf)
271
           IF(ALLOCATED(s % ResBuf(i) % ResVal)) DEALLOCATE(s % ResBuf(i) % ResVal)
272
           IF(ALLOCATED(s % ResBuf(i) % ResInd)) DEALLOCATE(s % ResBuf(i) % ResInd)
273
         END DO
274
         DEALLOCATE(s % ResBuf)
275
       END IF
276

277
       IF(ASSOCIATED(s % RHS)) THEN
278
         DO i=1,SIZE(s % RHS)
279
           IF(ASSOCIATED(s % RHS(i) % RHSVec)) DEALLOCATE(s % RHS(i) % RHSVec)
280
           IF(ASSOCIATED(s % RHS(i) % RHSInd)) DEALLOCATE(s % RHS(i) % RHSInd)
281
         END DO
282
         DEALLOCATE(s % RHS)
283
       END IF
284

285
       IF (ASSOCIATED(s % Work)) DEALLOCATE(s % Work)
286
       IF (ASSOCIATED(s % TmpXVec)) DEALLOCATE(s % TmpXVec)
287
       IF (ASSOCIATED(s % TmpRVec)) DEALLOCATE(s % TmpRVec)
288

289
       IF(ASSOCIATED(s % GlueTable)) THEN
290
         IF(ASSOCIATED(s % GlueTable % Rows)) DEALLOCATE(s % GlueTable % Rows)
291
         IF(ASSOCIATED(s % GlueTable % Cols)) DEALLOCATE(s % GlueTable % Cols)
292
         IF(ASSOCIATED(s % GlueTable % Inds)) DEALLOCATE(s % GlueTable % Inds)
293
         IF(ASSOCIATED(s % GlueTable % RowOwner)) DEALLOCATE(s % GlueTable % RowOwner)
294
         DEALLOCATE(s % GlueTable)
295
       END IF
296
       DEALLOCATE(s)
297

298
       IF(ASSOCIATED(p % ParEnv % Active)) THEN
299
         DEALLOCATE(p % ParEnv % Active)
300
       END IF
301

302
       IF(ASSOCIATED(p % ParEnv % Isneighbour)) THEN
303
         DEALLOCATE(p % ParEnv % Isneighbour)
304
       END IF
305

306
       DEALLOCATE(p)
307
     END IF
308

309
#ifdef HAVE_HYPRE
310
     IF (Matrix % Hypre /= 0) THEN
311
       i = ListGetInteger( CurrentModel % Simulation,'Max Output Level',Found ) 
312
       IF(ParEnv % MyPe /= 0) i = 0 
313
       CALL SolveHypre4(Matrix % Hypre, i )
314
       Matrix % Hypre = 0
315
     END IF
316
#endif
317

318
     DEALLOCATE( Matrix )
319
!------------------------------------------------------------------------------
320
   END SUBROUTINE FreeMatrix
321
!------------------------------------------------------------------------------
322

323

324
   
325
!------------------------------------------------------------------------------
326
!> Create a list matrix given the mesh, the active domains and the elementtype 
327
!> related to the solver. The list matrix is flexible since it can account 
328
!> for any entries. Also constraints and periodic BCs may give rise to entries
329
!> in the list matrix topology.
330
!------------------------------------------------------------------------------
331
  SUBROUTINE MakeListMatrix( Model,Solver,Mesh,List,Reorder, &
332
        LocalNodes, Equation, DGSolver, GlobalBubbles, &
333
        NodalDofsOnly, ProjectorDofs, CalcNonZeros )
334
!------------------------------------------------------------------------------
335
    TYPE(Model_t)  :: Model
336
    TYPE(Solver_t) :: Solver
337
    TYPE(Mesh_t)   :: Mesh
338
    TYPE(ListMatrix_t), POINTER :: List(:)
339
    INTEGER, OPTIONAL :: Reorder(:)
340
    INTEGER :: LocalNodes
341
    CHARACTER(LEN=*), OPTIONAL :: Equation
342
    LOGICAL, OPTIONAL :: DGSolver
343
    LOGICAL, OPTIONAL :: GlobalBubbles
344
    LOGICAL, OPTIONAL :: NodalDofsOnly
345
    LOGICAL, OPTIONAL :: ProjectorDofs
346
    LOGICAL, OPTIONAL :: CalcNonZeros
347
!------------------------------------------------------------------------------
348
    INTEGER :: t,i,j,k,l,m,k1,k2,n,p,q,e1,e2,f1,f2, EDOFs, FDOFs, BDOFs, This, istat
349
    INTEGER :: NDOFs, DOFsPerNode
350
    INTEGER, ALLOCATABLE :: InvPerm(:), IndirectPairs(:)
351
    LOGICAL :: Flag, FoundDG, GB, DB, Found, Radiation, DoProjectors, &
352
        DoNonZeros, DgIndirect
353
    TYPE(Matrix_t), POINTER :: Projector
354
    INTEGER :: IndexSize, NumberOfFactors
355
    INTEGER, ALLOCATABLE :: Indexes(:)
356
    TYPE(ListMatrixEntry_t), POINTER :: CList, Lptr
357
    TYPE(Matrix_t),POINTER :: PMatrix
358
    TYPE(Element_t), POINTER :: Element,Elm, Edge1, Edge2, Face1, Face2, Left, Right
359
    CHARACTER(:), ALLOCATABLE :: RadiationFlag
360
    LOGICAL :: GotIt, PSA
361
    CHARACTER(*), PARAMETER :: Caller = 'MakeListMatrix'
362
!------------------------------------------------------------------------------
363

364
    CALL Info(Caller,'Creating list matrix',Level=14)
365

366
    GB = .FALSE.
367
    IF ( PRESENT(GlobalBubbles) ) GB = GlobalBubbles
368

369
    IF( DgSolver ) THEN
370
      DB = ListGetLogical( Solver % Values,'DG Reduced Basis',Found) 
371
    ELSE
372
      DB = .FALSE.
373
    END IF
374

375
    ! When this has been checked properly the old can be removed
376
    PSA = ListGetLogical( Solver % Values,'PSA',Found ) 
377
    IF(.NOT. Found) PSA = .TRUE.
378
    
379
    List => List_AllocateMatrix(LocalNodes)
380

381
    NDOFs = Mesh % MaxNDOFs
382
    BDOFs = Mesh % MaxBDOFs
383
    EDOFs = Mesh % MaxEdgeDOFs
384
    FDOFs = Mesh % MaxFaceDOFs
385
    IF( PRESENT( NodalDofsOnly ) ) THEN
386
      IF( NodalDofsOnly ) THEN
387
        EDOFS = 0
388
        FDOFS = 0
389
      END IF
390
    END IF
391

392
    IF( PRESENT( ProjectorDofs ) ) THEN
393
      DoProjectors = ProjectorDofs 
394
    ELSE
395
      DoProjectors = .TRUE.
396
    END IF
397

398
    IF( EDOFS > 0 .AND. .NOT. ASSOCIATED(Mesh % Edges) ) THEN
399
      CALL Warn(Caller,'Edge dofs requested but not edges exist in mesh!')
400
      EDOFS = 0
401
    END IF
402

403
    IF( FDOFS > 0 .AND. .NOT. ASSOCIATED(Mesh % Faces) ) THEN
404
      CALL Warn(Caller,'Face dofs requested but not faces exist in mesh!')
405
      FDOFS = 0
406
    END IF
407

408
    IndexSize = 128
409
    ALLOCATE( Indexes(IndexSize), STAT=istat )
410
    IF( istat /= 0 ) THEN
411
      CALL Fatal(Caller,'Allocation error for Indexes')
412
    END IF
413

414
    ! Create sparse matrix for the Discontinuous Galerkin solver 
415
    ! Using either reduced or full basis
416
    !-------------------------------------------------------------------
417
    FoundDG = .FALSE.
418
    IF ( DGSolver ) THEN    
419
      ! Create the sparse matrix for the Discontinuous Bodies solver
420
      !-------------------------------------------------------------------
421
      IF ( DB ) THEN
422
        DGIndirect = ListGetLogical( Solver % Values,'DG Indirect Connections',Found )
423

424
        IF( DGIndirect ) THEN
425
          CALL Info(Caller,'Creating also indirect connections!',Level=12)
426
          ALLOCATE( IndirectPairs( LocalNodes ), STAT=istat )
427
          IF ( istat /= 0 ) CALL Fatal(Caller,'Memory allocation error for IndirectPairs work.')
428
          IndirectPairs = 0
429
        END IF
430
                        
431

432
        DO t=1,Mesh % NumberOfBulkElements
433
          n = 0
434
          Elm => Mesh % Elements(t)
435
          IF ( .NOT. CheckElementEquation(Model,Elm,Equation) ) CYCLE
436

437
          FoundDG = FoundDG .OR. Elm % DGDOFs > 0
438
          DO j=1,Elm % DGDOFs
439
            n = n + 1
440
            Indexes(n) = Elm % DGIndexes(j)
441
          END DO
442

443
          IF( Elm % DGDofs /= Elm % TYPE % NumberOfNodes ) THEN
444
            CALL Fatal(Caller,'Mismatch in sizes in reduced basis DG!')
445
          END IF
446

447
          IF( PSA ) THEN
448
            CALL PackSortAdd(n,Indexes,Reorder)
449
          ELSE
450
            DO i=1,n
451
              k1 = Reorder(Indexes(i))
452
              IF ( k1 <= 0 ) CYCLE
453
              DO j=1,n
454
                k2 = Reorder(Indexes(j))
455
                IF ( k2 <= 0 ) CYCLE
456
                Lptr => List_GetMatrixIndex( List,k1,k2 )
457
              END DO
458
            END DO
459
          END IF
460
        END DO
461

462
        IF( Mesh % NumberOfFaces == 0 ) THEN
463
          DO t=1,Mesh % NumberOfEdges
464
            n = 0
465
            Elm => Mesh % Edges(t)
466
            IF(.NOT. ASSOCIATED( Elm % BoundaryInfo ) ) CYCLE
467

468
            Left => Elm % BoundaryInfo % Left
469
            IF(.NOT. ASSOCIATED( Left ) ) CYCLE
470
            IF ( .NOT. CheckElementEquation(Model,Left,Equation) ) CYCLE
471

472
            Right => Elm % BoundaryInfo %  Right
473
            IF(.NOT. ASSOCIATED( Right ) ) CYCLE
474
            IF ( .NOT. CheckElementEquation(Model,Right,Equation) ) CYCLE
475

476
            IF( Left % BodyId == Right % BodyId ) CYCLE
477

478
            IF( DGIndirect ) THEN
479
              DO i=1,Left % DGDOFs
480
                k1 = ReOrder( Left % DgIndexes(i) )
481
                DO j=1,Right % DGDOFs
482
                  IF( Left % NodeIndexes(i) == Right % NodeIndexes(j) ) THEN
483
                    k2 = ReOrder( Right % DgIndexes(j) )
484
                    IF( k1 /= k2 ) THEN
485
                      IndirectPairs( k1 ) = k2
486
                      EXIT
487
                    END IF
488
                  END IF
489
                END DO
490
              END DO
491
            END IF
492

493
            FoundDG = FoundDG .OR. Left % DGDOFs > 0
494
            DO j=1,Left % DGDOFs
495
              n = n + 1
496
              Indexes(n) = Left % DGIndexes(j)
497
            END DO
498

499
            FoundDG = FoundDG .OR. Right % DGDOFs > 0
500
            DO j=1,Right % DGDOFs
501
              n = n + 1
502
              Indexes(n) = Right % DGIndexes(j)
503
            END DO
504

505
            IF( PSA ) THEN
506
              CALL PackSortAdd(n,Indexes,Reorder)
507
            ELSE
508
              DO i=1,n
509
                k1 = Reorder(Indexes(i))
510
                IF ( k1 <= 0 ) CYCLE
511
                DO j=1,n
512
                  k2 = Reorder(Indexes(j))
513
                  IF ( k2 <= 0 ) CYCLE
514
                  Lptr => List_GetMatrixIndex( List,k1,k2 )
515
                END DO
516
              END DO
517
            END IF
518
          END DO
519
        END IF
520

521
        
522
        DO t=1,Mesh % NumberOfFaces
523
          n = 0
524

525
          Elm => Mesh % Faces(t)
526
          IF(.NOT. ASSOCIATED( Elm % BoundaryInfo ) ) CYCLE
527
          
528
          Left => Elm  % BoundaryInfo % Left
529
          IF(.NOT. ASSOCIATED( Left ) ) CYCLE
530
          IF ( .NOT. CheckElementEquation(Model,Left,Equation) ) CYCLE
531

532
          Right => Elm % BoundaryInfo %  Right
533
          IF(.NOT. ASSOCIATED( Right ) ) CYCLE
534
          IF ( .NOT. CheckElementEquation(Model,Right,Equation) ) CYCLE
535

536
          IF( Left % BodyId == Right % BodyId ) CYCLE
537

538
          IF( DGIndirect ) THEN
539
            DO i=1,Left % DGDOFs
540
              k1 = ReOrder( Left % DgIndexes(i) )
541
              DO j=1,Right % DGDOFs
542
                IF( Left % NodeIndexes(i) == Right % NodeIndexes(j) ) THEN
543
                  k2 = ReOrder( Right % DgIndexes(j) )
544
                  IF( k1 /= k2 ) THEN
545
                    IF( IndirectPairs(k1) > 0 .AND. IndirectPairs(k1) /= k2 ) THEN
546
                      PRINT *,'Problematic node:',k1,IndirectPairs(k1),k2
547
                    END IF
548
                    IndirectPairs( k1 ) = k2
549
                    EXIT
550
                  END IF
551
                END IF
552
              END DO
553
            END DO
554
          END IF            
555

556
          FoundDG = FoundDG .OR. Left % DGDOFs > 0
557
          DO j=1,Left % DGDOFs
558
            n = n + 1
559
            Indexes(n) = Left % DGIndexes(j)
560
          END DO
561

562
          FoundDG = FoundDG .OR. Right % DGDOFs > 0
563
          DO j=1,Right % DGDOFs
564
            n = n + 1
565
            Indexes(n) = Right % DGIndexes(j)
566
          END DO
567

568
          IF( PSA ) THEN
569
            CALL PackSortAdd(n,Indexes,Reorder)
570
          ELSE
571
            DO i=1,n
572
              k1 = Reorder(Indexes(i))
573
              IF ( k1 <= 0 ) CYCLE
574
              DO j=1,n
575
                k2 = Reorder(Indexes(j))
576
                IF ( k2 <= 0 ) CYCLE
577
                Lptr => List_GetMatrixIndex( List,k1,k2 )
578
              END DO
579
            END DO
580
          END IF
581
        END DO
582

583
        IF( DGIndirect ) THEN
584
          DO k1 = 1, LocalNodes
585
            k2 = IndirectPairs(k1)
586
            IF( k2 == 0 ) CYCLE 
587
            !PRINT *,'Exchange structure between rows:',k1,k2
588
            CALL List_ExchangeRowStructure( List, k1, k2 ) 
589
          END DO
590
        END IF
591
          
592
        
593
      ELSE
594
        ! Classical DG solver
595
        !-------------------------------------        
596
        DO t=1,Mesh % NumberOfEdges
597
          n = 0
598
          Elm => Mesh % Edges(t) % BoundaryInfo % Left
599
          IF ( ASSOCIATED( Elm ) ) THEN
600
            IF ( CheckElementEquation(Model,Elm,Equation) ) THEN
601
              FoundDG = FoundDG .OR. Elm % DGDOFs > 0
602
              DO j=1,Elm % DGDOFs
603
                n = n + 1
604
                Indexes(n) = Elm % DGIndexes(j)
605
              END DO
606
            END IF
607
          END IF
608

609
          Elm => Mesh % Edges(t) % BoundaryInfo %  Right
610
          IF ( ASSOCIATED( Elm ) ) THEN
611
            IF ( CheckElementEquation(Model,Elm,Equation) ) THEN
612
              FoundDG = FoundDG .OR. Elm % DGDOFs > 0
613
              DO j=1,Elm % DGDOFs
614
                n = n + 1
615
                Indexes(n) = Elm % DGIndexes(j)
616
              END DO
617
            END IF
618
          END IF
619

620
          IF( PSA ) THEN
621
            CALL PackSortAdd(n,Indexes,Reorder)
622
          ELSE
623
            DO i=1,n
624
              k1 = Reorder(Indexes(i))
625
              IF ( k1 <= 0 ) CYCLE
626
              DO j=1,n
627
                k2 = Reorder(Indexes(j))
628
                IF ( k2 <= 0 ) CYCLE
629
                Lptr => List_GetMatrixIndex( List,k1,k2 )
630
              END DO
631
            END DO
632
          END IF          
633
        END DO
634
        
635
        DO t=1,Mesh % NumberOfFaces
636
          n = 0
637
          Elm => Mesh % Faces(t) % BoundaryInfo % Left
638
          IF ( ASSOCIATED( Elm ) ) THEN
639
            IF ( CheckElementEquation(Model,Elm,Equation) ) THEN
640
              FoundDG = FoundDG .OR. Elm % DGDOFs > 0
641
              DO j=1,Elm % DGDOFs
642
                n = n + 1
643
                Indexes(n) = Elm % DGIndexes(j)
644
              END DO
645
            END IF
646
          END IF
647

648
          Elm => Mesh % Faces(t) % BoundaryInfo %  Right
649
          IF ( ASSOCIATED( Elm ) ) THEN
650
            IF ( CheckElementEquation(Model,Elm,Equation) ) THEN
651
              FoundDG = FoundDG .OR. Elm % DGDOFs > 0
652
              DO j=1,Elm % DGDOFs
653
                n = n + 1
654
                Indexes(n) = Elm % DGIndexes(j)
655
              END DO
656
            END IF
657
          END IF
658

659
          IF( PSA ) THEN
660
            CALL PackSortAdd(n,Indexes,Reorder)
661
          ELSE
662
            DO i=1,n
663
              k1 = Reorder(Indexes(i))
664
              IF ( k1 <= 0 ) CYCLE
665
              DO j=1,n
666
                k2 = Reorder(Indexes(j))
667
                IF ( k2 <= 0 ) CYCLE
668
                Lptr => List_GetMatrixIndex( List,k1,k2 )
669
              END DO
670
            END DO
671
          END IF            
672
        END DO
673
      END IF
674
    END IF ! DGSolver
675

676
    
677
    !     Diffuse gray radiation condition:
678
    !     ---------------------------------
679
    Radiation = ListGetLogical( Solver % Values, 'Radiation Solver', Found )
680
    IF( Radiation ) THEN        
681
      Radiation = .FALSE.
682
      DO i=1,Model % NumberOfBCs
683
        RadiationFlag = ListGetString( Model % BCs(i) % Values, 'Radiation', GotIt )
684
        IF (GotIt) THEN
685
          IF ( RadiationFlag == 'diffuse gray' ) THEN
686
            Radiation = .TRUE.
687
            EXIT
688
          END IF
689
        END IF
690
      END DO
691
    END IF
692
    
693
    IF ( Radiation ) THEN
694
      BLOCK
695
        INTEGER, ALLOCATABLE :: Inds(:),ElemInds(:),ElemInds2(:)
696
        INTEGER :: cnt, maxnodes
697

698
        n = Mesh % MaxElementNodes
699
        ALLOCATE( ElemInds(n), ElemInds2(n), STAT=istat )
700
        IF ( istat /= 0 ) CALL Fatal(Caller,'Memory allocation error for ElemInds.')
701
        
702
        ElemInds = 0
703
        ElemInds2 = 0
704

705
        ! Check the maximum number of nodes in boundary element
706
        maxnodes = 0
707
        DO i = Mesh % NumberOfBulkElements+1, &
708
            Mesh % NumberOfBulkElements + Mesh % NumberOfBoundaryElements
709
          
710
          Element => Mesh % Elements(i)
711

712
          IF( .NOT. ASSOCIATED( Element % BoundaryInfo ) ) CYCLE
713
          IF ( .NOT. ASSOCIATED(Element % BoundaryInfo % RadiationFactors) ) CYCLE
714
          
715
          n = Element % Type % NumberOfNodes
716
          maxnodes = MAX( n, maxnodes )
717
        END DO
718
        
719
        
720
        CALL Info(Caller,'Adding radiation matrix',Level=14)      
721
        DO i = Mesh % NumberOfBulkElements+1, &
722
            Mesh % NumberOfBulkElements + Mesh % NumberOfBoundaryElements
723
          
724
          Element => Mesh % Elements(i)
725

726
          IF( .NOT. ASSOCIATED( Element % BoundaryInfo ) ) CYCLE
727
          IF ( .NOT. ASSOCIATED(Element % BoundaryInfo % RadiationFactors) ) CYCLE
728

729
          NumberOfFactors = Element % BoundaryInfo % &
730
              RadiationFactors % NumberOfImplicitFactors
731
          IF(NumberOfFactors == 0 ) CYCLE
732
          
733
          n = Element % Type % NumberOfNodes
734
          
735
          IF( DgSolver ) THEN
736
            CALL DgRadiationIndexes(Element,n,ElemInds,.TRUE.)
737
          END IF
738
            
739
          DO j=1,Element % TYPE % NumberOfNodes
740

741
            IF( DgSolver ) THEN
742
              k1 = Reorder(ElemInds(j))
743
            ELSE
744
              k1 = Reorder(Element % NodeIndexes(j))
745
            END IF
746
                        
747
            IF (.NOT.ALLOCATED(inds)) THEN
748
              ALLOCATE(inds(maxnodes*NumberOfFactors),STAT=istat)
749
            ELSE IF(SIZE(inds)<maxnodes*NumberOfFactors) THEN
750
              DEALLOCATE(inds)
751
              ALLOCATE(inds(maxnodes*NumberOfFactors),STAT=istat)
752
            END IF
753
            IF ( istat /= 0 ) CALL Fatal(Caller,'Memory allocation error for inds.')
754

755
            cnt = 0
756
            DO n=1,NumberOfFactors
757

758
              Elm => Mesh % Elements( Element % BoundaryInfo % &
759
                  RadiationFactors % Elements(n) )
760

761
              IF( DgSolver ) THEN
762
                CALL DgRadiationIndexes(Elm,Elm % TYPE % NumberOfNodes,ElemInds2,.TRUE.)
763
              END IF
764
              
765
              DO k=1,Elm % TYPE % NumberOfNodes
766
                cnt = cnt + 1
767

768
                IF( DGSolver ) THEN
769
                  inds(cnt) = Reorder(ElemInds2(k))
770
                ELSE
771
                  inds(cnt) = Reorder(Elm % NodeIndexes(k))
772
                END IF
773
              END DO
774
            END DO
775
            
776
            CALL Sort(cnt,inds)
777

778
            CALL List_AddMatrixIndexes(List,k1,cnt,inds)
779
          END DO
780
        END DO
781
      END BLOCK
782
      CALL Info(Caller,'Done Adding radiation matrix',Level=14)
783
    END IF
784

785
    
786
    ! If this is not a DG solver then create permutation considering 
787
    ! nodal, edge, face and bubble dofs. 
788
    !-------------------------------------------------------------------
789
    IF ( .NOT. FoundDG ) THEN
790
      t = 1
791
      DO WHILE( t<=Mesh % NumberOfBulkElements+Mesh % NumberOFBoundaryElements )
792
         Element => Mesh % Elements(t)
793

794
         IF ( PRESENT(Equation) ) THEN
795
           DO WHILE( t<=Mesh % NumberOfBulkElements+Mesh % NumberOfBoundaryElements )
796
             Element => Mesh % Elements(t)
797
             IF ( CheckElementEquation(Model,Element,Equation) ) EXIT
798
             t = t + 1
799
           END DO
800
           IF ( t > Mesh % NumberOfBulkElements+Mesh % NumberOfBoundaryElements ) EXIT
801
         END IF
802

803
         n = Element % NDOFs 
804
         IF( EDOFS > 0 ) n = n + Element % TYPE % NumberOfEdges * EDOFs 
805
         IF( FDOFS > 0 ) n = n + Element % TYPE % NumberOfFaces * FDOFs
806
         IF ( GB ) n = n + Element % BDOFs
807

808
         IF ( n > IndexSize ) THEN
809
            IndexSize = n
810
            IF ( ALLOCATED( Indexes ) ) DEALLOCATE( Indexes )
811
            ALLOCATE( Indexes(n), STAT=istat )
812
            IF( istat /= 0 ) CALL Fatal(Caller,'Allocation error for Indexes of size: '//I2S(n))
813
         END IF
814

815
         n = 0
816
         DOFsPerNode = Element % NDOFs / Element % TYPE % NumberOfNodes
817
         IF (DOFsPerNode > 0) THEN
818
           DO i=1,Element % TYPE % NumberOfNodes
819
             DO j=1,DOFsPerNode
820
               n = n + 1
821
               Indexes(n) = NDOFs * (Element % NodeIndexes(i)-1) + j
822
             END DO
823
           END DO
824
         END IF
825
!         DO i=1,Element % NDOFs
826
!            n = n + 1
827
!            Indexes(n) = Element % NodeIndexes(i)
828
!         END DO
829

830
         IF ( EDOFs > 0 ) THEN
831
            IF ( ASSOCIATED(Element % EdgeIndexes) ) THEN
832
              DO j=1,Element % TYPE % NumberOFEdges
833
                 DO i=1, Mesh % Edges(Element % EdgeIndexes(j)) % BDOFs
834
                   n = n + 1
835
                   Indexes(n) = EDOFs * (Element % EdgeIndexes(j)-1) + i &
836
                                + NDOFs * Mesh % NumberOfNodes
837
                END DO
838
             END DO
839

840
             IF ( GB ) THEN
841
               Edge1 => Mesh % Edges(Element % EdgeIndexes(1))
842
               IF(Element % Type % ElementCode==Edge1 % Type % ElementCode) THEN
843
                 DO i=1, Element % BDOFs
844
                   n = n + 1
845
                   Indexes(n) = EDOFs*(Element % EdgeIndexes(1)-1) + i + &
846
                       NDOFs * Mesh % NumberOfNodes
847
                 END DO
848
               END IF
849
             END IF
850
           END IF
851
         END IF
852

853
         IF ( FDOFS > 0 ) THEN
854
           IF ( ASSOCIATED(Element % FaceIndexes) ) THEN
855
             DO j=1,Element % TYPE % NumberOFFaces
856
               DO i=1, Mesh % Faces(Element % FaceIndexes(j)) % BDOFs
857
                 n = n + 1
858
                 Indexes(n) = FDOFs*(Element % FaceIndexes(j)-1) + i + &
859
                     NDOFs * Mesh % NumberOfNodes + EDOFs*Mesh % NumberOfEdges
860
               END DO
861
             END DO
862

863
             IF ( GB ) THEN
864
               Face1 => Mesh % Faces(Element % FaceIndexes(1))
865
               IF(Element % Type % ElementCode==Face1 % Type % ElementCode) THEN
866
                 DO i=1, Element % BDOFs
867
                   n = n + 1
868
                   Indexes(n) = FDOFs*(Element % FaceIndexes(1)-1) + i + &
869
                       NDOFs * Mesh % NumberOfNodes + EDOFs*Mesh % NumberOfEdges
870
                 END DO
871
               END IF
872
             END IF
873
           END IF
874
         END IF
875

876
         IF ( GB .AND. ASSOCIATED(Element % BubbleIndexes) ) THEN
877
            DO i=1,Element % BDOFs
878
              n = n + 1
879
              Indexes(n) = FDOFs*Mesh % NumberOfFaces + &
880
                   NDOFs * Mesh % NumberOfNodes + EDOFs*Mesh % NumberOfEdges + &
881
                        Element % BubbleIndexes(i)
882
            END DO
883
         END IF
884

885
         IF( PSA ) THEN
886
           CALL PackSortAdd(n,Indexes,Reorder)
887
         ELSE
888
           DO i=1,n
889
             k1 = Reorder(Indexes(i))
890
             IF ( k1 <= 0 ) CYCLE
891
             DO j=1,n
892
               k2 =  Reorder(Indexes(j))
893
               IF ( k2 <= 0 ) CYCLE
894
               Lptr => List_GetMatrixIndex( List,k1,k2 )
895
             END DO
896
           END DO
897
         END IF
898
         t = t + 1
899
      END DO
900

901
      IF ( ALLOCATED( Indexes ) ) DEALLOCATE( Indexes )
902

903
      
904
      DO i=Mesh % NumberOfBulkElements+1, Mesh % NumberOfBulkElements+ &
905
                     Mesh % NumberOfBoundaryElements
906
        IF ( Mesh % Elements(i) % TYPE % ElementCode <  102 .OR. &
907
             Mesh % Elements(i) % TYPE % ElementCode >= 200 ) CYCLE
908

909
        k1 = Reorder( Mesh % Elements(i) % NodeIndexes(1) )
910
        IF ( k1 > 0 ) THEN
911
          DO k=1,Mesh % Elements(i) % TYPE % NumberOFNodes
912
            k2 = Reorder( Mesh % Elements(i) % NodeIndexes(k) )
913
            IF ( k2 > 0 ) THEN
914
              Lptr => List_GetMatrixIndex( List,k1,k2 )
915
              Lptr => List_GetMatrixIndex( List,k2,k1 )
916
            END IF
917
          END DO
918
        END IF
919

920
        ! This is a connection element, make a matrix connection for that 
921
        IF ( Mesh % Elements(i) % TYPE % ElementCode == 102 ) THEN
922
          k2 = Reorder( Mesh % Elements(i) % NodeIndexes(2) )
923
          IF ( k2 > 0 ) Lptr => List_GetMatrixIndex( List,k2,k2 )
924
        END IF
925
      END DO
926

927

928
      ! Add connection from projectors. These are only needed if the projector is treated
929
      ! implicely. For explicit projectors or when using Lagrange coefficients the 
930
      ! connections are not needed. 
931
      !----------------------------------------------------------------------------------
932
      IF( DoProjectors ) THEN
933
        DO This=1,Model % NumberOfBCs
934
          Projector => Model % BCs(This) % PMatrix
935
          IF ( .NOT. ASSOCIATED(Projector) ) CYCLE
936

937
          IF( ListGetLogical( Model % BCs(This) % Values,&
938
              'Periodic BC Explicit',Found)) CYCLE
939
          IF( ListGetLogical( Model % BCs(This) % Values,&
940
              'Periodic BC Use Lagrange Coefficient',Found)) CYCLE
941

942
          CALL Info(Caller,'Adding matrix topology for BC: '//I2S(This),Level=10)
943

944
          DO i=1,Projector % NumberOfRows
945
            k = Reorder( Projector % InvPerm(i) )
946
            IF ( k > 0 ) THEN
947
              DO l=Projector % Rows(i),Projector % Rows(i+1)-1
948
                IF ( Projector % Cols(l) <= 0 ) CYCLE
949
                m = Reorder( Projector % Cols(l) )
950
                IF ( m > 0 ) THEN
951
                  Lptr => List_GetMatrixIndex( List,k,m )
952
                  Lptr => List_GetMatrixIndex( List,m,k ) ! keep structure symm.
953
                  CList => List(k) % Head
954
                  DO WHILE( ASSOCIATED( CList ) )
955
                    Lptr => List_GetMatrixIndex( List,m,CList % Index )
956
                    Lptr => List_GetMatrixIndex( List,CList % Index,m ) ! keep structure symm.
957
                    CList => CList % Next
958
                  END DO
959
                END IF
960
              END DO
961
            END IF
962
          END DO
963
        END DO
964
      END IF ! DoProjectors
965

966
    END IF
967
    
968
    Model % TotalMatrixElements = 0
969
    DO i=1,LocalNodes
970
      j = List(i) % Degree
971
      Model % TotalMatrixElements = Model % TotalMatrixElements + j
972
    END DO
973
    
974
      
975
    DoNonZeros = .TRUE.
976
    IF( PRESENT( CalcNonZeros ) ) DoNonZeros = CalcNonZeros
977
    
978
    IF( DoNonZeros ) THEN
979
      ALLOCATE( InvPerm(LocalNodes), STAT=istat )
980
      IF ( istat /= 0 ) CALL Fatal(Caller,'Memory allocation error for InvPerm.')
981
      
982
      InvPerm = 0
983
      k = 0
984
      DO i=1,SIZE(Reorder)
985
        IF (Reorder(i)>0) THEN
986
          k = k + 1
987
          InvPerm(Reorder(i)) = k
988
        END IF
989
      END DO
990
      
991
      Model % Rownonzeros = 0
992
      DO i=1,LocalNodes
993
        j = List(i) % Degree
994
        Model % RowNonzeros(InvPerm(i)) = j
995
      END DO      
996
      DEALLOCATE( InvPerm ) 
997
    END IF
998

999
  CONTAINS
1000

1001
    ! Takes elemental indexes (for nodes, edges, faces etc.) and using the initial permutation
1002
    ! pack the nonzeros and then sort them and finally add them to the list matrix structure.
1003
    ! Adding the whole row is computationally advantageous to adding just one entry at a time.
1004
    !-----------------------------------------------------------------------------------------
1005
    SUBROUTINE PackSortAdd(n,Ind,Perm)
1006
      INTEGER :: n
1007
      INTEGER :: Ind(:),Perm(:)
1008
      INTEGER :: i,j,m
1009
      
1010
      m = 0
1011
      DO i=1,n
1012
        j = Perm(Ind(i))
1013
        IF(j==0) CYCLE
1014
        m = m+1
1015
        Ind(m) = j
1016
      END DO
1017

1018
      CALL Sort(m,Ind)
1019

1020
      DO i=1,m
1021
        CALL List_AddMatrixIndexes(List,Ind(i),m,Ind)
1022
      END DO
1023

1024
    END SUBROUTINE PackSortAdd
1025
        
1026

1027

1028
#if 0    
1029
    SUBROUTINE DgRadiationIndexes(Element,n,ElemInds)
1030

1031
      TYPE(Element_t), POINTER :: Element
1032
      INTEGER :: n
1033
      INTEGER :: ElemInds(:)
1034

1035
      TYPE(Element_t), POINTER :: Left, Right, Parent
1036
      INTEGER :: i,i1,n1,mat_id
1037
      LOGICAL :: Found
1038
           
1039
      Left => Element % BoundaryInfo % Left
1040
      Right => Element % BoundaryInfo % Right
1041
      
1042
      IF(ASSOCIATED(Left) .AND. ASSOCIATED(Right)) THEN
1043
        DO i=1,2
1044
          IF(i==1) THEN
1045
            Parent => Left
1046
          ELSE
1047
            Parent => Right
1048
          END IF
1049
          
1050
          mat_id = ListGetInteger( CurrentModel % Bodies(Parent % BodyId) % Values, &
1051
              'Material', Found, minv=1,maxv=CurrentModel % NumberOfMaterials )
1052
          IF( ListCheckPresent( CurrentModel % Materials(mat_id) % Values,'Emissivity') ) EXIT
1053

1054
          IF( i==2) THEN
1055
            CALL Fatal(Caller,'DG boundary parents should have emissivity!')
1056
          END IF          
1057
        END DO                
1058
      ELSE IF(ASSOCIATED(Left)) THEN
1059
        Parent => Left
1060
      ELSE IF(ASSOCIATED(Right)) THEN
1061
        Parent => Right
1062
      ELSE
1063
        CALL Fatal(Caller,'DG boundary should have parents!')
1064
      END IF
1065

1066
      n1 = Parent % TYPE % NumberOfNodes
1067
      DO i=1,n
1068
        DO i1 = 1, n1
1069
          IF( Parent % NodeIndexes( i1 ) == Element % NodeIndexes(i) ) THEN
1070
            ElemInds(i) = Parent % DGIndexes( i1 )
1071
            EXIT
1072
          END IF
1073
        END DO
1074
      END DO
1075
      
1076
    END SUBROUTINE DgRadiationIndexes
1077
#endif
1078
    
1079
!------------------------------------------------------------------------------
1080
  END SUBROUTINE MakeListMatrix
1081
!------------------------------------------------------------------------------
1082

1083

1084
  ! Pick the correct indexes for radition when using discontinuous Galerkin.
1085
  ! In internal BCs we expect to find 'emissivity' given on either side.
1086
  !--------------------------------------------------------------------------
1087
  SUBROUTINE DgRadiationIndexes(Element,n,ElemInds,DiffuseGray)
1088

1089
    TYPE(Element_t), POINTER :: Element
1090
    INTEGER :: n
1091
    INTEGER :: ElemInds(:)
1092
    LOGICAL :: DiffuseGray
1093

1094
    TYPE(Element_t), POINTER :: Left, Right, Parent
1095
    INTEGER :: i,i1,n1,mat_id
1096
    LOGICAL :: Found
1097

1098
    Left => Element % BoundaryInfo % Left
1099
    Right => Element % BoundaryInfo % Right
1100

1101
    IF(ASSOCIATED(Left) .AND. ASSOCIATED(Right)) THEN
1102

1103
      IF( DiffuseGray ) THEN
1104
        DO i=1,2
1105
          IF(i==1) THEN
1106
            Parent => Left
1107
          ELSE
1108
            Parent => Right
1109
          END IF
1110

1111
          mat_id = ListGetInteger( CurrentModel % Bodies(Parent % BodyId) % Values, &
1112
              'Material', Found, minv=1,maxv=CurrentModel % NumberOfMaterials )
1113
          IF(.NOT. Found ) THEN
1114
            CALL Fatal('DGRadiationIndexes','Body '//I2S(Parent % BodyId)//' has no Material associated!')
1115
          END IF
1116
          IF( ListCheckPresent( CurrentModel % Materials(mat_id) % Values,'Emissivity') ) EXIT
1117

1118
          IF( i==2) THEN
1119
            CALL Fatal('DGRadiationIndexes','DG boundary parents should have emissivity!')
1120
          END IF
1121
        END DO
1122
      ELSE
1123
        IF( Left % BodyId > Right % BodyId ) THEN
1124
          Parent => Left
1125
        ELSE
1126
          Parent => Right
1127
        END IF
1128
      END IF
1129
    ELSE IF(ASSOCIATED(Left)) THEN
1130
      Parent => Left
1131
    ELSE IF(ASSOCIATED(Right)) THEN
1132
      Parent => Right
1133
    ELSE
1134
      CALL Fatal('DGRadiationIndexes','DG boundary should have parents!')
1135
    END IF
1136

1137
    n1 = Parent % TYPE % NumberOfNodes
1138
    DO i=1,n
1139
      DO i1 = 1, n1
1140
        IF( Parent % NodeIndexes( i1 ) == Element % NodeIndexes(i) ) THEN
1141
          ElemInds(i) = Parent % DGIndexes( i1 )
1142
          EXIT
1143
        END IF
1144
      END DO
1145
    END DO
1146

1147
  END SUBROUTINE DgRadiationIndexes
1148

1149
  
1150
!------------------------------------------------------------------------------
1151
!> Create a list matrix array given the mesh, the active domains and the elementtype 
1152
!> related to the solver. The list matrix is flexible since it can account 
1153
!> for any entries. Also constraints and periodic BCs may give rise to entries
1154
!> in the list matrix topology. Multithreaded version using ListMatrixArray for
1155
!> matrix storage.
1156
!------------------------------------------------------------------------------
1157
  SUBROUTINE MakeListMatrixArray( Model,Solver,Mesh,List,Reorder, &
1158
        LocalNodes,Equation, DGSolver, GlobalBubbles, &
1159
        NodalDofsOnly, ProjectorDofs, CalcNonZeros )
1160
!------------------------------------------------------------------------------
1161
    TYPE(Model_t)  :: Model
1162
    TYPE(Solver_t) :: Solver
1163
    TYPE(Mesh_t)   :: Mesh
1164
    TYPE(ListMatrixArray_t) :: List
1165
    INTEGER, OPTIONAL :: Reorder(:)
1166
    INTEGER :: LocalNodes
1167
    CHARACTER(LEN=*), OPTIONAL :: Equation
1168
    LOGICAL, OPTIONAL :: DGSolver
1169
    LOGICAL, OPTIONAL :: GlobalBubbles
1170
    LOGICAL, OPTIONAL :: NodalDofsOnly
1171
    LOGICAL, OPTIONAL :: ProjectorDofs
1172
    LOGICAL, OPTIONAL :: CalcNonZeros
1173
!------------------------------------------------------------------------------
1174
    INTEGER :: t,i,j,k,l,m,k1,k2,n,p,q,e1,e2,f1,f2, nReord, EDOFs, FDOFs, BDOFs, This, istat, nthr
1175
    INTEGER :: NDOFs, DOFsPerNode
1176
    LOGICAL :: Flag, FoundDG, GB, Found, Radiation, DoProjectors, DoNonZeros
1177
    INTEGER, ALLOCATABLE :: InvPerm(:)    
1178
    TYPE(Matrix_t), POINTER :: Projector
1179
    INTEGER :: IndexSize, NumberOfFactors
1180
    INTEGER, ALLOCATABLE :: Indexes(:), IndexReord(:), IPerm(:)
1181
    TYPE(ListMatrixEntry_t), POINTER :: CList, Lptr
1182
    TYPE(Matrix_t),POINTER :: PMatrix
1183
    TYPE(Element_t), POINTER :: Element,Elm, Edge1, Edge2, Face1, Face2, ElementsList(:)
1184
    TYPE(Graph_t), POINTER :: CurrentColourList
1185
    INTEGER :: CurrentColour, BoundaryColour, CurrentColourStart, &
1186
          CurrentColourEnd, NumberOfMeshColours
1187
    LOGICAL :: NeedLocking, GotIt
1188
    CHARACTER(:), ALLOCATABLE :: RadiationFlag
1189
    CHARACTER(*), PARAMETER :: Caller = 'MakeListMatrixArray'
1190
!------------------------------------------------------------------------------
1191

1192
    CALL Info(Caller,'Creating list matrix',Level=14)
1193

1194
    GB = .FALSE.
1195
    IF ( PRESENT(GlobalBubbles) ) GB = GlobalBubbles
1196

1197
    ! No colouring equals a single colour
1198
    NumberOfMeshColours = 1
1199
    IF (ASSOCIATED(Solver % ColourIndexList)) THEN
1200
       NumberOfMeshColours = Solver % ColourIndexList % n
1201
       ! If boundary mesh exists, it has been coloured as well
1202
       IF (ASSOCIATED(Solver % BoundaryColourIndexList)) &
1203
            NumberOfMeshColours = NumberOfMeshColours + Solver % BoundaryColourIndexList % n
1204
    END IF
1205

1206
    nthr=1
1207
    !$ nthr = omp_get_max_threads()
1208
    NeedLocking = (NumberOfMeshColours == 1) .AND. (nthr > 1)
1209
    CALL ListMatrixArray_Allocate(List, LocalNodes, Atomic=NeedLocking)
1210
    
1211
    NDOFs = Mesh % MaxNDOFs
1212
    BDOFs = Mesh % MaxBDOFs
1213
    EDOFs = Mesh % MaxEdgeDOFs
1214
    FDOFs = Mesh % MaxFaceDOFs
1215
    IF( PRESENT( NodalDofsOnly ) ) THEN
1216
      IF( NodalDofsOnly ) THEN
1217
        EDOFS = 0
1218
        FDOFS = 0
1219
      END IF
1220
    END IF
1221

1222
    IF( PRESENT( ProjectorDofs ) ) THEN
1223
      DoProjectors = ProjectorDofs 
1224
    ELSE
1225
      DoProjectors = .TRUE.
1226
    END IF
1227

1228
    IF( EDOFS > 0 .AND. .NOT. ASSOCIATED(Mesh % Edges) ) THEN
1229
      CALL Warn(Caller,'Edge dofs requested but not edges exist in mesh!')
1230
      EDOFS = 0
1231
    END IF
1232

1233
    IF( FDOFS > 0 .AND. .NOT. ASSOCIATED(Mesh % Faces) ) THEN
1234
      CALL Warn(Caller,'Face dofs requested but not faces exist in mesh!')
1235
      FDOFS = 0
1236
    END IF
1237

1238
    ! Create the permutation for the Discontinuous Galerkin solver 
1239
    !-------------------------------------------------------------------
1240
    FoundDG = .FALSE.
1241
    IF ( DGSolver ) THEN
1242
       IndexSize = 128
1243
       ALLOCATE( Indexes(IndexSize), STAT=istat )
1244
       IF( istat /= 0 ) CALL Fatal(Caller,'Allocation error for Indexes')
1245

1246
       ! TODO: Add multithreading
1247
       DO t=1,Mesh % NumberOfEdges
1248
         n = 0
1249
         Elm => Mesh % Edges(t) % BoundaryInfo % Left
1250
         IF ( ASSOCIATED( Elm ) ) THEN
1251
             IF ( CheckElementEquation(Model,Elm,Equation) ) THEN
1252
                FoundDG = FoundDG .OR. Elm % DGDOFs > 0
1253
                DO j=1,Elm % DGDOFs
1254
                   n = n + 1
1255
                   Indexes(n) = Elm % DGIndexes(j)
1256
                END DO
1257
             END IF
1258
         END IF
1259

1260
         Elm => Mesh % Edges(t) % BoundaryInfo %  Right
1261
         IF ( ASSOCIATED( Elm ) ) THEN
1262
             IF ( CheckElementEquation(Model,Elm,Equation) ) THEN
1263
                FoundDG = FoundDG .OR. Elm % DGDOFs > 0
1264
                DO j=1,Elm % DGDOFs
1265
                   n = n + 1
1266
                   Indexes(n) = Elm % DGIndexes(j)
1267
                END DO
1268
             END IF
1269
         END IF
1270

1271
         DO i=1,n
1272
            k1 = Reorder(Indexes(i))
1273
            IF ( k1 <= 0 ) CYCLE
1274
            DO j=1,n
1275
              k2 = Reorder(Indexes(j))
1276
              IF ( k2 <= 0 ) CYCLE
1277
              
1278
              CALL ListMatrixArray_AddEntry(List, k1, k2)
1279
            END DO
1280
         END DO
1281
      END DO
1282
      DO t=1,Mesh % NumberOfFaces
1283
         n = 0
1284
         Elm => Mesh % Faces(t) % BoundaryInfo % Left
1285
         IF ( ASSOCIATED( Elm ) ) THEN
1286
             IF ( CheckElementEquation(Model,Elm,Equation) ) THEN
1287
                FoundDG = FoundDG .OR. Elm % DGDOFs > 0
1288
                DO j=1,Elm % DGDOFs
1289
                   n = n + 1
1290
                   Indexes(n) = Elm % DGIndexes(j)
1291
                END DO
1292
             END IF
1293
         END IF
1294

1295
         Elm => Mesh % Faces(t) % BoundaryInfo %  Right
1296
         IF ( ASSOCIATED( Elm ) ) THEN
1297
             IF ( CheckElementEquation(Model,Elm,Equation) ) THEN
1298
                FoundDG = FoundDG .OR. Elm % DGDOFs > 0
1299
                DO j=1,Elm % DGDOFs
1300
                   n = n + 1
1301
                   Indexes(n) = Elm % DGIndexes(j)
1302
                END DO
1303
             END IF
1304
         END IF
1305

1306
         DO i=1,n
1307
            k1 = Reorder(Indexes(i))
1308
            IF ( k1 <= 0 ) CYCLE
1309
            DO j=1,n
1310
              k2 = Reorder(Indexes(j))
1311
              IF ( k2 <= 0 ) CYCLE
1312
              
1313
              CALL ListMatrixArray_AddEntry(List, k1, k2)
1314
            END DO
1315
         END DO
1316
      END DO
1317
      DEALLOCATE(Indexes)
1318
    END IF
1319

1320
    ! If this is not a GD solver then create permutation considering 
1321
    ! nodal, edge, face and bubble dofs. 
1322
    !-------------------------------------------------------------------
1323
    IF (.NOT. FoundDG) THEN
1324
      
1325
      !$OMP PARALLEL &
1326
      !$OMP SHARED(LocalNodes, List, Equation, NDOFs, EDOFS, FDOFS, GB, &
1327
      !$OMP        Reorder, Model, Mesh, NumberOfMeshColours, CurrentColourStart, NeedLocking, &
1328
      !$OMP        CurrentColourEnd, CurrentColourList, ElementsList, Solver, BoundaryColour) &
1329
      !$OMP PRIVATE(Element, Indexes, istat, IndexSize, IndexReord, &
1330
      !$OMP         IPerm, n, i, j, nReord, k1, k2, Lptr, DOFsPerNode, &
1331
      !$OMP         CurrentColour, Edge1, Face1) &
1332
      !$OMP DEFAULT(NONE)
1333

1334
      IndexSize = 0
1335

1336
      ! Loop over mesh colours
1337
      DO CurrentColour=1,NumberOfMeshColours
1338

1339
         !$OMP SINGLE
1340
         ! Test if matrix is actually coloured or not
1341
         IF (NumberOfMeshColours == 1) THEN
1342
           CurrentColourList => NULL()
1343
           ElementsList => Mesh % Elements(1:Mesh % NumberOfBulkElements+Mesh % NumberOFBoundaryElements)
1344
           CurrentColourStart = 1
1345
           CurrentColourEnd = Mesh % NumberOfBulkElements+Mesh % NumberOFBoundaryElements
1346
         ELSE IF (CurrentColour <= Solver % ColourIndexList % n) THEN
1347
           CALL Info(Caller,'ListMatrix add colour: '//I2S(CurrentColour),Level=10)
1348
           CurrentColourList => Solver % ColourIndexList
1349
           ElementsList => Mesh % Elements(1:Mesh % NumberOfBulkElements)
1350
           CurrentColourStart = CurrentColourList % Ptr(CurrentColour)
1351
           CurrentColourEnd = CurrentColourList % Ptr(CurrentColour+1)-1
1352
         ELSE
1353
           BoundaryColour = CurrentColour-Solver % ColourIndexList % n
1354
           CALL Info(Caller,'ListMatrix add boundary colour: '//I2S(BoundaryColour),Level=10)
1355

1356
           CurrentColourList => Solver % BoundaryColourIndexList
1357
           ! Boundary elements are stored after bulk elements in Mesh
1358
           ElementsList => Mesh % Elements(Mesh % NumberOfBulkElements+1:&
1359
                 Mesh % NumberOfBulkElements + Mesh % NumberOfBoundaryElements)
1360
           CurrentColourStart = CurrentColourList % Ptr(BoundaryColour)
1361
           CurrentColourEnd = CurrentColourList % Ptr(BoundaryColour+1)-1
1362
         END IF
1363
         !$OMP END SINGLE
1364
         
1365
         !$OMP DO
1366
         DO t=CurrentColourStart,CurrentColourEnd
1367
            IF (ASSOCIATED(CurrentColourList)) THEN
1368
               Element => ElementsList(CurrentColourList % ind(t))
1369
            ELSE
1370
               Element => ElementsList(t)
1371
            END IF
1372

1373

1374
            IF ( PRESENT(Equation) ) THEN
1375
               IF ( .NOT. CheckElementEquation(Model,Element,Equation) ) CYCLE
1376
            END IF
1377

1378
            n = Element % NDOFs 
1379
            IF( EDOFS > 0 ) n = n + Element % TYPE % NumberOfEdges * EDOFs 
1380
            IF( FDOFS > 0 ) n = n + Element % TYPE % NumberOfFaces * FDOFs
1381
            IF ( GB ) n = n + Element % BDOFs
1382
            
1383
            ! Reallocate index permutation if needed
1384
            IF ( n > IndexSize ) THEN
1385
               IF (ALLOCATED(Indexes)) DEALLOCATE(Indexes, IndexReord, IPerm)
1386
               
1387
               IndexSize = MAX(MAX(128, IndexSize*2), n)
1388
               ALLOCATE(Indexes(IndexSize), &
1389
                   IndexReord(IndexSize), &
1390
                   IPerm(IndexSize), STAT=istat )
1391
               IF( istat /= 0 ) CALL Fatal(Caller,'Allocation error for Indexes of size: '//I2S(n))
1392
            END IF
1393
            
1394
            n = 0
1395
            DOFsPerNode = Element % NDOFs / Element % TYPE % NumberOfNodes
1396
            IF (DOFsPerNode > 0) THEN
1397
              DO i=1,Element % TYPE % NumberOfNodes
1398
                DO j=1,DOFsPerNode
1399
                  n = n + 1
1400
                  Indexes(n) = NDOFs * (Element % NodeIndexes(i)-1) + j
1401
                END DO
1402
              END DO
1403
            END IF
1404
!            DO i=1,Element % NDOFs
1405
!               n = n + 1
1406
!               Indexes(n) = Element % NodeIndexes(i)
1407
!            END DO
1408
            
1409
            IF ( EDOFs > 0 ) THEN
1410
               IF ( ASSOCIATED(Element % EdgeIndexes) ) THEN
1411
                  DO j=1,Element % TYPE % NumberOFEdges
1412
                     DO i=1, Mesh % Edges(Element % EdgeIndexes(j)) % BDOFs
1413
                        n = n + 1
1414
                        Indexes(n) = EDOFs * (Element % EdgeIndexes(j)-1) + i &
1415
                             + NDOFs * Mesh % NumberOfNodes
1416
                     END DO
1417
                  END DO
1418
                  IF ( GB .AND. ASSOCIATED(Element % BoundaryInfo)) THEN
1419
                    Edge1 => Mesh % Edges(Element % EdgeIndexes(1))
1420
                    IF(Element % Type % ElementCode==Edge1 % Type % ElementCode) THEN
1421
                      DO i=1, Element % BDOFs
1422
                        n = n + 1
1423
                        Indexes(n) = EDOFs*(Element % EdgeIndexes(1)-1) + i + &
1424
                            NDOFs * Mesh % NumberOfNodes
1425
                      END DO
1426
                    END IF
1427
                  END IF
1428
               END IF
1429
            END IF
1430
            
1431
            IF ( FDOFS > 0 ) THEN
1432
               IF ( ASSOCIATED(Element % FaceIndexes) ) THEN
1433
                  DO j=1,Element % TYPE % NumberOFFaces
1434
                     DO i=1, Mesh % Faces(Element % FaceIndexes(j)) % BDOFs
1435
                        n = n + 1
1436
                        Indexes(n) = FDOFs*(Element % FaceIndexes(j)-1) + i + &
1437
                             NDOFs * Mesh % NumberOfNodes + EDOFs*Mesh % NumberOfEdges
1438
                     END DO
1439
                  END DO
1440
               END IF
1441

1442
               IF ( GB.AND. ASSOCIATED(Element % BoundaryInfo) ) THEN
1443
                 Face1 => Mesh % Faces(Element % FaceIndexes(1))
1444
                 IF(Element % Type % ElementCode==Face1 % Type % ElementCode) THEN
1445
                   DO i=1, Element % BDOFs
1446
                     n = n + 1
1447
                     Indexes(n) = FDOFs*(Element % FaceIndexes(1)-1) + i + &
1448
                         NDOFs * Mesh % NumberOfNodes + EDOFs*Mesh % NumberOfEdges
1449
                   END DO
1450
                 END IF
1451
               END IF
1452
            END IF
1453
            
1454
            IF ( GB .AND. ASSOCIATED(Element % BubbleIndexes) ) THEN
1455
               DO i=1,Element % BDOFs
1456
                  n = n + 1
1457
                  Indexes(n) = FDOFs*Mesh % NumberOfFaces + &
1458
                       NDOFs * Mesh % NumberOfNodes + EDOFs*Mesh % NumberOfEdges + &
1459
                       Element % BubbleIndexes(i)
1460
               END DO
1461
            END IF
1462

1463
            nReord=0
1464
            DO i=1,n
1465
               IF (Reorder(Indexes(i)) > 0) THEN
1466
                  nReord=nReord+1
1467
                  IndexReord(nReord)=Reorder(Indexes(i))
1468
               END IF
1469
            END DO
1470
            CALL InsertionSort(nReord, IndexReord, IPerm)
1471

1472
            DO i=1,nReord
1473
               k1 = IndexReord(IPerm(i))
1474
             
1475
               ! Bulk add all sorted nonzero indices to a matrix row in one go
1476
               CALL ListMatrixArray_AddEntries(List, k1, nReord, IndexReord, IPerm, NeedLocking)
1477
            END DO
1478
         END DO
1479
         !$OMP END DO
1480
      END DO ! Loop over mesh colours
1481

1482
      IF (ALLOCATED(Indexes)) DEALLOCATE(Indexes, IndexReord, IPerm)
1483
      !$OMP END PARALLEL
1484

1485
!
1486
!     Diffuse gray radiation condition:
1487
!     ---------------------------------
1488
      Radiation = ListGetLogical( Solver % Values, 'Radiation Solver', Found )
1489
      IF( Radiation ) THEN        
1490
        Radiation = .FALSE.
1491
        DO i=1,Model % NumberOfBCs
1492
          RadiationFlag = ListGetString( Model % BCs(i) % Values, 'Radiation', GotIt )
1493
          IF (GotIt) THEN
1494
            IF ( RadiationFlag == 'diffuse gray' ) THEN
1495
              Radiation = .TRUE.
1496
              EXIT
1497
            END IF
1498
          END IF
1499
        END DO
1500
      END IF
1501
      
1502
      IF ( Radiation ) THEN
1503
        CALL Info(Caller,'Adding radiation matrix entries',Level=14)
1504
        DO i = Mesh % NumberOfBulkElements+1, &
1505
          Mesh % NumberOfBulkElements + Mesh % NumberOfBoundaryElements
1506

1507
          Element => Mesh % Elements(i)
1508
          IF ( ASSOCIATED(Element % BoundaryInfo % RadiationFactors) ) THEN
1509
             DO j=1,Element % TYPE % NumberOfNodes
1510
                k1 = Reorder(Element % NodeIndexes(j))
1511

1512
                NumberOfFactors = Element % BoundaryInfo % &
1513
                  RadiationFactors % NumberOfImplicitFactors
1514

1515
                DO n=1,NumberOfFactors
1516

1517
                  Elm => Mesh % Elements( Element % BoundaryInfo % &
1518
                              RadiationFactors % Elements(n) )
1519

1520
                  DO k=1,Elm % TYPE % NumberOfNodes
1521
                     k2 = Reorder( Elm % NodeIndexes(k) )
1522
                     CALL ListMatrixArray_AddEntry(List, k1, k2)
1523
                  END DO
1524
                END DO
1525
             END DO
1526
          END IF
1527
        END DO
1528
        CALL Info(Caller,'Done Adding radiation matrix',Level=14)
1529
      END IF
1530

1531
      ! TODO: Add multithreading
1532
      DO i=Mesh % NumberOfBulkElements+1, Mesh % NumberOfBulkElements+ &
1533
                     Mesh % NumberOfBoundaryElements
1534
        IF ( Mesh % Elements(i) % TYPE % ElementCode <  102 .OR. &
1535
             Mesh % Elements(i) % TYPE % ElementCode >= 200 ) CYCLE
1536

1537
        k1 = Reorder( Mesh % Elements(i) % NodeIndexes(1) )
1538
        IF ( k1 > 0 ) THEN
1539
          DO k=1,Mesh % Elements(i) % TYPE % NumberOFNodes
1540
            k2 = Reorder( Mesh % Elements(i) % NodeIndexes(k) )
1541
            IF ( k2 > 0 ) THEN
1542
              CALL ListMatrixArray_AddEntry(List, k1, k2)
1543
              CALL ListMatrixArray_AddEntry(List, k2, k1)
1544
            END IF
1545
          END DO
1546
        END IF
1547

1548
        ! This is a connection element, make a matrix connection for that 
1549
        IF ( Mesh % Elements(i) % TYPE % ElementCode == 102 ) THEN
1550
          k2 = Reorder( Mesh % Elements(i) % NodeIndexes(2) )
1551
          IF ( k2 > 0 ) CALL ListMatrixArray_AddEntry(List, k2, k2)
1552
        END IF
1553
      END DO
1554

1555

1556
      ! Add connection from projectors. These are only needed if the projector is treated
1557
      ! implicely. For explicit projectors or when using Lagrange coefficients the 
1558
      ! connections are not needed. 
1559
      !----------------------------------------------------------------------------------
1560
      IF( DoProjectors ) THEN
1561
        DO This=1,Model % NumberOfBCs
1562
          Projector => Model % BCs(This) % PMatrix
1563
          IF ( .NOT. ASSOCIATED(Projector) ) CYCLE
1564

1565
          IF( ListGetLogical( Model % BCs(This) % Values,&
1566
              'Periodic BC Explicit',Found)) CYCLE
1567
          IF( ListGetLogical( Model % BCs(This) % Values,&
1568
              'Periodic BC Use Lagrange Coefficient',Found)) CYCLE
1569

1570
          CALL Info(Caller,'Adding matrix topology for BC: '//I2S(This),Level=10)
1571

1572
          ! TODO: Add multithreading
1573
          DO i=1,Projector % NumberOfRows
1574
            k = Reorder( Projector % InvPerm(i) )
1575
            IF ( k > 0 ) THEN
1576
              DO l=Projector % Rows(i),Projector % Rows(i+1)-1
1577
                IF ( Projector % Cols(l) <= 0 ) CYCLE
1578
                m = Reorder( Projector % Cols(l) )
1579
                IF ( m > 0 ) THEN
1580
                  CALL ListMatrixArray_AddEntry(List, k, m)
1581
                  CALL ListMatrixArray_AddEntry(List, m, k)
1582
                  CList => List % Rows(k) % Head
1583
                  DO WHILE( ASSOCIATED( CList ) )
1584
                    CALL ListMatrixArray_AddEntry(List, m, CList % Index)
1585
                    CALL ListMatrixArray_AddEntry(List, CList % Index, m)
1586
                    CList => CList % Next
1587
                  END DO
1588
                END IF
1589
              END DO
1590
            END IF
1591
          END DO
1592
        END DO
1593
      END IF ! DoProjectors
1594

1595
    END IF ! (.NOT. FoundDG)
1596

1597

1598

1599
    Model % TotalMatrixElements = 0
1600
    DO i=1,LocalNodes
1601
      j = List % Rows(i) % Degree
1602
      Model % TotalMatrixElements = Model % TotalMatrixElements + j
1603
    END DO
1604
    
1605
      
1606
    DoNonZeros = .TRUE.
1607
    IF( PRESENT( CalcNonZeros ) ) DoNonZeros = CalcNonZeros
1608
    
1609
    IF( DoNonZeros ) THEN
1610
      ALLOCATE( InvPerm(LocalNodes) )
1611
      InvPerm = 0
1612
      k = 0
1613
      DO i=1,SIZE(Reorder)
1614
        IF (Reorder(i)>0) THEN
1615
          k = k + 1
1616
          InvPerm(Reorder(i)) = k
1617
        END IF
1618
      END DO
1619
      
1620
      Model % Rownonzeros = 0
1621
      DO i=1,LocalNodes
1622
        j = List % Rows(i) % Degree
1623
        Model % RowNonzeros(InvPerm(i)) = j
1624
      END DO      
1625
      DEALLOCATE( InvPerm ) 
1626
    END IF
1627

1628
!------------------------------------------------------------------------------
1629
  END SUBROUTINE MakeListMatrixArray
1630
!------------------------------------------------------------------------------
1631

1632
  
1633

1634
!------------------------------------------------------------------------------
1635
!>    Initialize a CRS format matrix to the effect that it will be ready to
1636
!>    accept values when CRS_GlueLocalMatrix is called (build up the index
1637
!>    tables of a CRS format matrix)....
1638
!------------------------------------------------------------------------------
1639
  SUBROUTINE InitializeMatrix( Matrix, n, List, DOFs, Reorder, InvInitialReorder )
1640
!------------------------------------------------------------------------------
1641
    INTEGER :: DOFs, n
1642
    TYPE(Matrix_t),POINTER :: Matrix
1643
    TYPE(ListMatrix_t) :: List(:)
1644
    INTEGER, OPTIONAL :: Reorder(:), InvInitialReorder(:)
1645
!------------------------------------------------------------------------------
1646
    TYPE(ListMatrixEntry_t), POINTER :: Clist
1647
    INTEGER :: i,j,k,l,m,k1,k2
1648
    INTEGER, POINTER CONTIG :: Rows(:), Cols(:)
1649
    LOGICAL :: DoReorder 
1650

1651
    
1652
    Rows => Matrix % Rows
1653
    Cols => Matrix % Cols
1654
    
1655
    DoReorder = .FALSE.
1656
    IF( PRESENT( Reorder ) ) THEN
1657
      IF( .NOT. PRESENT( InvInitialReorder ) ) THEN
1658
        CALL Fatal('InitializeMatrix','Need both old and new numbering!')
1659
      END IF
1660
      DoReorder = .TRUE.
1661
    END IF
1662
    
1663
    
1664
    IF( DoReorder ) THEN
1665
      ! In case of reordering we need to compute the row offset in advance
1666
      Rows(1) = 1
1667
      DO i=1,n       
1668
        DO l=1,DOFs
1669
          j = Reorder( InvInitialReorder(i) )
1670
          k1 = DOFs * (j-1) + l
1671
          Rows(k1+1) = Dofs * List(i) % Degree
1672
        END DO
1673
      END DO
1674
      DO i=1,Dofs*n
1675
        Rows(i+1) = Rows(i) + Rows(i+1)
1676
      END DO
1677
      
1678
      !$OMP PARALLEL DO SHARED(Rows, Cols, List, n, DOFs, Reorder, InvInitialReorder) &
1679
      !$OMP PRIVATE(CList, l, j, k1, k2, k, m) &
1680
      !$OMP DEFAULT(NONE)
1681
      DO i=1,n
1682
        DO l=1,DOFs
1683
          CList => List(i) % Head
1684
          j = Reorder( InvInitialReorder(i) )
1685
          k1 = DOFs * (j-1) + l
1686
          k2 = Rows(k1)-1
1687
          
1688
          DO WHILE( ASSOCIATED( CList ) )
1689
            k = Reorder( InvInitialReorder(Clist % INDEX))
1690
            k = DOFs*(k-1)
1691
            DO m=k+1,k+DOFs
1692
              k2 = k2+1
1693
              Cols(k2) = m
1694
            END DO
1695
            CList => Clist % Next
1696
          END DO
1697
        END DO
1698
      END DO
1699
      !$OMP END PARALLEL DO      
1700
    ELSE
1701
      Rows(1) = 1
1702
      DO i=1,n       
1703
        DO l=1,DOFs
1704
          k1 = DOFs * (i-1) + l
1705
          Rows(k1+1) = Dofs * List(i) % Degree
1706
        END DO
1707
      END DO
1708
      DO i=1,Dofs*n
1709
        Rows(i+1) = Rows(i) + Rows(i+1)
1710
      END DO
1711

1712
      ! If there is no renumbering then the reordering is one-to-one mapping
1713
      !$OMP PARALLEL DO SHARED(Rows, Cols, List, n, DOFs) &
1714
      !$OMP PRIVATE(CList, l, j, k1, k2, k, m) &
1715
      !$OMP DEFAULT(NONE)
1716

1717
      DO i=1,n
1718
        DO l=1,DOFs
1719
          CList => List(i) % Head
1720
          j = i
1721
          k1 = DOFs * (j-1) + l
1722
          k2 = Rows(k1)-1
1723
          
1724
          DO WHILE( ASSOCIATED( CList ) )
1725
            k = Clist % index 
1726
            k = DOFs*(k-1)
1727
            DO m=k+1,k+DOFs
1728
              k2 = k2+1
1729
              Cols(k2) = m
1730
            END DO
1731
            CList => Clist % Next
1732
          END DO
1733
        END DO
1734
      END DO
1735
      !$OMP END PARALLEL DO      
1736
    END IF
1737
    
1738
    IF ( Matrix % FORMAT == MATRIX_CRS ) CALL CRS_SortMatrix( Matrix )
1739
!------------------------------------------------------------------------------
1740
  END SUBROUTINE InitializeMatrix
1741
!------------------------------------------------------------------------------
1742

1743

1744
!------------------------------------------------------------------------------
1745
   FUNCTION CreateMatrix( Model, Solver, Mesh, Perm, DOFs, MatrixFormat, &
1746
          OptimizeBW, Equation, DGSolver, GlobalBubbles, &
1747
          NodalDofsOnly, ProjectorDofs, ThreadedStartup, &
1748
          UseGivenPerm ) RESULT(Matrix)
1749
!------------------------------------------------------------------------------
1750
     IMPLICIT NONE
1751
     TYPE(Model_t) :: Model
1752
     TYPE(Mesh_t) :: Mesh
1753
     TYPE(Solver_t), TARGET :: Solver
1754
     INTEGER :: DOFs, MatrixFormat
1755
     INTEGER, TARGET :: Perm(:)
1756
     LOGICAL :: OptimizeBW
1757
     LOGICAL, OPTIONAL :: DGSolver, GlobalBubbles
1758
     LOGICAL, OPTIONAL :: NodalDofsOnly, ProjectorDofs
1759
     LOGICAL, OPTIONAL :: ThreadedStartup
1760
     LOGICAL, OPTIONAL :: USeGivenPerm
1761

1762
     CHARACTER(LEN=*), OPTIONAL :: Equation
1763

1764
     TYPE(Matrix_t),POINTER :: Matrix
1765
!------------------------------------------------------------------------------
1766
     TYPE(ListMatrix_t), POINTER :: ListMatrix(:)
1767
     TYPE(ListMatrixArray_t) :: ListMatrixArray
1768
     TYPE(Matrix_t), POINTER :: A
1769
     TYPE(Element_t), POINTER :: Element
1770
     TYPE(ListMatrixEntry_t), POINTER :: CList
1771
     CHARACTER(:), ALLOCATABLE :: Eq,str
1772
     LOGICAL :: GotIt, DG, GB, UseOptimized, Found, UseGiven
1773
     INTEGER i,j,k,l,k1,t,n, p,m, minEdgeDOFs, maxEdgeDOFs, &
1774
           minFaceDOFs, maxFaceDOFs, BDOFs, cols, istat, &
1775
           NDOFs
1776
     INTEGER, POINTER :: Ivals(:)
1777
     INTEGER, ALLOCATABLE, SAVE :: InvInitialReorder(:)
1778
     INTEGER :: nthr
1779
     LOGICAL :: UseThreads
1780
     LOGICAL, ALLOCATABLE :: ConstrainedNode(:)
1781
     CHARACTER(*), PARAMETER :: Caller = 'CreateMatrix'
1782
   
1783
!------------------------------------------------------------------------------
1784

1785
     NULLIFY( Matrix )
1786

1787
     DG = .FALSE.
1788
     IF ( PRESENT(DGSolver) )  DG = DGSolver
1789

1790
     GB = .FALSE.
1791
     IF ( PRESENT(GlobalBubbles) ) GB = GlobalBubbles
1792

1793
     UseGiven = .FALSE.
1794
     IF( PRESENT( UseGivenPerm ) ) THEN
1795
       IF( UseGivenPerm ) THEN
1796
         CALL Info(Caller,'Using given Perm table to create the matrix',Level=6)
1797
         UseGiven = .TRUE.
1798
         OptimizeBW = .FALSE.
1799
       END IF
1800
     END IF
1801
       
1802
     IF( OptimizeBW ) THEN
1803
       IF( ListGetLogical( Solver % Values,'DG Reduced Basis',Found ) ) THEN
1804
         CALL Info(Caller,'Suppressing bandwidth optimization for discontinuous bodies',Level=8)
1805
         OptimizeBW = .FALSE.
1806
       END IF
1807
       IF( ListGetLogical( Solver % Values,'Apply Conforming BCs',Found ) ) THEN
1808
         CALL Info(Caller,'Suppressing bandwidth optimization for conforming bcs',Level=8)
1809
         OptimizeBW = .FALSE.
1810
       END IF
1811
     END IF
1812

1813
     
1814
     UseThreads = .FALSE.
1815
     nthr = 1
1816
     IF ( PRESENT(ThreadedStartup) ) THEN
1817
#ifdef _OPENMP
1818
       IF (ThreadedStartup) THEN
1819
          UseThreads = .TRUE.
1820
          nthr = omp_get_max_threads()
1821
       END IF
1822
#endif
1823
     END IF
1824

1825
     minEdgeDOFs = HUGE(minEdgeDOFs)
1826
     maxEdgeDOFs = 0
1827
     minFaceDOFs = HUGE(minFaceDOFs)
1828
     maxFaceDOFs = 0
1829
     BDOFs = 0
1830
     NDOFs = 0
1831

1832
     !$OMP PARALLEL SHARED(Mesh, GB) PRIVATE(j, Element) &
1833
     !$OMP          REDUCTION(min:minEdgeDOFs) REDUCTION(max:maxEdgeDOFs) &
1834
     !$OMP          REDUCTION(min:minFaceDOFs) REDUCTION(max:maxFaceDOFs) &
1835
     !$OMP          REDUCTION(max:BDOFs) &
1836
     !$OMP          DEFAULT(NONE) NUM_THREADS(nthr)
1837
     
1838
     !$OMP DO 
1839
     DO i=1,Mesh % NumberOfEdges
1840
        minEdgeDOFs = MIN( minEdgeDOFs, Mesh % Edges(i) % BDOFs )
1841
        maxEdgeDOFs = MAX( maxEdgeDOFs, Mesh % Edges(i) % BDOFs )
1842
     END DO
1843
     !$OMP END DO NOWAIT
1844
     !$OMP DO
1845
     DO i=1,Mesh % NumberOfFaces
1846
        minFaceDOFs = MIN( minFaceDOFs, Mesh % Faces(i) % BDOFs )
1847
        maxFaceDOFs = MAX( maxFaceDOFs, Mesh % Faces(i) % BDOFs )
1848
     END DO
1849
     !$OMP END DO NOWAIT
1850

1851
     IF ( GB ) THEN
1852
       DO i=1,Mesh % NumberOfBoundaryElements
1853
          j = i + Mesh % NumberOfBulkElements
1854
          Element => Mesh % Elements(j)
1855
          IF(Element % Type % ElementCode >= 300) THEN
1856
            minFaceDOFs = MIN( minFaceDOFs, Element % BDOFs )
1857
            maxFaceDOFs = MAX( maxFaceDOFs, Element % BDOFs )
1858
          ELSE
1859
            minEdgeDOFs = MIN( minEdgeDOFs, Element % BDOFs )
1860
            maxEdgeDOFs = MAX( maxEdgeDOFs, Element % BDOFs )
1861
          END IF
1862
       END DO
1863
     END IF
1864

1865
     !$OMP DO
1866
     DO i=1,Mesh % NumberOfBulkElements
1867
        BDOFs = MAX( BDOFs, Mesh % Elements(i) % BDOFs )
1868
     END DO
1869
     !$OMP END DO NOWAIT
1870
     !$OMP END PARALLEL
1871
     DO i=1,Mesh % NumberOfBulkElements
1872
        NDOFs = MAX( NDOFs, Mesh % Elements(i) % NDOFs )
1873
     END DO
1874
     
1875
     Mesh % MaxEdgeDOFs = maxEdgeDOFs
1876
     IF(minEdgeDOFs <= maxEdgeDOFs ) THEN
1877
       Mesh % MinEdgeDOFs = minEdgeDOFs
1878
     ELSE
1879
       Mesh % MinEdgeDOFs = maxEdgeDOFs
1880
     END IF
1881

1882
     Mesh % MaxFaceDOFs = maxFaceDOFs
1883
     IF(minFaceDOFs <= maxFaceDOFs ) THEN
1884
       Mesh % MinFaceDOFs = minFaceDOFs
1885
     ELSE
1886
       Mesh % MinFaceDOFs = maxFaceDOFs
1887
     END IF
1888
     
1889
     Mesh % MaxBDOFs = BDOFs
1890
     
1891
     IF (PRESENT( Equation)) THEN
1892
       n = LEN(Equation)
1893
       ALLOCATE(CHARACTER(n)::Eq)
1894
       n=StringToLowerCase(Eq,Equation)
1895
     ELSE
1896
       Eq = ' '
1897
     END IF
1898

1899
#if 0
1900
     ! This is tentative code for p edge multigrid 
1901
     IF( ListGetLogical( Solver % Values,'quadratic edge ordering',Found ) ) THEN
1902
       n = Mesh % NumberOfNodes
1903
       m = Mesh % NumberOfEdges
1904
       p = Mesh % NumberOfFaces
1905

1906
       k = 0
1907
       ! Order the linear edges
1908
       DO i=1,m
1909
         k = k+1
1910
         Perm(n+2*i-1) = k
1911
       END DO
1912
       ! quadratic edges
1913
       DO i=1,m
1914
         k = k+1
1915
         Perm(n+2*i) = k
1916
       END DO
1917
       ! and finally two dofs for each face
1918
       DO i=1,p
1919
         k = k+1
1920
         Perm(n+2*m+2*i-1) = k
1921
         k = k+1
1922
         Perm(n+2*m+2*i) = k
1923
       END DO
1924

1925
       UseGiven = .TRUE.
1926
       OptimizeBW = .FALSE.
1927
     END IF     
1928
#endif
1929
     
1930
     IF( UseGiven ) THEN
1931
       k = MAXVAL( Perm ) 
1932
       ALLOCATE( InvInitialReorder(k), STAT=istat )
1933
       IF( istat /= 0 ) THEN
1934
         CALL Fatal(Caller,'Allocation error for InvInitialReorder of size: '//I2S(k))
1935
       END IF
1936
       InvInitialReorder = 0
1937
       DO i=1,SIZE(Perm)
1938
         IF (Perm(i)>0) InvInitialReorder(Perm(i)) = i
1939
       END DO
1940
       GOTO 10
1941
     END IF
1942

1943
       
1944
     Perm = 0
1945
     IF ( PRESENT(Equation) ) THEN
1946
       CALL Info(Caller,'Creating initial permutation',Level=14)
1947
       k = InitialPermutation( Perm,Model,Solver,Mesh,Eq,DG,GB )
1948
       IF ( k <= 0 ) THEN
1949
         IF(ALLOCATED(InvInitialReorder)) DEALLOCATE(InvInitialReorder)
1950
         RETURN
1951
       END IF
1952
     ELSE
1953
       k = SIZE( Perm )
1954
     END IF
1955
     
1956
     IF ( k == SIZE(Perm) ) THEN
1957
       IF(PRESENT(NodalDofsOnly)) THEN
1958
         IF(NodalDofsOnly) k=Mesh % NumberOfNodes
1959
       END IF       
1960
       DO i=1,k 
1961
         Perm(i) = i
1962
       END DO
1963
     END IF
1964

1965
     
1966
     IF( ParEnv % PEs > 1 .AND. &
1967
         ListGetLogical( Solver % Values,'Skip Pure Halo Nodes',Found ) ) THEN
1968
       CALL Info(Caller,'Skipping pure halo nodes',Level=14)
1969
       j = 0
1970
       DO i=1,Mesh % NumberOfNodes 
1971
         ! These are pure halo nodes that need not be communicated. They are created only 
1972
         ! for sufficient geometric information on the boundaries.
1973
         IF( .NOT. ANY( ParEnv % Mype == Mesh % ParallelInfo % NeighbourList(i) % Neighbours ) ) THEN
1974
           Perm(i) = 0
1975
         ELSE IF( Perm(i) > 0 ) THEN
1976
           j = j + 1
1977
           Perm(i) = j
1978
         END IF
1979
       END DO
1980
       PRINT *,'Eliminating '//I2S(k-j)//' halo nodes out of '&
1981
           //I2S(k)//' in partition '//I2S(ParEnv % MyPe)
1982
       k = j
1983
     END IF
1984

1985
     
1986
     IF( OptimizeBW ) THEN
1987
       CALL Info(Caller,'Creating inverse of initial order of size: '//I2S(k),Level=14)
1988
       ALLOCATE( InvInitialReorder(k), STAT=istat )
1989
       IF( istat /= 0 ) THEN
1990
         CALL Fatal(Caller,'Allocation error for InvInitialReorder of size: '//I2S(k))
1991
       END IF
1992

1993
       ! We need to keep the initial numbering only in case we optimize the bandwidth!
1994
       InvInitialReorder = 0
1995
       DO i=1,SIZE(Perm)
1996
         IF (Perm(i)>0) InvInitialReorder(Perm(i)) = i
1997
       END DO
1998
     END IF
1999
     
2000
     UseOptimized = ListGetLogical( Solver % Values, &
2001
         'Optimize Bandwidth Use Always', GotIt )
2002
          
2003
10   Matrix => NULL()
2004

2005
     ! check if matrix structures really need to be created:
2006
     ! -----------------------------------------------------
2007
     IF ( ListGetLogical( Solver % Values, 'No matrix',GotIt)) THEN
2008
       IF(ALLOCATED(InvInitialReorder)) DEALLOCATE(InvInitialReorder)
2009
       RETURN
2010
     END IF
2011

2012
     !------------------------------------------------------------------------------
2013
     IF (UseThreads) THEN
2014
       CALL Info(Caller,'Creating threaded list matrix array for equation',Level=14)
2015
       IF ( PRESENT(Equation) ) THEN
2016
         CALL MakeListMatrixArray( Model, Solver, Mesh, ListMatrixArray, Perm, k, Eq, DG, GB,&
2017
               NodalDofsOnly, ProjectorDofs, CalcNonZeros = .FALSE. )
2018
         n = OptimizeBandwidth( ListMatrixArray % Rows, Perm, InvInitialReorder, &
2019
               k, OptimizeBW, UseOptimized, Eq )
2020
       ELSE
2021
         CALL MakeListMatrixArray( Model, Solver, Mesh, ListMatrixArray, Perm, k, &
2022
               DGSolver=DG, GlobalBubbles=GB, NodalDofsOnly=NodalDofsOnly, &
2023
               ProjectorDofs=ProjectorDofs, CalcNonZeros = .FALSE. )
2024
         n = OptimizeBandwidth( ListMatrixArray % Rows, Perm, InvInitialReorder, &
2025
               k, OptimizeBW,UseOptimized, ' ' )
2026
       END IF
2027
       
2028
       !------------------------------------------------------------------------------
2029
       ! Initialize the matrix. Multithreading only supports CRS.
2030
       !------------------------------------------------------------------------------
2031
       CALL Info(Caller,'Initializing list matrix array for equation',Level=14)
2032
       IF ( MatrixFormat == MATRIX_CRS) THEN
2033
         Matrix => CRS_CreateMatrix( DOFs*k, Model % TotalMatrixElements, Ndeg=DOFs, &
2034
             Reorder=Perm, AllocValues=.TRUE., SetRows = .FALSE.)
2035
         Matrix % FORMAT = MatrixFormat
2036
         IF( OptimizeBW ) THEN
2037
           CALL InitializeMatrix( Matrix, k, ListMatrixArray % Rows, &
2038
               DOFs, Perm, InvInitialReorder )
2039
         ELSE
2040
           CALL InitializeMatrix( Matrix, k, ListMatrixArray % Rows, DOFs )          
2041
         END IF
2042
       ELSE
2043
         CALL Fatal(Caller,'Multithreaded startup only supports CRS matrix format')
2044
       END IF
2045
       
2046
       CALL Info(Caller,'Sparse atrix created',Level=14)
2047

2048
       CALL ListMatrixArray_Free( ListMatrixArray )       
2049
     ELSE
2050
       NULLIFY( ListMatrix )
2051
       CALL Info(Caller,'Creating list matrix for equation: '//TRIM(Eq),Level=14)
2052

2053
       IF ( PRESENT(Equation) ) THEN
2054
         CALL MakeListMatrix( Model, Solver, Mesh, ListMatrix, Perm, k, Eq, DG, GB,&
2055
               NodalDofsOnly, ProjectorDofs, CalcNonZeros = .FALSE.)
2056
         n = OptimizeBandwidth( ListMatrix, Perm, InvInitialReorder, &
2057
               k, OptimizeBW, UseOptimized, Eq )
2058
       ELSE
2059
         CALL MakeListMatrix( Model, Solver, Mesh, ListMatrix, Perm, k, &
2060
               DGSolver=DG, GlobalBubbles=GB, NodalDofsOnly=NodalDofsOnly, &
2061
               ProjectorDofs=ProjectorDofs, CalcNonZeros = .FALSE.)
2062
         n = OptimizeBandwidth( ListMatrix, Perm, InvInitialReorder, &
2063
               k, OptimizeBW,UseOptimized, ' ' )
2064
       END IF
2065
       
2066
       !------------------------------------------------------------------------------
2067
       ! Initialize the matrix. 
2068
       !------------------------------------------------------------------------------
2069
       CALL Info(Caller,'Initializing list matrix for equation',Level=14)
2070
       SELECT CASE( MatrixFormat )
2071
       CASE( MATRIX_CRS )
2072
         
2073
         Matrix => CRS_CreateMatrix( DOFs*k, Model % TotalMatrixElements, Ndeg=DOFs, &
2074
             Reorder=Perm, AllocValues=.TRUE., SetRows = .FALSE.)
2075
         Matrix % FORMAT = MatrixFormat
2076
         IF( OptimizeBW ) THEN
2077
           CALL InitializeMatrix( Matrix, k, ListMatrix, &
2078
               DOFs, Perm, InvInitialReorder )
2079
         ELSE
2080
           CALL InitializeMatrix( Matrix, k, ListMatrix, DOFs )
2081
         END IF
2082
      CASE( MATRIX_BAND )
2083
        Matrix => Band_CreateMatrix( DOFs*k, DOFs*n,.FALSE.,.TRUE. )
2084
         
2085
       CASE( MATRIX_SBAND )
2086
         Matrix => Band_CreateMatrix( DOFs*k, DOFs*n,.TRUE.,.TRUE. )
2087
       END SELECT
2088
       CALL Info(Caller,'Sparse matrix created',Level=14)
2089

2090
       CALL List_FreeMatrix( k, ListMatrix )
2091
     END IF
2092
     
2093
     NULLIFY( Matrix % MassValues, Matrix % DampValues, Matrix % Force, Matrix % RHS_im )
2094

2095
     IF (ListGetLogical(Solver % Values, 'Allocate Preconditioning Matrix', GotIt)) THEN
2096
       CALL Info(Caller, 'Allocating for separate preconditioning matrix!', Level=20)
2097
       ALLOCATE(Matrix % PrecValues(SIZE(Matrix % Values)) )
2098
       Matrix % PrecValues = 0.0d0
2099
     END IF
2100

2101
     !------------------------------------------------------------------------------
2102
     Matrix % Solver => Solver
2103
     Matrix % DGMatrix = DG
2104
     Matrix % Subband = DOFs * n
2105
     Matrix % COMPLEX = .FALSE.
2106
     Matrix % FORMAT  = MatrixFormat
2107
!------------------------------------------------------------------------------
2108

2109
     n = ListGetInteger( Solver % Values, 'Constraint DOFs', Found )
2110
     IF ( n>0 ) THEN       
2111
       Matrix % ConstraintMatrix => AllocateMatrix()
2112
       A => Matrix % ConstraintMatrix
2113
       A % NumberOfRows = n
2114
       ALLOCATE( A % Rows(n+1), A % Diag(n), A % RHS(n), &
2115
           ConstrainedNode(Mesh % NumberOfNodes), STAT=istat )
2116
       IF( istat /= 0 ) THEN
2117
         CALL Fatal(Caller,'Allocation error for CRS matrix topology: '//I2S(n))
2118
       END IF
2119

2120
       DO i=1,n
2121
         A % RHS(i:i) = ListGetConstReal( Solver % Values,  &
2122
           'Constraint DOF ' // i2s(i) // ' Value' )
2123
       END DO
2124

2125
       Cols = 0
2126
       A % Rows(1) = 1
2127
       DO i=1,n
2128
         str = 'Constraint DOF '//i2s(i)//' Body' 
2129
         ivals => ListGetIntegerArray( Solver % Values, str, Found )
2130
         IF ( ASSOCIATED(ivals) ) THEN
2131
           ConstrainedNode = .FALSE.
2132
           DO k=1,Solver % Mesh % NumberOfBulkElements
2133
             Element => Solver % Mesh % Elements(k)
2134
             IF ( ALL(ivals /= Element % Bodyid) ) CYCLE
2135
             IF ( ALL( Perm(Element % NodeIndexes) > 0 ) ) &
2136
               ConstrainedNode(Element % NodeIndexes) = .TRUE.
2137
           END DO
2138
           Cols = Cols+DOFs*COUNT(ConstrainedNode)
2139
         END IF
2140

2141
         str = 'Constraint DOF ' // i2s(i) // ' BC'
2142
         Ivals => ListGetIntegerArray( Solver % Values, str, Found )
2143
         IF ( ASSOCIATED(Ivals) ) THEN
2144
           ConstrainedNode = .FALSE.
2145
           DO k=Solver % Mesh % NumberOfBulkElements+1, &
2146
                Solver % Mesh % NumberOfBulkElements+Solver % Mesh % NumberOfBoundaryElements
2147
             Element => Solver % Mesh % Elements(k)
2148
             IF ( ALL(Element % Boundaryinfo % Constraint /= ivals) ) CYCLE
2149
             IF ( ALL( Perm(Element % NodeIndexes) > 0 ) ) &
2150
                 ConstrainedNode(Element % NodeIndexes) = .TRUE.
2151
           END DO
2152
           Cols = Cols+DOFs*COUNT(ConstrainedNode)
2153
         END IF
2154
         A % Rows(i+1) = A % Rows(i)+Cols
2155
       END DO
2156

2157
       ALLOCATE( A % Cols(cols), A % Values(cols), STAT=istat )
2158
       IF( istat /= 0 ) THEN
2159
         CALL Fatal(Caller,'Allocation error for CRS cols and values: '//I2S(cols))
2160
       END IF
2161
       A % Cols = 0
2162
       A % Values = 0
2163

2164
       DO i=1,n
2165
         str = 'Constraint DOF ' // i2s(i) // ' Body' 
2166
         ivals => ListGetIntegerArray( Solver % Values, str, Found )
2167
         IF ( ASSOCIATED(ivals) ) THEN
2168
           DO k=1,Solver % Mesh % NumberOfBulkElements
2169
             Element => Solver % Mesh % Elements(k)
2170
             IF ( ALL(ivals /= Element % Bodyid) ) CYCLE
2171

2172
             IF ( ALL( Perm(Element % NodeIndexes) > 0 ) ) THEN
2173
               DO p=1,Element % TYPE % NumberOfNodes
2174
                 l = Perm(Element % NodeIndexes(p))
2175
                 DO m=1,DOFs 
2176
                   k1 = DOFs*(l-1)+m
2177
                   CALL CRS_MakeMatrixIndex( A,i,k1 )
2178
                 END DO
2179
               END DO
2180
             END IF
2181
           END DO
2182
         END IF
2183

2184
         str = 'Constraint DOF ' // i2s(i) // ' BC'
2185
         ivals => ListGetIntegerArray( Solver % Values, str, Found )
2186
         IF ( ASSOCIATED(ivals) ) THEN
2187
           DO k=Solver % Mesh % NumberOfBulkElements+1, &
2188
                Solver % Mesh % NumberOfBulkElements+Solver % Mesh % NumberOfBoundaryElements
2189
             Element => Solver % Mesh % Elements(k)
2190
             IF ( ALL(ivals /= Element % Boundaryinfo % Constraint) ) CYCLE
2191
             IF ( ALL(Perm(Element % NodeIndexes) > 0) ) THEN
2192
               DO p=1,Element % TYPE % NumberOfNodes
2193
                 l = Perm(Element % NodeIndexes(p))
2194
                 DO m=1,DOFs 
2195
                   k1 = DOFs*(l-1)+m
2196
                   CALL CRS_MakeMatrixIndex( A,i,k1 )
2197
                 END DO
2198
               END DO
2199
             END IF
2200
           END DO
2201
         END IF
2202
       END DO
2203
       CALL CRS_SortMatrix(A)
2204
     END IF
2205

2206
!     DEALLOCATE( Model % RowNonZeros )
2207
     IF( ALLOCATED(InvInitialReorder) ) DEALLOCATE( InvInitialReorder )
2208
!------------------------------------------------------------------------------
2209
   END FUNCTION CreateMatrix
2210
!------------------------------------------------------------------------------
2211

2212

2213
!------------------------------------------------------------------------------
2214
   FUNCTION CreateOdeMatrix( Model, Solver, Dofs ) RESULT(Matrix)
2215
!------------------------------------------------------------------------------
2216
     TYPE(Model_t) :: Model
2217
     TYPE(Solver_t), TARGET :: Solver
2218
     INTEGER :: DOFs
2219
     TYPE(Matrix_t), POINTER :: Matrix
2220
!------------------------------------------------------------------------------
2221
     LOGICAL :: Found
2222
     INTEGER i,j,k
2223
!------------------------------------------------------------------------------
2224

2225
     Matrix => NULL()
2226

2227
     IF ( ListGetLogical( Solver % Values, 'No matrix',Found)) RETURN
2228
     
2229
     ! Create a list matrix that allows for unspecified entries in the matrix 
2230
     ! structure to be introduced.
2231
     Matrix => AllocateMatrix()
2232
     Matrix % FORMAT = MATRIX_LIST
2233
     
2234
     ! Initialize matrix indices
2235
     DO i = 1, Dofs
2236
       DO j = 1, Dofs
2237
          CALL List_AddMatrixIndex(Matrix % ListMatrix, i, j) 
2238
       END DO
2239
     END DO
2240

2241
     CALL List_ToCRSMatrix(Matrix)
2242
     CALL CRS_SortMatrix(Matrix,.TRUE.)
2243
     
2244
     CALL Info('CreateOdeMatrix','Number of rows in ode matrix: '//&
2245
         I2S(Matrix % NumberOfRows), Level=9)
2246
     CALL Info('CreateOdeMatrix','Number of entries in ode matrix: '//&
2247
         I2S(SIZE(Matrix % Cols) ), Level=9)
2248
     
2249
     Matrix % Solver => Solver
2250
     Matrix % DGMatrix = .FALSE.
2251
     Matrix % Subband = DOFs
2252
     Matrix % COMPLEX = .FALSE.
2253
     ! Matrix % FORMAT  = MatrixFormat
2254

2255
!------------------------------------------------------------------------------
2256
   END FUNCTION CreateOdeMatrix
2257
!------------------------------------------------------------------------------
2258

2259

2260
!------------------------------------------------------------------------------
2261
   FUNCTION CreateDiagMatrix( Model, Solver, Dofs, TimeOrder ) RESULT(Matrix)
2262
!------------------------------------------------------------------------------
2263
     TYPE(Model_t) :: Model
2264
     TYPE(Solver_t), TARGET :: Solver
2265
     INTEGER :: DOFs
2266
     TYPE(Matrix_t), POINTER :: Matrix
2267
     INTEGER, OPTIONAL :: TimeOrder
2268
!------------------------------------------------------------------------------
2269
     LOGICAL :: Found
2270
     INTEGER i,j,k
2271
!------------------------------------------------------------------------------
2272

2273
     Matrix => NULL()
2274

2275
     !IF ( ListGetLogical( Solver % Values, 'No matrix',Found)) RETURN
2276
     
2277
     ! Create a list matrix that allows for unspecified entries in the matrix 
2278
     ! structure to be introduced.
2279
     Matrix => AllocateMatrix()
2280
     Matrix % FORMAT = MATRIX_LIST
2281
     
2282
     ! Initialize matrix indices
2283
     DO i = 1, Dofs
2284
       CALL List_AddMatrixIndex(Matrix % ListMatrix, i, i) 
2285
     END DO
2286

2287
     CALL List_ToCRSMatrix(Matrix)
2288
     CALL CRS_SortMatrix(Matrix,.TRUE.)
2289
     
2290
     CALL Info('CreateOdeMatrix','Number of rows in diag matrix: '//&
2291
         I2S(Matrix % NumberOfRows), Level=9)
2292

2293
     IF( PRESENT( TimeOrder ) ) THEN
2294
       IF( TimeOrder >= 1 ) THEN
2295
         ALLOCATE( Matrix % MassValues( SIZE( Matrix % Values ) ) )
2296
         Matrix % MassValues = 0.0_dp
2297
       END IF
2298
       IF( TimeOrder >= 2 ) THEN
2299
         ALLOCATE( Matrix % DampValues( SIZE( Matrix % Values ) ) )
2300
         Matrix % DampValues = 0.0_dp
2301
       END IF
2302
     END IF
2303
     
2304
     Matrix % Solver => Solver
2305
     Matrix % DGMatrix = .FALSE.
2306
     Matrix % Subband = 1
2307
     Matrix % COMPLEX = .FALSE.
2308
     ! Matrix % FORMAT  = MatrixFormat
2309

2310
!------------------------------------------------------------------------------
2311
   END FUNCTION CreateDiagMatrix
2312
!------------------------------------------------------------------------------
2313

2314

2315
#if 1
2316
!------------------------------------------------------------------------------
2317
   SUBROUTINE RotateMatrix( Matrix,Vector,n,DIM,DOFs,NodeIndexes,  &
2318
       Normals,Tangent1,Tangent2 )
2319
!------------------------------------------------------------------------------
2320

2321
    REAL(KIND=dp) :: Matrix(:,:),Vector(:)
2322
    REAL(KIND=dp), POINTER :: Normals(:,:), Tangent1(:,:),Tangent2(:,:)
2323
    INTEGER :: n,DIM,DOFs,NodeIndexes(:)
2324
!------------------------------------------------------------------------------
2325
    INTEGER :: i,j,k,l,ii
2326
    REAL(KIND=dp) :: s,R(n*DOFs,n*DOFs),Q(n*DOFs,n*DOFs),N1(3),T1(3),T2(3)
2327
    LOGICAL :: Found
2328
!------------------------------------------------------------------------------
2329
    DO i=1,MIN(n,SIZE(NodeIndexes))
2330
      ii = NodeIndexes(i)
2331
      IF ( ii<=0 .OR. ii>SIZE(Normals,1) ) CYCLE
2332

2333
      IF(ASSOCIATED(CurrentModel % Mesh % PeriodicPerm)) THEN
2334
        j = CurrentModel % Mesh % PeriodicPerm(i)
2335
        IF(j>0) THEN
2336
          IF( ListGetLogical( CurrentModel % Solver % Values, &
2337
              'Apply Conforming BCs',Found ) ) ii = NodeIndexes(j)
2338
        END IF
2339
      END IF
2340
      
2341
      R = 0.0d0
2342
      DO j=1,n*DOFs
2343
        R(j,j) = 1.0d0
2344
      END DO
2345

2346
      N1 = Normals( ii,: )
2347

2348
      SELECT CASE(DIM)
2349
      CASE (2)
2350
        R(DOFs*(i-1)+1,DOFs*(i-1)+1) =  N1(1)
2351
        R(DOFs*(i-1)+1,DOFs*(i-1)+2) =  N1(2)
2352

2353
        R(DOFs*(i-1)+2,DOFs*(i-1)+1) = -N1(2)
2354
        R(DOFs*(i-1)+2,DOFs*(i-1)+2) =  N1(1)
2355
      CASE (3)
2356
        T1 = Tangent1( ii,: )
2357
        T2 = Tangent2( ii,: )
2358

2359
        R(DOFs*(i-1)+1,DOFs*(i-1)+1) = N1(1)
2360
        R(DOFs*(i-1)+1,DOFs*(i-1)+2) = N1(2)
2361
        R(DOFs*(i-1)+1,DOFs*(i-1)+3) = N1(3)
2362

2363
        R(DOFs*(i-1)+2,DOFs*(i-1)+1) = T1(1)
2364
        R(DOFs*(i-1)+2,DOFs*(i-1)+2) = T1(2)
2365
        R(DOFs*(i-1)+2,DOFs*(i-1)+3) = T1(3)
2366

2367
        R(DOFs*(i-1)+3,DOFs*(i-1)+1) = T2(1)
2368
        R(DOFs*(i-1)+3,DOFs*(i-1)+2) = T2(2)
2369
        R(DOFs*(i-1)+3,DOFs*(i-1)+3) = T2(3)
2370
      END SELECT
2371

2372
            
2373
      DO j=1,n*DOFs
2374
        DO k=1,n*DOFs
2375
          s = 0.0D0
2376
          DO l=1,n*DOFs
2377
            s = s + R(j,l) * Matrix(l,k)
2378
          END DO
2379
          Q(j,k) = s
2380
        END DO
2381
      END DO
2382

2383
      DO j=1,n*DOFs
2384
        DO k=1,n*DOFs
2385
          s = 0.0D0
2386
          DO l=1,n*DOFs
2387
            s = s + Q(j,l) * R(k,l)
2388
          END DO
2389
          Matrix(j,k) = s
2390
        END DO
2391
      END DO
2392

2393
      DO j=1,n*DOFs
2394
        s = 0.0D0
2395
        DO k=1,n*DOFs
2396
          s = s + R(j,k) * Vector(k)
2397
        END DO
2398
        Q(j,1) = s
2399
      END DO
2400
      Vector(1:n*DOFs) = Q(:,1)
2401
    END DO
2402
!------------------------------------------------------------------------------
2403
  END SUBROUTINE RotateMatrix
2404
!------------------------------------------------------------------------------
2405
#else
2406

2407
! This should be the same as above but more economical but it does not work...
2408
!------------------------------------------------------------------------------
2409
  SUBROUTINE RotateMatrix( Matrix,Vector,n,DIM,DOFs,NodeIndexes,  &
2410
                   Normals,Tangent1,Tangent2 )
2411
!------------------------------------------------------------------------------
2412

2413
    REAL(KIND=dp) :: Matrix(:,:),Vector(:)
2414
    REAL(KIND=dp), POINTER :: Normals(:,:), Tangent1(:,:),Tangent2(:,:)
2415
    INTEGER :: n,DIM,DOFs,NodeIndexes(:)
2416
!------------------------------------------------------------------------------
2417

2418
    INTEGER :: i,ii,j,k,l
2419
    REAL(KIND=dp) :: s,R(DOFs,DOFs),Force0(Dofs),Force(Dofs),SubMat(Dofs,Dofs), &
2420
        SubMat0(Dofs,Dofs),N1(dofs),T1(dofs),T2(dofs)
2421
    INTEGER :: iInds(n),jInds(n)
2422
    LOGICAL :: Found
2423
!------------------------------------------------------------------------------
2424
    DO i=1,MIN(n,SIZE(NodeIndexes))
2425
      ii = NodeIndexes(i)
2426
      IF ( ii <= 0 .OR. ii > SIZE(Normals,1) ) CYCLE
2427

2428
      IF(ASSOCIATED(CurrentModel % Mesh % PeriodicPerm)) THEN
2429
        j = CurrentModel % Mesh % PeriodicPerm(i)
2430
        IF(j>0) THEN
2431
          IF( ListGetLogical( CurrentModel % Solver % Values, &
2432
              'Apply Conforming BCs',Found ) ) ii = NodeIndexes(j)
2433
        END IF
2434
      END IF
2435
      
2436
      SELECT CASE(DIM)
2437
      CASE (2)
2438
        R(1,1:2) = Normals(ii,1:2)
2439
        R(2,1) = -R(1,2)
2440
        R(2,2) = R(1,1)
2441
      CASE (3)
2442
        R(1,1:3) = Normals(ii,:)
2443
        R(2,1:3) = Tangent1(ii,:)
2444
        R(3,1:3) = Tangent2(ii,:)
2445
      END SELECT
2446

2447
      DO k=1,Dofs
2448
        iInds(k) = Dofs*(i-1)+k
2449
      END DO
2450

2451
      DO j=1,n
2452
        DO k=1,Dofs
2453
          jInds(k) = Dofs*(j-1)+k
2454
        END DO
2455

2456
        SubMat0 = Matrix(iInds,jInds)
2457
        SubMat = MATMUL(R,SubMat0)
2458
        Matrix(iInds,jInds) = SubMat
2459
        
2460
        SubMat0 = Matrix(jInds,iInds)
2461
        SubMat = MATMUL(SubMat0,TRANSPOSE(R))
2462
        Matrix(jInds,iInds) = SubMat
2463
      END DO
2464

2465
      Force0 = Vector(iInds)
2466
      Force = MATMUL(R,Force0)
2467
      Vector(iInds) = Force
2468
      
2469
    END DO
2470
!------------------------------------------------------------------------------
2471
  END SUBROUTINE RotateMatrix
2472
!------------------------------------------------------------------------------
2473
#endif
2474

2475

2476
!------------------------------------------------------------------------------
2477
!> Given the normal return the tangent directions. The 
2478
!> First tangent direction will always be on the xy-plane if
2479
!> also the normal is in the xy-plane.
2480
!------------------------------------------------------------------------------
2481
  SUBROUTINE TangentDirections( Normal,Tangent1,Tangent2 )
2482
!------------------------------------------------------------------------------
2483
   REAL(KIND=dp) :: Normal(3),Tangent1(3)
2484
   REAL(KIND=dp), OPTIONAL :: Tangent2(3)
2485
!------------------------------------------------------------------------------
2486
   REAL(KIND=dp) :: n1,n2,n3
2487
!------------------------------------------------------------------------------
2488
   n1 = ABS(Normal(1))
2489
   n2 = ABS(Normal(2))
2490
   n3 = ABS(Normal(3))
2491

2492
   IF( PRESENT( Tangent2 ) ) THEN   
2493
     IF ( n1 <= n3 .AND. n2 <= n3 ) THEN
2494
       Tangent1(1) =  0.0_dp
2495
       Tangent1(2) = -Normal(3)
2496
       Tangent1(3) =  Normal(2)
2497
     ELSE
2498
       Tangent1(1) = -Normal(2)
2499
       Tangent1(2) =  Normal(1)
2500
       Tangent1(3) =  0.0_dp
2501
     END IF
2502

2503
     Tangent1 = Tangent1 / SQRT(SUM(Tangent1**2))
2504
     Tangent2(1) = Normal(2)*Tangent1(3) - Normal(3)*Tangent1(2)
2505
     Tangent2(2) = Normal(3)*Tangent1(1) - Normal(1)*Tangent1(3)
2506
     Tangent2(3) = Normal(1)*Tangent1(2) - Normal(2)*Tangent1(1)
2507
     Tangent2 = Tangent2 / SQRT(SUM(Tangent2**2))
2508
   ELSE
2509
     ! This is a 2D tangent only
2510
     Tangent1(1) = Normal(2)
2511
     Tangent1(2) = -Normal(1)
2512
     Tangent1(3) = 0.0_dp
2513
     Tangent1 = Tangent1 / SQRT(SUM(Tangent1**2))
2514
   END IF
2515
     
2516
!------------------------------------------------------------------------------
2517
 END SUBROUTINE TangentDirections
2518
!------------------------------------------------------------------------------
2519

2520

2521
!------------------------------------------------------------------------------
2522
!> Given one or two tangents return the normal direction.
2523
!> It is assumed that if one tangent is given the normal is in xy-plane otherwise
2524
!> in 3D. Note that the sign of normal vector is not unique.
2525
!------------------------------------------------------------------------------
2526
  FUNCTION NormalDirection( Tangent1,Tangent2 ) RESULT ( Normal ) 
2527
!------------------------------------------------------------------------------
2528
   REAL(KIND=dp) :: Normal(3),Tangent1(3)
2529
   REAL(KIND=dp), OPTIONAL :: Tangent2(3)
2530
!------------------------------------------------------------------------------
2531
   REAL(KIND=dp) :: t1(3),t2(3)
2532
!------------------------------------------------------------------------------
2533

2534
   IF(PRESENT(Tangent2)) THEN
2535
     Normal = CrossProduct(Tangent1,Tangent2)
2536
   ELSE
2537
     Normal(1) = Tangent1(2)
2538
     Normal(2) = -Tangent1(1)
2539
   END IF
2540
   Normal = Normal/SQRT(SUM(Normal**2))       
2541
!------------------------------------------------------------------------------
2542
 END FUNCTION NormalDirection
2543
!------------------------------------------------------------------------------
2544

2545

2546
!------------------------------------------------------------------------------
2547
!>    Integrates a user-defined function over the specified bulk elements.
2548
!> \deprecated Is this used for something?
2549
!------------------------------------------------------------------------------
2550
   FUNCTION VolumeIntegrate( Model, ElementList, IntegrandFunctionName ) &
2551
       RESULT(Integral)
2552
!------------------------------------------------------------------------------
2553
   TYPE(Model_t) :: Model                     !< All model information (mesh, materials, BCs, etc...)
2554
   INTEGER, DIMENSION(:) :: ElementList       !< List of elements that belong to the integration volume
2555
   CHARACTER(LEN=*) :: IntegrandFunctionName  !< Name the function has in the .sif file
2556
   REAL(KIND=dp) :: Integral                  !< The value of the volume integral
2557
!------------------------------------------------------------------------------
2558
     INTEGER :: n
2559
     TYPE(Element_t), POINTER :: CurrentElement
2560
     INTEGER, POINTER :: NodeIndexes(:)
2561
     TYPE(Nodes_t)   :: ElementNodes
2562

2563
     TYPE(GaussIntegrationPoints_t), TARGET :: IntegStuff
2564
     REAL(KIND=dp), DIMENSION(:), POINTER :: &
2565
         U_Integ,V_Integ,W_Integ,S_Integ
2566

2567
     REAL(KIND=dp), DIMENSION(Model % MaxElementNodes) :: IntegrandFunction
2568
     REAL(KIND=dp) :: s,ug,vg,wg
2569
     REAL(KIND=dp) :: Basis(Model % MaxElementNodes)
2570
     REAL(KIND=dp) :: dBasisdx(Model % MaxElementNodes,3),SqrtElementMetric
2571
     REAL(KIND=dp) :: IntegrandAtGPt, dV
2572
     INTEGER :: N_Integ, t, tg, i, istat
2573
     LOGICAL :: stat
2574

2575
! Need MaxElementNodes only in allocation
2576
     n = Model % MaxElementNodes
2577
     ALLOCATE( ElementNodes % x( n ),   &
2578
               ElementNodes % y( n ),   &
2579
               ElementNodes % z( n ), STAT=istat )
2580
     IF( istat /= 0 ) THEN
2581
       CALL Fatal('VolumeIntegrate','Allocation error for ElementNodes')
2582
     END IF
2583

2584
     Integral = 0.0d0
2585

2586
! Loop over all elements in the list
2587
     DO i=1,SIZE(ElementList)
2588

2589
       t = ElementList(i)
2590

2591
       IF ( t < 1 .OR. t > Model % NumberOfBulkElements ) THEN
2592
! do something
2593
       END IF
2594

2595
       CurrentElement => Model % Elements(t)
2596
       n = CurrentElement % TYPE % NumberOfNodes
2597
       NodeIndexes => CurrentElement % NodeIndexes
2598

2599
!------------------------------------------------------------------------------
2600
! Get element nodal coordinates
2601
!------------------------------------------------------------------------------
2602
       ElementNodes % x(1:n) = Model % Nodes % x(NodeIndexes)
2603
       ElementNodes % y(1:n) = Model % Nodes % y(NodeIndexes)
2604
       ElementNodes % z(1:n) = Model % Nodes % z(NodeIndexes)
2605

2606
! Read this from Simulation block in the .sif file
2607
       IntegrandFunction(1:n) = ListGetReal( Model % Simulation, &
2608
           IntegrandFunctionName, n, NodeIndexes )
2609

2610
!------------------------------------------------------------------------------
2611
!    Gauss integration stuff
2612
!------------------------------------------------------------------------------
2613
       IntegStuff = GaussPoints( CurrentElement )
2614
       U_Integ => IntegStuff % u
2615
       V_Integ => IntegStuff % v
2616
       W_Integ => IntegStuff % w
2617
       S_Integ => IntegStuff % s
2618
       N_Integ =  IntegStuff % n
2619

2620
!------------------------------------------------------------------------------
2621
! Loop over Gauss integration points
2622
!------------------------------------------------------------------------------
2623
       DO tg=1,N_Integ
2624

2625
         ug = U_Integ(tg)
2626
         vg = V_Integ(tg)
2627
         wg = W_Integ(tg)
2628

2629
!------------------------------------------------------------------------------
2630
! Need SqrtElementMetric and Basis at the integration point
2631
!------------------------------------------------------------------------------
2632
         stat = ElementInfo( CurrentElement,ElementNodes,ug,vg,wg, &
2633
             SqrtElementMetric,Basis,dBasisdx )
2634

2635
         s = SqrtElementMetric * S_Integ(tg)
2636

2637
! Calculate the function to be integrated at the Gauss point
2638
         IntegrandAtGPt = SUM( IntegrandFunction(1:n) * Basis )
2639

2640
! Use general coordinate system for dV
2641
         dV = CoordinateSqrtMetric( SUM( ElementNodes % x(1:n) * Basis), &
2642
             SUM( ElementNodes % y(1:n) * Basis), &
2643
             SUM( ElementNodes % z(1:n) * Basis) )
2644

2645
         Integral = Integral + s*IntegrandAtGPt*dV
2646

2647
       END DO! of the Gauss integration points
2648

2649
     END DO! of the bulk elements
2650

2651
     DEALLOCATE( ElementNodes % x, &
2652
         ElementNodes % y, &
2653
         ElementNodes % z )
2654
!------------------------------------------------------------------------------
2655
   END FUNCTION VolumeIntegrate
2656
!------------------------------------------------------------------------------
2657

2658

2659

2660
!------------------------------------------------------------------------------
2661
!>    Integrates the normal component of a user-defined vector function
2662
!>    over the specified boundary elements. 
2663
!>  \deprecated Is this used for something?
2664
!------------------------------------------------------------------------------
2665
   FUNCTION FluxIntegrate( Model, ElementList, IntegrandFunctionName ) &
2666
!------------------------------------------------------------------------------
2667
       RESULT(Integral)
2668
!------------------------------------------------------------------------------
2669
   TYPE(Model_t) :: Model                     !< All model information (mesh, materials, BCs, etc...)
2670
   INTEGER, DIMENSION(:) :: ElementList       !< List of elements that belong to the integration boundary
2671
   CHARACTER(LEN=*) :: IntegrandFunctionName  !< Name the function has in the .sif file
2672
   REAL(KIND=dp) :: Integral                  !< The value of the boundary integral
2673
!------------------------------------------------------------------------------
2674
     INTEGER :: n
2675
     TYPE(Element_t), POINTER :: CurrentElement
2676
     INTEGER, POINTER :: NodeIndexes(:)
2677
     TYPE(Nodes_t)    :: ElementNodes
2678

2679
     TYPE(GaussIntegrationPoints_t), TARGET :: IntegStuff
2680
     REAL(KIND=dp), DIMENSION(:), POINTER :: &
2681
         U_Integ,V_Integ,W_Integ,S_Integ
2682

2683
     REAL(KIND=dp), DIMENSION(Model % MaxElementNodes,3) :: IntegrandFunction
2684
!     REAL(KIND=dp), POINTER :: IntegrandFunction(:,:)
2685
     CHARACTER(LEN=2) :: Component
2686
     CHARACTER(:), ALLOCATABLE :: IntegrandFunctionComponent
2687
     REAL(KIND=dp) :: s,ug,vg,wg
2688
     REAL(KIND=dp) :: Basis(Model % MaxElementNodes)
2689
     REAL(KIND=dp) :: dBasisdx(Model % MaxElementNodes,3),SqrtElementMetric
2690
     REAL(KIND=dp) :: Normal(3)
2691
     REAL(KIND=dp) :: IntegrandAtGPt(3), FluxAtGPt, dS
2692
     INTEGER :: N_Integ, t, tg, i, j, DIM, istat
2693
     LOGICAL :: stat
2694

2695
     DIM = CoordinateSystemDimension()
2696

2697
! Need MaxElementNodes only in allocation
2698
     n = Model % MaxElementNodes
2699
     ALLOCATE( ElementNodes % x( n ),   &
2700
               ElementNodes % y( n ),   &
2701
               ElementNodes % z( n ), STAT=istat )
2702
     IF( istat /= 0 ) THEN
2703
       CALL Fatal('FluxIntegrate','Allocation error for ElementNodes')
2704
     END IF
2705

2706
     Integral = 0.0d0
2707

2708
! Loop over all elements in the list
2709
     DO i=1,SIZE(ElementList)
2710

2711
       t = ElementList(i)
2712

2713
       IF ( t < 1 .OR. t > Model % NumberOfBulkElements ) THEN
2714
! do something
2715
       END IF
2716

2717
       CurrentElement => Model % Elements(t)
2718
       n = CurrentElement % TYPE % NumberOfNodes
2719
       NodeIndexes => CurrentElement % NodeIndexes
2720

2721
       Model % CurrentElement => Model % Elements(t)
2722
!------------------------------------------------------------------------------
2723
! Get element nodal coordinates
2724
!------------------------------------------------------------------------------
2725
       ElementNodes % x(1:n) = Model % Nodes % x(NodeIndexes)
2726
       ElementNodes % y(1:n) = Model % Nodes % y(NodeIndexes)
2727
       ElementNodes % z(1:n) = Model % Nodes % z(NodeIndexes)
2728

2729
! Read the integrand from Simulation block in the .sif file
2730
! It is assumed to be a contravariant vector, but
2731
! ListGetRealArray doesn t exist, so we READ it component by component
2732
! naming them with suffixes " 1" etc.
2733
       DO j=1,DIM
2734
         IntegrandFunctionComponent = TRIM(IntegrandFunctionName)//' '//I2S(j)
2735
         IntegrandFunction(1:n,j) = ListGetReal( Model % Simulation, &
2736
             IntegrandFunctionComponent, n, NodeIndexes )
2737
       END DO
2738

2739
!------------------------------------------------------------------------------
2740
!    Gauss integration stuff
2741
!------------------------------------------------------------------------------
2742
       IntegStuff = GaussPoints( CurrentElement )
2743
       U_Integ => IntegStuff % u
2744
       V_Integ => IntegStuff % v
2745
       W_Integ => IntegStuff % w
2746
       S_Integ => IntegStuff % s
2747
       N_Integ =  IntegStuff % n
2748

2749
!------------------------------------------------------------------------------
2750
! Loop over Gauss integration points
2751
!------------------------------------------------------------------------------
2752
       DO tg=1,N_Integ
2753

2754
         ug = U_Integ(tg)
2755
         vg = V_Integ(tg)
2756
         wg = W_Integ(tg)
2757

2758
!------------------------------------------------------------------------------
2759
! Need SqrtElementMetric and Basis at the integration point
2760
!------------------------------------------------------------------------------
2761
         stat = ElementInfo( CurrentElement,ElementNodes,ug,vg,wg, &
2762
             SqrtElementMetric,Basis,dBasisdx )
2763

2764
! If we want to allow covariant integrand vectors given
2765
! we would also need Metric, and read it as follows...
2766
!     IF ( CurrentCoordinateSystem() /= Cartesian ) THEN
2767
!       x = SUM( nodes % x(1:n) * Basis )
2768
!       y = SUM( nodes % y(1:n) * Basis )
2769
!       z = SUM( nodes % z(1:n) * Basis )
2770
!     END IF
2771
!
2772
!     CALL CoordinateSystemInfo( Metric,SqrtMetric,Symb,dSymb,x,y,z )
2773
!     s = SqrtMetric * SqrtElementMetric * S_Integ(t)
2774
! And also FluxAtGpt...
2775
! Now we won t need Metric, so get SqrtMetric later...
2776

2777
         s = SqrtElementMetric * S_Integ(tg)
2778

2779
! Get normal to the element boundary at the integration point
2780
! N.B. NormalVector returns covariant normal vector
2781
         Normal = NormalVector( CurrentElement,ElementNodes,ug,vg,.TRUE. )
2782

2783
! Calculate the contravariant vector function to be integrated
2784
! at the Gauss point
2785
         DO j=1,DIM
2786
           IntegrandAtGPt(j) = SUM( IntegrandFunction(1:n,j) * Basis )
2787
         END DO
2788
! Calculate the normal component of the vector function
2789
         FluxAtGPt = SUM( IntegrandAtGPt * Normal )
2790

2791
! Use general coordinate system for dS
2792
! Would be included in s by SqrtMetric
2793
         dS = CoordinateSqrtMetric( SUM( ElementNodes % x(1:n) * Basis), &
2794
             SUM( ElementNodes % y(1:n) * Basis), &
2795
             SUM( ElementNodes % z(1:n) * Basis) )
2796

2797
         Integral = Integral + s*FluxAtGPt*dS
2798

2799
       END DO! of the Gauss integration points
2800

2801
     END DO! of the boundary elements
2802

2803
     DEALLOCATE( ElementNodes % x, &
2804
         ElementNodes % y, &
2805
         ElementNodes % z )
2806
!------------------------------------------------------------------------------
2807
   END FUNCTION FluxIntegrate
2808
!------------------------------------------------------------------------------
2809

2810

2811

2812
!------------------------------------------------------------------------------
2813
!>    Integrates A user-defined vector function 
2814
!>    over the specified boundary elements.
2815
!> \deprecated Is this used for something?
2816
!------------------------------------------------------------------------------
2817
   FUNCTION SurfaceIntegrate( Model, ElementList, IntegrandFunctionName ) &
2818
       RESULT(Integral)
2819
!------------------------------------------------------------------------------
2820
   TYPE(Model_t) :: Model                     !< All model information (mesh, materials, BCs, etc...)
2821
   INTEGER, DIMENSION(:) :: ElementList       !< List of elements that belong to the integration boundary
2822
   CHARACTER(LEN=*) :: IntegrandFunctionName  !< Name the function has in the .sif file
2823
   REAL(KIND=dp), DIMENSION(3) :: Integral    !< The vector value of the integral
2824
!------------------------------------------------------------------------------
2825
     INTEGER :: n
2826
     TYPE(Element_t), POINTER :: CurrentElement
2827
     INTEGER, POINTER :: NodeIndexes(:)
2828
     TYPE(Nodes_t)    :: ElementNodes
2829

2830
     TYPE(GaussIntegrationPoints_t), TARGET :: IntegStuff
2831
     REAL(KIND=dp), DIMENSION(:), POINTER :: &
2832
         U_Integ,V_Integ,W_Integ,S_Integ
2833

2834
     REAL(KIND=dp), DIMENSION(Model % MaxElementNodes,3) :: IntegrandFunction
2835
!     REAL(KIND=dp), POINTER :: IntegrandFunction(:,:)
2836
     CHARACTER(LEN=2) :: Component
2837
     CHARACTER(:), ALLOCATABLE :: IntegrandFunctionComponent
2838
     REAL(KIND=dp) :: s,ug,vg,wg
2839
     REAL(KIND=dp) :: Basis(Model % MaxElementNodes)
2840
     REAL(KIND=dp) :: dBasisdx(Model % MaxElementNodes,3),SqrtElementMetric
2841
!     REAL(KIND=dp) :: Normal(3)
2842
     REAL(KIND=dp) :: IntegrandAtGPt(3), dS
2843
     INTEGER :: N_Integ, t, tg, i, j, DIM, istat
2844
     LOGICAL :: stat
2845

2846
     DIM = CoordinateSystemDimension()
2847

2848
! Need MaxElementNodes only in allocation
2849
     n = Model % MaxElementNodes
2850
     ALLOCATE( ElementNodes % x( n ),   &
2851
               ElementNodes % y( n ),   &
2852
               ElementNodes % z( n ), STAT=istat )
2853
     IF( istat /= 0 ) THEN
2854
       CALL Fatal('SurfaceIntegrate','Allocation error for ElementNodes')
2855
     END IF
2856

2857
     Integral = 0.0d0
2858

2859
! Loop over all elements in the list
2860
     DO i=1,SIZE(ElementList)
2861

2862
       t = ElementList(i)
2863

2864
       IF ( t < 1 .OR. t > Model % NumberOfBulkElements ) THEN
2865
! do something
2866
       END IF
2867

2868
       CurrentElement => Model % Elements(t)
2869
       Model % CurrentElement => CurrentElement
2870
       n = CurrentElement % TYPE % NumberOfNodes
2871
       NodeIndexes => CurrentElement % NodeIndexes
2872

2873
!------------------------------------------------------------------------------
2874
! Get element nodal coordinates
2875
!------------------------------------------------------------------------------
2876
       ElementNodes % x(1:n) = Model % Nodes % x(NodeIndexes)
2877
       ElementNodes % y(1:n) = Model % Nodes % y(NodeIndexes)
2878
       ElementNodes % z(1:n) = Model % Nodes % z(NodeIndexes)
2879

2880
! Read the integrand from Simulation block in the .sif file
2881
! It is assumed to be a contravariant vector, but
2882
! ListGetRealArray doesn t exist, so we READ it component by component
2883
! naming them with suffixes " 1" etc.
2884
       DO j=1,DIM
2885
         IntegrandFunctionComponent = TRIM(IntegrandFunctionName)//' '//I2S(j)
2886
         IntegrandFunction(1:n,j) = ListGetReal( Model % Simulation, &
2887
             IntegrandFunctionComponent, n, NodeIndexes )
2888
       END DO
2889

2890
!------------------------------------------------------------------------------
2891
!    Gauss integration stuff
2892
!------------------------------------------------------------------------------
2893
       IntegStuff = GaussPoints( CurrentElement )
2894
       U_Integ => IntegStuff % u
2895
       V_Integ => IntegStuff % v
2896
       W_Integ => IntegStuff % w
2897
       S_Integ => IntegStuff % s
2898
       N_Integ =  IntegStuff % n
2899

2900
!------------------------------------------------------------------------------
2901
! Loop over Gauss integration points
2902
!------------------------------------------------------------------------------
2903
       DO tg=1,N_Integ
2904

2905
         ug = U_Integ(tg)
2906
         vg = V_Integ(tg)
2907
         wg = W_Integ(tg)
2908

2909
!------------------------------------------------------------------------------
2910
! Need SqrtElementMetric and Basis at the integration point
2911
!------------------------------------------------------------------------------
2912
         stat = ElementInfo( CurrentElement,ElementNodes,ug,vg,wg, &
2913
             SqrtElementMetric,Basis,dBasisdx )
2914

2915
! If we want to allow covariant integrand vectors given
2916
! we would also need Metric, and read it as follows...
2917
!     IF ( CurrentCoordinateSystem() /= Cartesian ) THEN
2918
!       x = SUM( nodes % x(1:n) * Basis )
2919
!       y = SUM( nodes % y(1:n) * Basis )
2920
!       z = SUM( nodes % z(1:n) * Basis )
2921
!     END IF
2922
!
2923
!     CALL CoordinateSystemInfo( Metric,SqrtMetric,Symb,dSymb,x,y,z )
2924
!     s = SqrtMetric * SqrtElementMetric * S_Integ(t)
2925
! Now we won t need Metric, so get SqrtMetric later...
2926

2927
         s = SqrtElementMetric * S_Integ(tg)
2928

2929
! If you need normal directly at the integration point
2930
!         Normal = NormalVector( CurrentElement,ElementNodes,ug,vg,.TRUE. )
2931

2932
! Calculate the contravariant vector function to be integrated
2933
! at the Gauss point
2934
         DO j=1,DIM
2935
           IntegrandAtGPt(j) = SUM( IntegrandFunction(1:n,j) * Basis )
2936
         END DO
2937

2938
! Use general coordinate system for dS
2939
! Would be included in s by SqrtMetric
2940
         dS = CoordinateSqrtMetric( SUM( ElementNodes % x(1:n) * Basis), &
2941
             SUM( ElementNodes % y(1:n) * Basis), &
2942
             SUM( ElementNodes % z(1:n) * Basis) )
2943

2944
         DO j=1,DIM
2945
           Integral(j) = Integral(j) + s*IntegrandAtGPt(j)*dS
2946
         END DO
2947

2948
       END DO! of the Gauss integration points
2949

2950
     END DO! of the boundary elements
2951

2952
     DEALLOCATE( ElementNodes % x, &
2953
         ElementNodes % y, &
2954
         ElementNodes % z )
2955
!------------------------------------------------------------------------------
2956
   END FUNCTION SurfaceIntegrate
2957
!------------------------------------------------------------------------------
2958

2959

2960

2961
!------------------------------------------------------------------------------
2962
!>    Integrates the normal component of a user-defined vector function
2963
!>    over a specified line element.
2964
!> \deprecated Is this used for something?
2965
!------------------------------------------------------------------------------
2966
   FUNCTION LineIntegrate( Model, LineElement, LineElementNodes, &
2967
       IntegrandFunctionName, QuadrantTreeExists, RootQuadrant ) &
2968
       RESULT(Integral)
2969
!------------------------------------------------------------------------------
2970
!
2971
!  ARGUMENTS:
2972
!
2973
!  TYPE(Model_t), POINTER :: Model,  
2974
!     INPUT: All model information (mesh, materials, BCs, etc...)
2975
!
2976
!  TYPE(Element_t) :: LineElement
2977
!     INPUT: Line element that belongs to the line of integration
2978
!
2979
!  REAL(KIND=dp), DIMENSION(LineElement % Type % NumberOfNodes,3) ::
2980
!     LineElementNodes
2981
!     INPUT: List of nodal point coordinates
2982
!
2983
!  CHARACTER(LEN=*) :: IntegrandFunctionName
2984
!     INPUT: Name the function has in the .sif file or somewhere else
2985
!
2986
!  LOGICAL :: QuadrantTreeExists
2987
!     INPUT: QuadrantTree has been built, use it in element search
2988
!
2989
!  TYPE(Quadrant_t), POINTER :: RootQuadrant
2990
!     OUTPUT: Quadrant tree structure root
2991
!
2992
!  FUNCTION RETURN VALUE:
2993
!    REAL(KIND=dp) :: Integral
2994
!     The value of the flux integral
2995
!      
2996
!------------------------------------------------------------------------------
2997
   TYPE(Model_t) :: Model
2998
! TARGET only for CurrentElement
2999
   TYPE(Element_t), TARGET :: LineElement
3000
   REAL(KIND=dp), DIMENSION(:,:), TARGET CONTIG :: LineElementNodes
3001
   CHARACTER(LEN=*) :: IntegrandFunctionName
3002
   REAL(KIND=dp) :: Integral
3003
   LOGICAL :: QuadrantTreeExists
3004
   TYPE(Quadrant_t), POINTER :: RootQuadrant
3005
!------------------------------------------------------------------------------
3006
     INTEGER :: n
3007
! Only one element at a time, need CurrentElement only for NormalVector!
3008
     TYPE(Element_t), POINTER :: CurrentElement
3009
! LineElement nodes don t belong to global node structure
3010
! Need the structure for the function calls
3011
     TYPE(Nodes_t) :: ElementNodes
3012

3013
     TYPE(GaussIntegrationPoints_t), TARGET :: IntegStuff
3014
     REAL(KIND=dp), DIMENSION(:), POINTER :: &
3015
         U_Integ,V_Integ,W_Integ,S_Integ
3016

3017
! IntegrandFunction at the bulk element nodal points
3018
     REAL(KIND=dp), DIMENSION(Model % MaxElementNodes,3) :: IntegrandAtNodes
3019
! IntegrandFunction at the Gauss points
3020
     REAL(KIND=dp), DIMENSION(LineElement % TYPE % GaussPoints,3) :: IntegrandFunction
3021
     CHARACTER(LEN=2) :: Component
3022
     CHARACTER(:), ALLOCATABLE :: IntegrandFunctionComponent
3023
     REAL(KIND=dp) :: s,ug,vg,wg
3024
     REAL(KIND=dp) :: Basis(LineElement % TYPE % NumberOfNodes)
3025
     REAL(KIND=dp) :: dBasisdx(LineElement % TYPE % NumberOfNodes,3),SqrtElementMetric
3026
     REAL(KIND=dp) :: Normal(3)
3027
! IntegrandFunction already at GPts
3028
     REAL(KIND=dp) :: FluxAtGPt, dS
3029
     INTEGER :: N_Integ, t, tg, i, j, DIM
3030
     LOGICAL :: stat
3031
! Search for the bulk element each Gauss point belongs to
3032
     TYPE(Element_t), POINTER :: BulkElement
3033
     TYPE(Nodes_t) :: BulkElementNodes
3034
     INTEGER, POINTER :: NodeIndexes(:)
3035
     REAL(KIND=dp), DIMENSION(3) :: LocalCoordsInBulkElem
3036
     REAL(KIND=dp), DIMENSION(3) :: Point
3037
     INTEGER :: nBulk, maxlevel=10, k, Quadrant
3038
     TYPE(Quadrant_t), POINTER :: LeafQuadrant
3039

3040
     DIM = CoordinateSystemDimension()
3041

3042
     n = LineElement % TYPE % NumberOfNodes
3043

3044
     Integral = 0.0d0
3045

3046
!------------------------------------------------------------------------------
3047
! Get element nodal coordinates
3048
!------------------------------------------------------------------------------
3049
! Move from LineElementNodes to Nodes_t structure
3050
     ElementNodes % x => LineElementNodes(1:n,1)
3051
     ElementNodes % y => LineElementNodes(1:n,2)
3052
     ElementNodes % z => LineElementNodes(1:n,3)
3053

3054
!------------------------------------------------------------------------------
3055
!    Gauss integration stuff
3056
!------------------------------------------------------------------------------
3057
     IntegStuff = GaussPoints( LineElement )
3058
     U_Integ => IntegStuff % u
3059
     V_Integ => IntegStuff % v
3060
     W_Integ => IntegStuff % w
3061
     S_Integ => IntegStuff % s
3062
     N_Integ =  IntegStuff % n
3063

3064
!------------------------------------------------------------------------------
3065
! Loop over Gauss integration points
3066
!------------------------------------------------------------------------------
3067
     DO tg=1,N_Integ
3068

3069
       ug = U_Integ(tg)
3070
       vg = V_Integ(tg)
3071
       wg = W_Integ(tg)
3072

3073
!------------------------------------------------------------------------------
3074
! Need SqrtElementMetric and Basis at the integration point
3075
!------------------------------------------------------------------------------
3076
       stat = ElementInfo( LineElement,ElementNodes,ug,vg,wg, &
3077
           SqrtElementMetric,Basis,dBasisdx )
3078

3079
! If we want to allow covariant integrand vectors given... (see FluxIntegrate)
3080

3081
       s = SqrtElementMetric * S_Integ(tg)
3082

3083
! Find in which bulk element the Gauss point belongs to
3084
       IF ( QuadrantTreeExists ) THEN
3085
! Find the last existing quadrant that the point belongs to
3086
         Point = [ SUM( ElementNodes % x(1:n) * Basis), &
3087
             SUM( ElementNodes % y(1:n) * Basis), &
3088
             SUM( ElementNodes % z(1:n) * Basis) ]
3089
!         PRINT*,'Point:', Point
3090
         CALL FindLeafElements(Point, DIM, RootQuadrant, LeafQuadrant)         
3091
!         PRINT*,'Elems in LeafQuadrant',LeafQuadrant % NElemsInQuadrant
3092
! Go through the bulk elements in the last ChildQuadrant only
3093
         nBulk = Model % MaxElementNodes
3094
         ALLOCATE( BulkElementNodes % x( nBulk ),   &
3095
             BulkElementNodes % y( nBulk ),   &
3096
             BulkElementNodes % z( nBulk ) )
3097
!         PRINT*,'Elements:', LeafQuadrant % Elements
3098
         DO k=1, LeafQuadrant % NElemsInQuadrant
3099
           BulkElement => Model % Elements( &
3100
               LeafQuadrant % Elements(k) )
3101
           nBulk = BulkElement % TYPE % NumberOfNodes
3102
           NodeIndexes => BulkElement % NodeIndexes
3103
           BulkElementNodes % x(1:nBulk) = Model % Nodes % x(NodeIndexes)
3104
           BulkElementNodes % y(1:nBulk) = Model % Nodes % y(NodeIndexes)
3105
           BulkElementNodes % z(1:nBulk) = Model % Nodes % z(NodeIndexes)
3106
           IF ( PointInElement( BulkElement,BulkElementNodes, &
3107
               [ SUM( ElementNodes % x(1:n) * Basis), &
3108
               SUM( ElementNodes % y(1:n) * Basis), &
3109
               SUM( ElementNodes % z(1:n) * Basis) ], &
3110
               LocalCoordsInBulkElem) ) EXIT
3111
         END DO
3112
!         PRINT*,'Point in Element: ', LeafQuadrant % Elements(k)
3113
        ELSE
3114
! Go through all BulkElements
3115
! Need MaxElementNodes only in allocation
3116
         nBulk = Model % MaxElementNodes
3117
         ALLOCATE( BulkElementNodes % x( nBulk ),   &
3118
             BulkElementNodes % y( nBulk ),   &
3119
             BulkElementNodes % z( nBulk ) )
3120
         DO k=1,Model % NumberOfBulkElements
3121
           BulkElement => Model % Elements(k)
3122
           nBulk = BulkElement % TYPE % NumberOfNodes
3123
           NodeIndexes => BulkElement % NodeIndexes
3124
           BulkElementNodes % x(1:nBulk) = Model % Nodes % x(NodeIndexes)
3125
           BulkElementNodes % y(1:nBulk) = Model % Nodes % y(NodeIndexes)
3126
           BulkElementNodes % z(1:nBulk) = Model % Nodes % z(NodeIndexes)
3127
           IF ( PointInElement(BulkElement,BulkElementNodes, &
3128
               [ SUM(ElementNodes % x(1:n) * Basis), &
3129
                 SUM(ElementNodes % y(1:n) * Basis), &
3130
                 SUM(ElementNodes % z(1:n) * Basis) ], &
3131
               LocalCoordsInBulkElem) ) EXIT
3132
         END DO
3133
!         PRINT*,'Point in Element: ', k
3134
       END IF
3135
! Calculate value of the function in the bulk element
3136
! Read the integrand from Simulation block in the .sif file
3137
! It is assumed to be a contravariant vector, but
3138
! ListGetRealArray doesn t exist, so we read it component by component
3139
! naming them with suffixes " 1" etc.
3140

3141
      DO j=1,DIM
3142
        WRITE (Component, '(" ",I1.1)') j
3143
        IntegrandFunctionComponent = IntegrandFunctionName(1: &
3144
            LEN_TRIM(IntegrandFunctionName))
3145
        IntegrandFunctionComponent(LEN_TRIM(IntegrandFunctionName)+1: &
3146
            LEN_TRIM(IntegrandFunctionName)+2) = Component
3147
        IntegrandAtNodes(1:nBulk,j) = ListGetReal( Model % Simulation, &
3148
            IntegrandFunctionComponent(1:LEN_TRIM(IntegrandFunctionComponent)), &
3149
            nBulk, NodeIndexes )
3150
      END DO
3151

3152
      DO j=1,DIM
3153
        IntegrandFunction(tg,j) = InterpolateInElement( BulkElement, &
3154
            IntegrandAtNodes(1:nBulk,j),LocalCoordsInBulkElem(1), &
3155
            LocalCoordsInBulkElem(2),LocalCoordsInBulkElem(3) )
3156
      END DO
3157

3158
      DEALLOCATE( BulkElementNodes % x, &
3159
          BulkElementNodes % y, &
3160
          BulkElementNodes % z )
3161

3162
! Get normal to the element boundary at the integration point
3163
! N.B. NormalVector returns covariant normal vector
3164
! NormalVector defined weirdly, doesn t accept LineElement as an argument,
3165
! but wants a pointer
3166
         CurrentElement => LineElement
3167
         Normal = NormalVector( CurrentElement,ElementNodes,ug,vg,.FALSE. )
3168
! Might be consistently in wrong direction, since no check
3169
!         Normal = NormalVector( CurrentElement,ElementNodes,ug,vg,.TRUE. )
3170

3171
! Contravariant vector function to be integrated is already
3172
! at the Gauss point
3173

3174
! Calculate the normal component of the vector function
3175
         FluxAtGPt = SUM( IntegrandFunction(tg,1:DIM) * Normal(1:DIM) )
3176

3177
! Use general coordinate system for dS
3178
! Would be included in s by SqrtMetric
3179
         dS = CoordinateSqrtMetric( SUM( ElementNodes % x(1:n) * Basis), &
3180
             SUM( ElementNodes % y(1:n) * Basis), &
3181
             SUM( ElementNodes % z(1:n) * Basis) )
3182

3183
         Integral = Integral + s*FluxAtGPt*dS
3184

3185
       END DO! of the Gauss integration points
3186

3187
!------------------------------------------------------------------------------
3188
   END FUNCTION LineIntegrate
3189
!------------------------------------------------------------------------------
3190

3191

3192
!------------------------------------------------------------------------------
3193
   FUNCTION ElementArea( Mesh,Element,N ) RESULT(A)
3194
!------------------------------------------------------------------------------
3195
     TYPE(Mesh_t), POINTER :: Mesh
3196
     INTEGER :: N
3197
     TYPE(Element_t) :: Element
3198
!------------------------------------------------------------------------------
3199

3200
     REAL(KIND=dp), TARGET :: NX(N),NY(N),NZ(N)
3201

3202
     REAL(KIND=dp) :: A,R1,R2,Z1,Z2,S,U,V,W,X,Y,Z
3203

3204
     TYPE(GaussIntegrationPoints_t), TARGET :: IntegStuff
3205
     INTEGER :: N_Integ,t
3206

3207
     REAL(KIND=dp) :: Metric(3,3),Symb(3,3,3),dSymb(3,3,3,3), &
3208
              SqrtMetric,SqrtElementMetric
3209

3210
     TYPE(Nodes_t) :: Nodes
3211

3212
     LOGICAL :: stat
3213

3214
     REAL(KIND=dp) :: Basis(n)
3215
     REAL(KIND=dp) :: dBasisdx(n,3)
3216

3217
     REAL(KIND=dp), DIMENSION(:), POINTER :: U_Integ,V_Integ,W_Integ,S_Integ
3218
!------------------------------------------------------------------------------
3219
 
3220
#if 0
3221
     IF ( ( CurrentCoordinateSystem() == AxisSymmetric .OR. &
3222
            CurrentCoordinateSystem() == CylindricSymmetric ) .AND. Element % TYPE % ELementCode / 100 == 2 ) THEN
3223
       R1 = Mesh % Nodes % x(Element % NodeIndexes(1))
3224
       R2 = Mesh % Nodes % x(Element % NodeIndexes(2))
3225

3226
       Z1 = Mesh % Nodes % y(Element % NodeIndexes(1))
3227
       Z2 = Mesh % Nodes % y(Element % NodeIndexes(2))
3228

3229
       A = PI*ABS(R1+R2)*SQRT((Z1-Z2)*(Z1-Z2)+(R1-R2)*(R1-R2))
3230
     ELSE 
3231
#endif
3232
       Nodes % x => NX
3233
       Nodes % y => NY
3234
       Nodes % z => NZ
3235

3236
       Nodes % x = Mesh % Nodes % x(Element % NodeIndexes)
3237
       Nodes % y = Mesh % Nodes % y(Element % NodeIndexes)
3238
       Nodes % z = Mesh % Nodes % z(Element % NodeIndexes)
3239

3240
       IntegStuff = GaussPoints( element )
3241
       U_Integ => IntegStuff % u
3242
       V_Integ => IntegStuff % v
3243
       W_Integ => IntegStuff % w
3244
       S_Integ => IntegStuff % s
3245
       N_Integ  = IntegStuff % n
3246
!
3247
!------------------------------------------------------------------------------
3248
!   Now we start integrating
3249
!------------------------------------------------------------------------------
3250
!
3251
       A = 0.0
3252
       DO t=1,N_Integ
3253
!
3254
!        Integration stuff
3255
!
3256
         u = U_Integ(t)
3257
         v = V_Integ(t)
3258
         w = W_Integ(t)
3259
!
3260
!------------------------------------------------------------------------------
3261
!        Basis function values & derivatives at the integration point
3262
!------------------------------------------------------------------------------
3263
         stat = ElementInfo( Element,Nodes,u,v,w,SqrtElementMetric, &
3264
                    Basis,dBasisdx )
3265
!------------------------------------------------------------------------------
3266
!        Coordinatesystem dependent info
3267
!------------------------------------------------------------------------------
3268
         IF ( CurrentCoordinateSystem() /= Cartesian ) THEN
3269
           X = SUM( Nodes % x(1:n)*Basis )
3270
           Y = SUM( Nodes % y(1:n)*Basis )
3271
           Z = SUM( Nodes % z(1:n)*Basis )
3272

3273
           SqrtMetric = CoordinateSqrtMetric( x,y,z )
3274
           A =  A + SqrtMetric * SqrtElementMetric * S_Integ(t)
3275
         ELSE
3276
           A =  A + SqrtElementMetric * S_Integ(t)
3277
         END IF
3278
       END DO
3279
#if 0
3280
     END IF
3281
#endif
3282

3283
   END FUNCTION ElementArea
3284
!------------------------------------------------------------------------------
3285

3286

3287
   !------------------------------------------------------------------------------
3288
   !> If element has two of the same indexes regard the element as degenerate.
3289
   !------------------------------------------------------------------------------
3290
   FUNCTION DegenerateElement( Element ) RESULT ( Stat ) 
3291
     TYPE(Element_t), POINTER :: Element
3292
     LOGICAL Stat
3293

3294
     INTEGER :: i,n
3295
     INTEGER, POINTER :: Indexes(:)
3296
     
3297
     Stat = .FALSE.
3298

3299
     n = Element % TYPE % NumberOfNodes
3300
     Indexes => Element % NodeIndexes
3301
     
3302
     DO i = 1, n
3303
       IF( ANY( Indexes(i+1:n) == Indexes(i) ) ) THEN
3304
         Stat = .TRUE.           
3305
         EXIT
3306
       END IF
3307
     END DO
3308
     
3309
   END FUNCTION DegenerateElement
3310

3311
   !------------------------------------------------------------------------------
3312
   !> Return the aspect ratio of an element 
3313
   !------------------------------------------------------------------------------
3314
   FUNCTION ElementAspectRatio(Model, Element ) RESULT ( AspectRatio ) 
3315
     IMPLICIT NONE
3316
     TYPE(Model_t) Model
3317
     TYPE(Element_t) :: Element
3318
     REAL(KIND=dp) :: AspectRatio
3319
     REAL(KIND=dp) :: CharLen(2)
3320

3321
     CharLen = ElementCharacteristicLengths(Model, Element)
3322
     IF (CharLen(1) <= 0) THEN
3323
       AspectRatio = HUGE(AspectRatio)
3324
     ELSE
3325
       AspectRatio = CharLen(2)/CharLen(1)
3326
     END IF
3327
   END FUNCTION ElementAspectRatio
3328

3329
   !------------------------------------------------------------------------------
3330
   !> Return the characteristic lengths of an element 
3331
   !------------------------------------------------------------------------------
3332
   FUNCTION ElementCharacteristicLengths(Model, Element ) RESULT ( Charlengths ) 
3333
     IMPLICIT NONE
3334
     TYPE(Model_t) :: Model
3335
     TYPE(Element_t) :: Element
3336
     REAL(KIND=dp) :: Charlengths(2)
3337
     REAL(KIND=dp) :: Dist
3338

3339
     TYPE(Nodes_t) :: en
3340
     INTEGER :: i,j,n
3341

3342
     INTEGER :: istat
3343
     
3344
     n = Element % TYPE % NumberOfNodes
3345

3346
     ALLOCATE( en % x( n ),   &
3347
               en % y( n ),   &
3348
               en % z( n ), STAT=istat )
3349

3350
     IF( istat /= 0 ) THEN
3351
       CALL Fatal('ElementCharacteristicLengths','Allocation error for ElementNodes')
3352
     END IF
3353

3354
     en % x(1:n) = Model % Nodes % x(Element % NodeIndexes)
3355
     en % y(1:n) = Model % Nodes % y(Element % NodeIndexes)
3356
     en % z(1:n) = Model % Nodes % z(Element % NodeIndexes)
3357
     
3358
     Charlengths = 0._dp
3359
     DO i = 1, n
3360
       DO j = 1, n
3361
         IF (i /= j) THEN
3362
           Dist = SQRT((en % x(i)-en % x(j))**2. + (en % y(i)-en % y(j))**2. + (en % z(i)-en % z(j))**2.)
3363
           IF (Dist < Charlengths(1)) THEN
3364
             Charlengths(1) = Dist
3365
           ELSE IF (Dist > Charlengths(2)) THEN
3366
             Charlengths(2) = Dist
3367
           END IF
3368
         END IF
3369
       END DO
3370
     END DO
3371
     
3372
   END FUNCTION ElementCharacteristicLengths
3373
   
3374
   !------------------------------------------------------------------------------
3375
   !> Return normal of degenerate Element 
3376
   !------------------------------------------------------------------------------
3377
   FUNCTION NormalOfDegenerateElement(Model, Element ) RESULT ( Normal ) 
3378
     IMPLICIT NONE
3379
     TYPE(Model_t) :: Model
3380
     TYPE(Element_t) :: Element
3381
     REAL(KIND=dp) :: a(3), b(3), c(3), Normal(3)
3382

3383
     TYPE(Nodes_t) :: en
3384
     INTEGER :: i,n
3385

3386
     INTEGER :: istat
3387
     
3388
     n = Element % TYPE % NumberOfNodes
3389

3390
     ALLOCATE( en % x( n ),   &
3391
               en % y( n ),   &
3392
               en % z( n ), STAT=istat )
3393

3394
     IF( istat /= 0 ) THEN
3395
       CALL Fatal('NormalOfDegenerateElement','Allocation error for ElementNodes')
3396
     END IF
3397

3398
     en % x(1:n) = Model % Nodes % x(Element % NodeIndexes)
3399
     en % y(1:n) = Model % Nodes % y(Element % NodeIndexes)
3400
     en % z(1:n) = Model % Nodes % z(Element % NodeIndexes)
3401

3402
     a = (/ en % x(1), en % y(1), en % z(1) /)
3403
     b = (/ en % x(2), en % y(2), en % z(2) /)
3404
     c = (/ en % x(n), en % y(n), en % z(n) /)
3405

3406
     Normal = crossproduct(a-b, a-c)
3407

3408
     Normal = Normal / SQRT(SUM(c**2))
3409
     
3410
   END FUNCTION NormalOfDegenerateElement
3411

3412

3413
!------------------------------------------------------------------------------
3414
   FUNCTION FindBoundaryEdgeIndex(Mesh,Boundary,nedge) RESULT(n)
3415
!------------------------------------------------------------------------------
3416
     IMPLICIT NONE
3417
     INTEGER :: n,nedge
3418
     TYPE(Mesh_t), POINTER :: Mesh
3419
     TYPE(Element_t) :: Boundary
3420
!------------------------------------------------------------------------------
3421
     INTEGER :: i,j,k,jb1,jb2,je1,je2
3422
     TYPE(Element_t), POINTER :: Parent, Edge, Face
3423
!------------------------------------------------------------------------------
3424
     n = 0
3425
     SELECT CASE(Boundary % TYPE % ElementCode / 100 )
3426
     CASE(1)
3427
       RETURN
3428
     CASE(2)
3429
       IF ( nedge==1 ) THEN
3430
         Parent => Boundary % BoundaryInfo % Left
3431
         IF ( .NOT. ASSOCIATED(Parent) ) &
3432
             Parent => Boundary % BoundaryInfo % Right
3433

3434
         jb1 = Boundary % NodeIndexes(1)
3435
         jb2 = Boundary % NodeIndexes(2)
3436
         DO i=1,Parent % TYPE % NumberOfEdges
3437
           Edge => Mesh % Edges(Parent % EdgeIndexes(i))
3438
           je1 = Edge % NodeIndexes(1)
3439
           je2 = Edge % NodeIndexes(2)
3440
           IF ( jb1==je1.AND.jb2==je2 .OR. jb1==je2.AND.jb2==je1) EXIT
3441
         END DO
3442
         n = Parent % EdgeIndexes(i)
3443
       END IF
3444
     CASE(3,4)
3445
       j = FindBoundaryFaceIndex(Mesh,Boundary)
3446
       Face => Mesh % Faces(j)
3447
       IF ( nedge>0.AND.nedge<=Face % TYPE % NumberOfEdges ) &
3448
           n = Face % EdgeIndexes(nedge) 
3449
     END SELECT
3450
!------------------------------------------------------------------------------
3451
   END FUNCTION FindBoundaryEdgeIndex
3452
!------------------------------------------------------------------------------
3453

3454

3455
!------------------------------------------------------------------------------
3456
   FUNCTION FindBoundaryFaceIndex(Mesh,Boundary) RESULT(n)
3457
!------------------------------------------------------------------------------
3458
     IMPLICIT NONE
3459
     INTEGER :: n
3460
     TYPE(Element_t) :: Boundary
3461
     TYPE(Mesh_t), POINTER :: Mesh
3462
!------------------------------------------------------------------------------
3463
     INTEGER :: i,j,k,m
3464
     TYPE(Element_t), POINTER :: Parent, Face
3465
!------------------------------------------------------------------------------
3466
     Parent => Boundary % BoundaryInfo % Left
3467
     IF ( .NOT. ASSOCIATED(Parent) ) &
3468
         Parent => Boundary % BoundaryInfo % Right
3469

3470
     DO i=1,Parent % TYPE % NumberOfFaces
3471
       Face => Mesh % Faces(Parent % FaceIndexes(i))
3472
       m = 0
3473
       DO j=1,Face % TYPE % NumberOfNodes
3474
         DO k=1,Boundary % TYPE % NumberOfNodes
3475
           IF ( Face % NodeIndexes(j)==Boundary % NodeIndexes(k)) m=m+1
3476
         END DO
3477
       END DO
3478
       IF ( m==Face % TYPE % NumberOfNodes) EXIT
3479
     END DO
3480
     n = Parent % FaceIndexes(i)
3481
!------------------------------------------------------------------------------
3482
   END FUNCTION FindBoundaryFaceIndex
3483
!------------------------------------------------------------------------------
3484

3485
   
3486
!-----------------------------------------------------------------------------   
3487
!> Given basis function values at surface element find the corresponding local
3488
!> coordinate in the parent element. 
3489
!------------------------------------------------------------------------------
3490
   SUBROUTINE FindParentUVW( Element, n, Parent, np, U, V, W, Basis ) 
3491
!------------------------------------------------------------------------------
3492
     IMPLICIT NONE
3493
     TYPE( Element_t ), POINTER :: Element
3494
     TYPE( Element_t ), POINTER :: Parent
3495
     INTEGER :: n, np
3496
     REAL( KIND=dp ) :: U, V, W, Basis(:)
3497
!------------------------------------------------------------------------------
3498
    INTEGER :: i, j, nParent, check 
3499
    REAL(KIND=dp) :: NodalParentU(n), NodalParentV(n), NodalParentW(n)
3500
!------------------------------------------------------------------------------
3501

3502
    Check = 0
3503

3504
    DO i = 1,n
3505
      DO j = 1,np
3506
        IF( Element % NodeIndexes(i) == Parent % NodeIndexes(j) ) THEN
3507
          Check = Check + 1
3508
          NodalParentU(i) = Parent % Type % NodeU(j)
3509
          NodalParentV(i) = Parent % Type % NodeV(j)
3510
          NodalParentW(i) = Parent % Type % NodeW(j)
3511
        END IF
3512
      END DO
3513
    END DO
3514

3515
    IF( Check /= n ) THEN
3516
      IF(n /= Element % TYPE % NumberOfNodes ) THEN
3517
        CALL Warn('FindParentUVW','Inconsistent size for "n"!')
3518
      END IF
3519
      IF(np /= Parent % TYPE % NumberOfNodes ) THEN
3520
        CALL Warn('FindParentUVW','Inconsistent size for "np"!')
3521
      END IF
3522
      CALL Fatal('FindParentUVW','Could not find all nodes in parent!') 
3523
    END IF
3524
    
3525
    U = SUM( Basis(1:n) * NodalParentU(1:n) )
3526
    V = SUM( Basis(1:n) * NodalParentV(1:n) )
3527
    W = SUM( Basis(1:n) * NodalParentW(1:n) )
3528
!------------------------------------------------------------------------------      
3529
  END SUBROUTINE FindParentUVW
3530
!------------------------------------------------------------------------------      
3531

3532

3533
!------------------------------------------------------------------------------------------   
3534
!> Routine for obtaining values at a given point with as many features as required in
3535
!> SavaScalars, SaveLine etc. solvers. Idea is that the same info is not needed in many
3536
!> places.
3537
!------------------------------------------------------------------------------      
3538
  SUBROUTINE EvaluateVariableAtGivenPoint(No,Values,Mesh,Var,Var2,Var3,Element,LocalCoord,&
3539
      LocalBasis,LocalNode,LocalDGNode,DoGrad,DoDiv,GotEigen,GotModes,GotEdge,Parent)
3540

3541
    INTEGER :: No
3542
    REAL(KIND=dp) :: Values(:)
3543
    TYPE(Mesh_t), POINTER :: Mesh
3544
    TYPE(Variable_t), POINTER :: Var
3545
    TYPE(Variable_t), POINTER, OPTIONAL :: Var2, Var3
3546
    TYPE(Element_t), POINTER, OPTIONAL :: Element
3547
    REAL(KIND=dp), OPTIONAL :: LocalCoord(3)
3548
    INTEGER, OPTIONAL :: LocalNode
3549
    INTEGER, OpTIONAL :: LocalDGNode
3550
    REAL(KIND=dp), OPTIONAL, TARGET :: LocalBasis(:)
3551
    LOGICAL, OPTIONAL :: DoGrad, DoDiv
3552
    LOGICAL, OPTIONAL :: GotEigen, GotModes, GotEdge
3553
    TYPE(Element_t), POINTER, OPTIONAL :: Parent
3554
    
3555
    LOGICAL :: Found, EdgeBasis, AVBasis, DgVar, IpVar, ElemVar, DoEigen, &
3556
        PiolaVersion, PElem, NeedDerBasis, Stat, UseGivenNode, IsGrad, IsDiv, &
3557
        IsEigen, IsModes
3558
    INTEGER :: i1,i2,ii,i,j,k,l,comps, n, n2, nd, np, NoEigenValues, NoConstraintModes, iMode
3559
    INTEGER, TARGET :: DGIndexes(27), Indexes(100), DofIndexes(100), NodeIndex(1)
3560
    INTEGER, POINTER :: pToIndexes(:)
3561
    REAL(KIND=dp) :: u,v,w,detJ
3562
    TYPE(Nodes_t), SAVE :: Nodes
3563
    INTEGER :: lr
3564
    REAL(KIND=dp), TARGET :: Basis(100),dBasisdx(100,3),NodeBasis(1)
3565
    REAL(KIND=dp) :: WBasis(54,3),RotWBasis(54,3)
3566
    REAL(KIND=dp), ALLOCATABLE, SAVE :: fdg(:), fip(:)
3567
    REAL(KIND=dp), POINTER :: pToBasis(:)
3568
    TYPE(Variable_t), POINTER :: pVar
3569
    TYPE(Element_t), POINTER :: Element2    
3570
    REAL(KIND=dp), POINTER :: rValues(:)
3571
    COMPLEX(KIND=dp), POINTER :: cValues(:)
3572

3573
    INTERFACE 
3574
      SUBROUTINE Ip2DgFieldInElement( Mesh, Parent, nip, fip, np, fdg )
3575
        USE Types
3576
        TYPE(Mesh_t), POINTER :: Mesh
3577
        TYPE(Element_t), POINTER :: Parent
3578
        INTEGER :: nip, np
3579
        REAL(KIND=dp) :: fip(:), fdg(:)
3580
      END SUBROUTINE Ip2DgFieldInElement
3581
    END INTERFACE
3582

3583
    IF(PRESENT(GotEdge)) GotEdge = .FALSE.
3584
    IF(PRESENT(GotEigen)) GotEIgen = .FALSE.
3585
    IF(PRESENT(GotModes)) GotModes = .FALSE.
3586
        
3587
    IF(.NOT. ASSOCIATED(Var)) RETURN
3588
    IF(.NOT. ASSOCIATED(Var % Values)) RETURN
3589

3590
    ! Check that a vector field was not given by components 
3591
    comps = 1
3592
    IF(PRESENT(Var2)) THEN
3593
      IF(ASSOCIATED(Var2)) THEN
3594
        comps = 2
3595
        IF( PRESENT(Var3)) THEN
3596
          IF(ASSOCIATED(Var3)) THEN
3597
            comps = 3
3598
          END IF
3599
        END IF
3600
      END IF
3601
    END IF
3602

3603
    IF(.NOT. PRESENT(Element) ) THEN
3604
      CALL Fatal('EvaluteVariableAtGivenPoint','This routine most often really likes to have Element provided!')
3605
    END IF
3606
    
3607
    EdgeBasis = ( Var % TYPE == variable_on_edges )
3608

3609
    DGVar = ( Var % TYPE == variable_on_nodes_on_elements ) 
3610
    IpVar = ( Var % TYPE == variable_on_gauss_points )
3611
    ElemVar = ( Var % TYPE == Variable_on_elements )       
3612

3613
    ! We do not need P-elements if the value is to be found in node since
3614
    ! the higher order p-basis does not have any effect there.
3615
    pElem = .FALSE.
3616
    IF(PRESENT(LocalCoord)) THEN
3617
!    IF(.NOT. PRESENT(LocalNode)) THEN
3618
      IF(ASSOCIATED(Var % Solver)) THEN
3619
        pElem = isActivePElement(Element, Var % Solver) 
3620
      END IF
3621
    END IF
3622
      
3623
    PiolaVersion = .FALSE.
3624
    n = Element % TYPE % NumberOfNodes                  
3625

3626
    ! Elemental variable is a special case which is constant on the element
3627
    ! This does not depend on any other things. 
3628
    IF( ElemVar ) THEN
3629
      l = Element % ElementIndex 
3630
      IF( SIZE( Var % Perm ) >= l ) THEN
3631
        l = Var % Perm(l)
3632
      END IF
3633
      IF( l > 0 ) THEN
3634
        IF( Var % Dofs > 1 ) THEN
3635
          DO ii=1,Var % Dofs
3636
            Values(No+ii) = Var % Values(Var%Dofs*(l-1)+ii)
3637
          END DO
3638
        ELSE
3639
          Values(No+1) = Var % Values(l)              
3640
          IF( comps >= 2 ) Values(No+2) = Var % Values(l)
3641
          IF( comps >= 3 ) Values(No+3) = Var % Values(l)
3642
        END IF
3643
      END IF
3644
      No = No + MAX( Var % Dofs, comps )       
3645
      RETURN
3646
    END IF
3647

3648
    ! We may need to compute dBasisdx if we want derivative or divergence at the point. 
3649
    NeedDerBasis = .FALSE.
3650
    IsGrad = .FALSE.
3651
    IsDiv = .FALSE.
3652
    IF(PRESENT(DoGrad)) IsGrad = DoGrad
3653
    IF(PRESENT(DoDiv)) IsDiv = DoDiv
3654
    NeedDerBasis = IsGrad .OR. IsDiv
3655
    
3656
    pToBasis => NULL()
3657
    pToIndexes => NULL()
3658
    
3659
    IsEigen = .FALSE.
3660
    IF( PRESENT( GotEigen ) ) THEN
3661
      IsEigen = ASSOCIATED( Var % EigenValues )
3662
      IF(IsEigen) NoEigenValues = SIZE( Var % EigenValues )
3663
      IF( comps > 1 .AND. IsEigen ) THEN
3664
        CALL Warn('EvaluteVariableAtGivenPoint','Eigenmode cannot be given in components!')
3665
        IsEigen = .FALSE.
3666
      END IF
3667
      GotEigen = IsEigen
3668
    END IF
3669

3670
    IsModes = .FALSE.
3671
    IF( PRESENT( GotModes ) ) THEN
3672
      IsModes = ASSOCIATED( Var % ConstraintModes )
3673
      IF(IsModes) NoConstraintModes = Var % NumberOfConstraintModes
3674
      IF( comps > 1 .AND. IsEigen ) THEN
3675
        CALL Fatal('EvaluteVariableAtGivenPoint','Constraint modes cannot be given in components!')
3676
        IsModes = .FALSE.
3677
      END IF
3678
      GotModes = IsModes
3679
    END IF
3680
    
3681
    ! Given node is the quickest way to estimate the values at nodes.
3682
    ! However, it only applies to nodal fields. 
3683
    NodeIndex(1) = 0
3684
    IF(.NOT. (EdgeBasis .OR. IpVar .OR. ElemVar .OR. NeedDerBasis .OR. pElem) ) THEN      
3685
      IF( DgVar ) THEN        
3686
        IF(PRESENT(LocalDGNode)) NodeIndex(1) = LocalDGNode
3687
        IF(NodeIndex(1) == 0 ) THEN
3688
          IF( PRESENT(LocalNode)) THEN
3689
            PToIndexes => PickDGIndexes(Element,Parent)
3690
            DO i=1, n
3691
              IF( Element % NodeIndexes(i) == LocalNode ) THEN
3692
                NodeIndex(1) = pToIndexes(i)
3693
                EXIT
3694
              END IF
3695
            END DO
3696
          END IF
3697
        END IF
3698
      ELSE IF( PRESENT(LocalNode) ) THEN
3699
        NodeIndex(1) = LocalNode
3700
      END IF
3701

3702
      IF( NodeIndex(1) > 0 ) THEN
3703
        pToIndexes => NodeIndex        
3704
        NodeBasis(1) = 1.0_dp
3705
        pToBasis => NodeBasis
3706
        n = 1
3707
        nd = 1        
3708
      ELSE IF( PRESENT( LocalBasis ) ) THEN
3709
        ! The 2nd quickest is to use existing local nodal basis functions
3710
        PtoBasis => LocalBasis
3711
        nd = n        
3712
        IF( DgVar ) THEN
3713
          PtoIndexes => PickDGIndexes(Element,Parent)
3714
        ELSE
3715
          PtoIndexes => Element % NodeIndexes 
3716
        END IF
3717
      END IF
3718
    END IF    
3719

3720
    ! We don't have basis to evaluate the the fields. Continue to find suitable basis. 
3721
    IF(ASSOCIATED(PtoIndexes) ) THEN
3722
      Element2 => Element
3723
      n2 = n
3724
    ELSE      
3725
      IF(.NOT. PRESENT(LocalCoord) ) THEN
3726
        CALL Fatal('EvaluteVariableAtGivenPoint',&
3727
            'No recipe to evaluate variable without local coordinates: '//TRIM(Var % Name))
3728
      END IF
3729

3730
      u = LocalCoord(1)
3731

3732
      IF( PRESENT(Parent)) THEN
3733
        Element2 => Parent
3734
        n2 = Element2 % TYPE % NumberOfNodes
3735
        IF(.NOT. PRESENT(LocalBasis)) THEN
3736
          CALL Fatal('EvaluateVariableAtGivenPoint','"LocalBasis" is needed when using "Parent" element!')
3737
        END IF
3738
        CALL FindParentUVW( Element, n, Parent, n2, U, V, W, LocalBasis )
3739
      ELSE
3740
        Element2 => Element
3741
        n2 = n
3742
        u = LocalCoord(1)
3743
        v = LocalCoord(2)
3744
        w = LocalCoord(3)
3745
      END IF
3746

3747
      IF( ASSOCIATED( Var % Solver ) ) THEN
3748
        nd = mGetElementDOFs( Indexes, Element2, Var % Solver)
3749
      ELSE
3750
        nd = mGetElementDOFs( Indexes, Element2 )
3751
      END IF          
3752

3753
      IF( EdgeBasis ) THEN
3754
        IF( ListGetLogical(Var % Solver % Values, 'Quadratic Approximation', Found) ) THEN
3755
          PiolaVersion = .TRUE.
3756
        ELSE
3757
          PiolaVersion = ListGetLogical(Var % Solver % Values,'Use Piola Transform', Found )   
3758
        END IF
3759
        np = n2 * Var % Solver % Def_Dofs(Element2 % type % ElementCode/100,Element2 % BodyId,1)
3760
        AVBasis = (np > 0 )
3761
        IF( PRESENT( GotEdge ) ) GotEdge = .TRUE.
3762
      END IF
3763
            
3764
      IF( EdgeBasis ) THEN
3765
        CALL CopyElementNodesFromMesh(Nodes,Mesh,n2,Element2 % NodeIndexes)
3766
        stat = ElementInfo( Element2, Nodes, u, v, w, &
3767
            detJ, Basis, dBasisdx,  EdgeBasis = WBasis, &
3768
            RotBasis = RotWBasis, USolver = Var % Solver)
3769
      ELSE          
3770
        IF( pElem ) THEN
3771
          ! The standard element of SaveLine, SaveScalars etc. is most likely standard nodal element.
3772
          ! If the user gives something else we may be trouble...
3773
          !CALL CopyElementNodesFromMesh(Nodes,Mesh,n2,Element2 % NodeIndexes)
3774
          PtoIndexes => Indexes
3775
          CALL CopyElementNodesFromMesh(Nodes,Mesh,nd,pToIndexes)
3776
          CALL ConvertToPReference(Element2 % TYPE % ElementCode,u,v,w)            
3777
          IF( NeedDerBasis ) THEN
3778
            stat = ElementInfo( Element2, Nodes, u, v, w, detJ, Basis, dBasisdx, &
3779
                USolver = Var % Solver )
3780
          ELSE
3781
            stat = ElementInfo( Element2, Nodes, u, v, w, detJ, Basis, &
3782
                USolver = Var % Solver )
3783
          END IF
3784
        ELSE
3785
          IF( NeedDerBasis ) THEN
3786
            CALL CopyElementNodesFromMesh(Nodes,Mesh,n2,Element2 % NodeIndexes)
3787
            stat = ElementInfo( Element2, Nodes, u, v, w, detJ, Basis, dBasisdx )
3788
          ELSE
3789
            CALL CopyElementNodesFromMesh(Nodes,Mesh,n2,Element2 % NodeIndexes)          
3790
            stat = ElementInfo( Element2, Nodes, u, v, w, detJ, Basis )
3791
          END IF
3792
        END IF
3793

3794
        IF( DgVar ) THEN
3795
          PtoIndexes => PickDGIndexes(Element2)
3796
        END IF
3797
      END IF
3798
      IF(.NOT. ASSOCIATED(pToIndexes)) PtoIndexes => Indexes
3799
      IF(.NOT. ASSOCIATED(pToBasis)) PtoBasis => Basis
3800
    END IF
3801
    
3802
    IF( EdgeBasis ) THEN
3803
      DofIndexes(1:nd) = Var % Perm(PToIndexes(1:nd))
3804
      IF( IsEigen ) THEN
3805
        DO iMode = 1, NoEigenValues
3806
          cValues => Var % EigenVectors(iMode,:)
3807
          DO j=1,3          
3808
            No = No + 1
3809
            IF( ALL(DofIndexes(np+1:nd) > 0 ) ) THEN            
3810
              Values(No) = SUM( WBasis(1:nd-np,j) * cValues(DofIndexes(np+1:nd)))
3811
            ELSE
3812
              Values(No) = 0.0_dp
3813
            END IF
3814
          END DO
3815
        END DO
3816
      ELSE IF( IsModes ) THEN
3817
        DO iMode = 1, NoConstraintModes
3818
          rValues => Var % ConstraintModes(iMode,:)
3819
          DO j=1,3          
3820
            No = No + 1
3821
            IF( ALL(DofIndexes(np+1:nd) > 0 ) ) THEN            
3822
              Values(No) = SUM( WBasis(1:nd-np,j) * rValues(DofIndexes(np+1:nd)))
3823
            ELSE
3824
              Values(No) = 0.0_dp
3825
            END IF
3826
          END DO
3827
        END DO
3828
      ELSE
3829
        DO j=1,3
3830
          No = No + 1
3831
          IF( ALL(DofIndexes(np+1:nd) > 0 ) ) THEN
3832
            Values(No) = SUM( WBasis(1:nd-np,j) * Var % Values(DofIndexes(np+1:nd)))
3833
          ELSE
3834
            Values(No) = 0.0_dp
3835
          END IF
3836
        END DO
3837
      END IF
3838
      IF( AVBasis ) THEN
3839
        No = No + 1
3840
        Values(No) = 0.0_dp
3841
        IF( ALL(DofIndexes(1:np) > 0 ) ) THEN
3842
          Values(No) = SUM( Basis(1:np) * Var % Values(DofIndexes(1:np)))
3843
        END IF
3844
      END IF
3845

3846
    ELSE IF ( IpVar ) THEN
3847
      l = 0 
3848
      IF( SIZE( Var % Perm ) > Element2 % ElementIndex ) THEN 
3849
        i1 = Var % Perm(Element2 % ElementIndex)
3850
        i2 = Var % Perm(Element2 % ElementIndex+1)
3851
        l = i2-i1
3852
      END IF
3853

3854
      IF(l>0) THEN
3855
        IF( .NOT. ALLOCATED(fip)) THEN
3856
          ALLOCATE( fip(l) )
3857
        ELSE IF (SIZE(fip) < l) THEN
3858
          DEALLOCATE( fip )
3859
          ALLOCATE( fip(l) )
3860
        END IF
3861

3862
        IF(.NOT. ALLOCATED(fdg)) THEN
3863
           ALLOCATE( fdg(n) )
3864
        ELSE IF(SIZE(fdg) < n) THEN
3865
           DEALLOCATE( fdg )
3866
           ALLOCATE( fdg(n) )
3867
        END IF
3868

3869
        DO ii=1,MAX(Var % Dofs,comps)
3870
          IF( Var % Dofs > 1 ) THEN
3871
            CONTINUE
3872
          ELSE          
3873
            IF( ii == 1 ) THEN
3874
              pVar => Var
3875
            ELSE IF( ii == 2 ) THEN
3876
              pVar => Var2
3877
            ELSE IF( ii == 3 ) THEN
3878
              pVar => Var3
3879
            END IF
3880
            fip(1:l) = pVar % Values(i1+1:i2)
3881
          END IF
3882

3883
          CALL Ip2DgFieldInElement( Mesh, Element2, l, fip, n2, fdg )              
3884
          Values(No+ii) = SUM( PtoBasis(1:n) * fdg(1:n) )
3885
        END DO
3886
      END IF
3887
      
3888
      No = No + MAX(Var % Dofs, Comps )      
3889
    ELSE
3890
      IF(.NOT. ASSOCIATED(pToBasis)) THEN
3891
        CALL Fatal('EvaluteVariableAtGivenPoint',&
3892
            'pToBasis not associated for variable: '//TRIM(Var % Name))
3893
      END IF
3894
      IF(.NOT. ASSOCIATED(pToIndexes)) THEN
3895
        CALL Fatal('EvaluteVariableAtGivenPoint',&
3896
            'pToIndexes not associated for variable: '//TRIM(Var % Name))
3897
      END IF
3898

3899
      ! This maybe should be checked earlier, but better late than never perhaps ?
3900
      IF (DGVar .AND. pElem ) nd = Element2 % Type % NumberOfNodes
3901
      
3902
      IF( ASSOCIATED(Var % Perm) ) THEN
3903
        DofIndexes(1:nd) = Var % Perm(PToIndexes(1:nd))
3904
      ELSE
3905
        DofIndexes(1:nd) = PtoIndexes(1:nd)
3906
      END IF
3907
      
3908
        IF( IsGrad ) THEN          
3909
          IF( Var % Dofs /= 1 ) THEN
3910
            CALL Fatal('EvaluteVariableAtGivenPoint','Gradient only possible for one dof!')
3911
          END IF
3912
          IF( ALL(Dofindexes(1:nd) > 0 ) ) THEN
3913
            DO ii=1,3
3914
              Values(No+ii) = SUM( dBasisdx(1:nd,ii) * Var % Values(DofIndexes(1:nd)))
3915
            END DO
3916
          ELSE
3917
            Values(No+1:No+3) = 0._dp
3918
          END IF
3919
          No = No + 3
3920
        ELSE IF( IsEigen ) THEN          
3921
          DO iMode = 1, NoEigenValues
3922
            cValues => Var % EIgenVectors(iMode,:)            
3923
            IF( ALL(Dofindexes(1:nd) > 0 ) ) THEN
3924
              DO ii=1,Var % DOfs              
3925
                Values(No+ii) = Values(No+ii) + SUM( PtoBasis(1:nd) * &
3926
                    cValues(Var%Dofs*(DofIndexes(1:nd)-1)+ii))
3927
              END DO
3928
            ELSE
3929
              Values(No+1:No+Var % Dofs)=0._dp      
3930
            END IF
3931
            No = No + Var % Dofs
3932
          END DO
3933
        ELSE IF( IsModes ) THEN          
3934
          DO iMode = 1, NoConstraintModes
3935
            rValues => Var % ConstraintModes(iMode,:)            
3936
            IF( ALL(Dofindexes(1:nd) > 0 ) ) THEN
3937
              DO ii=1,Var % DOfs              
3938
                Values(No+ii) = Values(No+ii) + SUM( PtoBasis(1:nd) * &
3939
                    rValues(Var%Dofs*(DofIndexes(1:nd)-1)+ii))
3940
              END DO
3941
            ELSE
3942
              Values(No+1:No+Var % Dofs)=0._dp      
3943
            END IF
3944
            No = No + Var % Dofs
3945
          END DO
3946
        ELSE
3947
          IF( ALL(Dofindexes(1:nd) > 0 ) ) THEN
3948
            IF( Var % Dofs > 1 ) THEN
3949
              DO ii=1,Var % Dofs
3950
                Values(No+ii) = SUM( PtoBasis(1:nd) * &
3951
                    Var % Values(Var%Dofs*(DofIndexes(1:nd)-1)+ii))
3952
              END DO
3953
            ELSE
3954
              Values(No+1) = SUM(PToBasis(1:nd) * Var % Values(DofIndexes(1:nd)))            
3955
              IF( comps >= 2 ) Values(No+2) = SUM(PToBasis(1:nd) * Var2 % Values(DofIndexes(1:nd)))
3956
              IF( comps >= 3 ) Values(No+3) = SUM(PToBasis(1:nd) * Var3 % Values(DofIndexes(1:nd)))            
3957
            END IF
3958
          ELSE
3959
            Values(No+1:No+MAX( Var % Dofs, comps ))=0._dp
3960
          ENDIF
3961
          No = No + MAX( Var % Dofs, comps )
3962
        END IF
3963
      END IF
3964

3965
#if 0
3966
    ! Debugging code
3967
    IF( .NOT. EdgeBasis .AND. ASSOCIATED(Var % Perm) .AND. No > 0 .AND. ASSOCIATED( PToIndexes ) ) THEN
3968
      BLOCK
3969
        REAL(KIND=dp) :: vmin, vmax
3970
        INTEGER :: kmax
3971
        PRINT *,'VarInfo:',TRIM(Var % Name), nd, n, pToIndexes(1:nd)
3972

3973
        kmax = 1
3974
        DO i=1,nd
3975
          IF( PtoBasis(i) > pToBasis(kmax)) kmax = i
3976
        END DO
3977

3978
        IF( ASSOCIATED(pToIndexes) ) THEN
3979
          DofIndexes(1:nd) = Var % Perm(PToIndexes(1:nd))
3980
        ELSE
3981
          DofIndexes(1:nd) = PtoIndexes(1:nd)
3982
        END IF
3983
          
3984
        vmin = MINVAL(Var % Values(DofIndexes(1:nd)),DofIndexes(1:nd)>0)
3985
        vmax = MAXVAL(Var % Values(DofIndexes(1:nd)),DofIndexes(1:nd)>0)
3986

3987
        PRINT *,'PtoIndexes:',PtoIndexes(1:nd)
3988
        PRINT *,'PtoBasis',SUM(PToBasis(1:nd)),PToBasis(1:nd)
3989
        PRINT *,'ElemRange:',vmin< Values(1) .AND. vmax > Values(1), &
3990
            vmin, vmax, Var % Values(Var % Perm(PToIndexes(kmax))), Values(1:No)
3991
      END BLOCK
3992
    END IF
3993
#endif
3994
      
3995
  CONTAINS
3996

3997
    FUNCTION PickDgIndexes(Element,PParent) RESULT ( PToInds) 
3998
      TYPE(Element_t), POINTER :: Element
3999
      INTEGER, POINTER :: PtoInds(:)
4000
      TYPE(Element_t), POINTER, OPTIONAL :: PParent
4001

4002
      TYPE(Element_t), POINTER :: Parent
4003
      INTEGER :: i,j,lr
4004
            
4005
      IF( ASSOCIATED( Element % DgIndexes ) ) THEN
4006
        PToInds => Element % DGIndexes
4007
      ELSE IF( ASSOCIATED( Element % BoundaryInfo ) ) THEN
4008
        DO lr=1,2
4009
          IF(lr==1) THEN
4010
            Parent => Element % BoundaryInfo % Left
4011
          ELSE
4012
            Parent => Element % BoundaryInfo % Right
4013
          END IF
4014
          IF(.NOT. ASSOCIATED( Parent ) ) CYCLE
4015
          IF(.NOT. ASSOCIATED( Parent % DGIndexes ) ) CYCLE
4016
          IF( PRESENT(PParent)) THEN
4017
            IF(.NOT. ASSOCIATED(Parent, PParent) ) CYCLE
4018
          END IF          
4019
          IF( ALL( Var % Perm( Parent % DGIndexes ) /= 0) ) THEN                  
4020
            DO i=1,Element % TYPE % NumberOfNodes
4021
              DO j=1,Parent % TYPE % NumberOfNodes
4022
                IF( Element % NodeIndexes(i) == Parent % NodeIndexes(j) ) THEN
4023
                  DGIndexes(i) = Parent % DGIndexes(j)
4024
                  EXIT
4025
                END IF
4026
              END DO
4027
            END DO
4028
            EXIT
4029
          END IF
4030
        END DO
4031
        PtoInds => DGIndexes
4032
      END IF
4033
    END FUNCTION PickDgIndexes
4034
    
4035
  END SUBROUTINE EvaluateVariableAtGivenPoint
4036

4037

4038

4039

4040
!------------------------------------------------------------------------------
4041
   SUBROUTINE mGetBoundaryIndexesFromParent( Mesh, Element, Indexes, indSize )
4042
!------------------------------------------------------------------------------
4043
     IMPLICIT NONE
4044

4045
     ! Parameters
4046
     TYPE(Mesh_t) :: Mesh
4047
     TYPE(Element_t), POINTER :: Element
4048
     INTEGER :: indSize, Indexes(:)
4049
     
4050
     ! Variables
4051
     TYPE(Element_t), POINTER :: Edge, Face
4052
     INTEGER :: i,j,n
4053
     TYPE(Element_t), POINTER :: Parent
4054
     
4055
     ! Clear indexes
4056
     Indexes = 0
4057
     indSize = 0
4058

4059
     Parent => Element % pDefs % localParent
4060
     IF ( .NOT. ASSOCIATED(Parent) ) RETURN
4061
             
4062
     n = Element % TYPE % NumberOfNodes
4063

4064
     ! Nodal indexes
4065
     Indexes(1:n) = Element % NodeIndexes(1:n)
4066
     indSize = n
4067
     
4068
     
4069
     ! Assign rest of indexes if necessary
4070
     SELECT CASE(Parent % TYPE % DIMENSION)
4071
     CASE (1)
4072
       CONTINUE
4073

4074
     CASE (2)
4075
       IF(.NOT. ASSOCIATED(Mesh % Edges) ) RETURN
4076
         
4077
       ! Add index for each bubble dof in edge
4078
        DO i=1,Element % BDOFs
4079
           n = n+1
4080
           
4081
           IF (SIZE(Indexes) < n) THEN
4082
              CALL Warn('mGetBoundaryIndexes','Not enough space reserved for indexes')
4083
              RETURN
4084
           END IF
4085

4086
           Indexes(n) = Mesh % NumberOfNodes + &
4087
                (Parent % EdgeIndexes(Element % PDefs % localNumber)-1) * Mesh % MaxEdgeDOFs + i
4088
        END DO
4089
     
4090
        indSize = n 
4091
      CASE (3)
4092
        IF(.NOT. ( ASSOCIATED(Mesh % Faces) .AND. ASSOCIATED(Mesh % Edges) ) ) RETURN
4093
        IF(.NOT. ASSOCIATED(Element % PDefs) ) RETURN        
4094
        IF(Element % PDefs % LocalNumber == 0 ) RETURN
4095
                
4096
        ! Get boundary face
4097
        Face => Mesh % Faces( Parent % FaceIndexes(Element % PDefs % localNumber) )
4098
        
4099
        ! Add indexes of faces edges 
4100
        DO i=1, Face % TYPE % NumberOfEdges
4101
           Edge => Mesh % Edges( Face % EdgeIndexes(i) )
4102
           
4103
           ! If edge has no dofs jump to next edge
4104
           IF (Edge % BDOFs <= 0) CYCLE
4105

4106
           DO j=1,Edge % BDOFs
4107
              n = n + 1
4108
              
4109
              IF (SIZE(Indexes) < n) THEN
4110
                 CALL Warn('mGetBoundaryIndexes','Not enough space reserved for indexes')
4111
                 RETURN
4112
              END IF
4113
              
4114
              Indexes(n) = Mesh % NumberOfNodes +&
4115
                  ( Face % EdgeIndexes(i)-1)*Mesh % MaxEdgeDOFs + j
4116
           END DO
4117
        END DO
4118
               
4119
        ! Add indexes of faces bubbles
4120
        DO i=1,Face % BDOFs
4121
           n = n + 1
4122

4123
           IF (SIZE(Indexes) < n) THEN
4124
              CALL Warn('mGetBoundaryIndexes','Not enough space reserved for indexes')
4125
              RETURN
4126
           END IF
4127

4128
           Indexes(n) = Mesh % NumberOfNodes + &
4129
                Mesh % NumberOfEdges * Mesh % MaxEdgeDOFs + &
4130
                (Parent % FaceIndexes( Element % PDefs % localNumber )-1) * Mesh % MaxFaceDOFs + i
4131
        END DO        
4132

4133
        indSize = n
4134
     CASE DEFAULT
4135
        CALL Fatal('mGetBoundaryIndexes','Unsupported dimension')
4136
     END SELECT
4137
!------------------------------------------------------------------------------
4138
   END SUBROUTINE mGetBoundaryIndexesFromParent
4139
!------------------------------------------------------------------------------
4140

4141
!------------------------------------------------------------------------------
4142
  FUNCTION Find_Edge(Mesh,Parent,Element) RESULT(ptr)
4143
!------------------------------------------------------------------------------
4144
    TYPE(Element_t), POINTER :: Ptr
4145
    TYPE(Mesh_t) :: Mesh
4146
    TYPE(Element_t) :: Parent, Element
4147

4148
    INTEGER :: i,j,k,n
4149

4150
    Ptr => NULL()
4151
    DO i=1,Parent % TYPE % NumberOfEdges
4152
      Ptr => Mesh % Edges(Parent % EdgeIndexes(i))
4153
      n=0
4154
      DO j=1,Ptr % TYPE % NumberOfNodes
4155
        DO k=1,Element % TYPE % NumberOfNodes
4156
          IF (Ptr % NodeIndexes(j) == Element % NodeIndexes(k)) n=n+1
4157
        END DO
4158
      END DO
4159
      IF (n==Ptr % TYPE % NumberOfNodes) EXIT
4160
    END DO
4161
!------------------------------------------------------------------------------
4162
  END FUNCTION Find_Edge
4163
!------------------------------------------------------------------------------
4164

4165
!------------------------------------------------------------------------------
4166
  FUNCTION Find_Face(Mesh,Parent,Element) RESULT(ptr)
4167
!------------------------------------------------------------------------------
4168
    TYPE(Element_t), POINTER :: Ptr
4169
    TYPE(Mesh_t) :: Mesh
4170
    TYPE(Element_t) :: Parent, Element
4171

4172
    INTEGER :: i,j,k,n
4173

4174
    Ptr => NULL()
4175
    DO i=1,Parent % TYPE % NumberOfFaces
4176
      Ptr => Mesh % Faces(Parent % FaceIndexes(i))
4177
      n=0
4178
      DO j=1,Ptr % TYPE % NumberOfNodes
4179
        DO k=1,Element % TYPE % NumberOfNodes
4180
          IF (Ptr % NodeIndexes(j) == Element % NodeIndexes(k)) n=n+1
4181
        END DO
4182
      END DO
4183
      IF (n==Ptr % TYPE % NumberOfNodes) EXIT
4184
    END DO
4185
!------------------------------------------------------------------------------
4186
  END FUNCTION Find_Face
4187
!------------------------------------------------------------------------------
4188

4189

4190
   
4191
END MODULE ElementUtils
4192

4193
!> \} ElmerLib
4194

4195
Product

Resources

Company
