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
#ifdef HAVE_TRILINOS
91
     INTERFACE
92
      !! destroy the data structures (should be called when the matrix has
93
      !! to be updated and SolveTrilinos1 has to be called again).
94
      SUBROUTINE SolveTrilinos4(triliContainer) BIND(C,name='SolveTrilinos4')
95
        USE, INTRINSIC :: iso_c_binding
96
        INTEGER(KIND=C_INTPTR_T) :: triliContainer
97
      END SUBROUTINE SolveTrilinos4
98

99
     END INTERFACE
100
#endif
101

102
!------------------------------------------------------------------------------
103

104
     IF ( .NOT. ASSOCIATED( Matrix ) ) RETURN
105

106
     CALL DirectSolver( Matrix,x,b,Solver,Free_Fact=.TRUE.)
107

108
     IF ( ASSOCIATED( Matrix % Perm ) )        DEALLOCATE( Matrix % Perm )
109
     IF ( ASSOCIATED( Matrix % InvPerm ) )     DEALLOCATE( Matrix % InvPerm )
110

111
     IF ( ASSOCIATED( Matrix % Rows ) ) THEN
112
        IF ( ASSOCIATED(Matrix % Rows, Matrix % ILURows) .OR. SIZE(Matrix % Rows)==0 )  &
113
           Matrix % ILURows => Null()
114
        DEALLOCATE( Matrix % Rows )
115
     END IF
116

117
     IF ( ASSOCIATED( Matrix % Cols ) ) THEN
118
        IF ( ASSOCIATED(Matrix % Cols, Matrix % ILUCols) .OR. SIZE(Matrix % Cols)==0 ) &
119
           Matrix % ILUCols => Null()
120
        DEALLOCATE( Matrix % Cols )
121
     END IF
122

123
     IF ( ASSOCIATED( Matrix % Diag ) ) THEN
124
        IF ( ASSOCIATED(Matrix % Diag, Matrix % ILUDiag) .OR. SIZE(Matrix % Diag)==0 ) &
125
          Matrix % ILUDiag => Null()
126
        DEALLOCATE( Matrix % Diag )
127
     END IF
128

129
     IF ( ASSOCIATED( Matrix % RHS   ) )       DEALLOCATE( Matrix % RHS )
130
     IF ( ASSOCIATED( Matrix % Force ) )       DEALLOCATE( Matrix % Force )
131
     IF ( ASSOCIATED( Matrix % RHS_im ) )      DEALLOCATE( Matrix % RHS_im )
132

133
     IF ( ALLOCATED( Matrix % ExtraVals ) )    DEALLOCATE( Matrix % ExtraVals )
134
     IF ( ASSOCIATED( Matrix % Values ) )      DEALLOCATE( Matrix % Values )
135
     IF ( ASSOCIATED( Matrix % Values_im ) )   DEALLOCATE( Matrix % Values_im )
136
     IF ( ASSOCIATED( Matrix % MassValues ) )  DEALLOCATE( Matrix % MassValues )
137
     IF ( ASSOCIATED( Matrix % DampValues ) )  DEALLOCATE( Matrix % DampValues )
138
     IF ( ASSOCIATED( Matrix % PrecValues ) )  DEALLOCATE( Matrix % PrecValues )
139
     IF ( ASSOCIATED( Matrix % BulkValues ) )  DEALLOCATE( Matrix % BulkValues )
140
     IF ( ASSOCIATED( Matrix % BulkRHS   ) )   DEALLOCATE( Matrix % BulkRHS )
141

142
     IF ( ASSOCIATED( Matrix % ILUValues ) )   DEALLOCATE( Matrix % ILUValues )
143
     IF ( ASSOCIATED( Matrix % ILURows ) )     DEALLOCATE( Matrix % ILURows )
144
     IF ( ASSOCIATED( Matrix % ILUCols ) )     DEALLOCATE( Matrix % ILUCols )
145
     IF ( ASSOCIATED( Matrix % ILUDiag ) )     DEALLOCATE( Matrix % ILUDiag )
146

147
     IF ( ASSOCIATED( Matrix % CRHS   ) )      DEALLOCATE( Matrix % CRHS )
148
     IF ( ASSOCIATED( Matrix % CForce ) )      DEALLOCATE( Matrix % CForce )
149

150
     IF ( ASSOCIATED( Matrix % CValues ) )     DEALLOCATE( Matrix % CValues )
151
     IF ( ASSOCIATED( Matrix % CILUValues ) )  DEALLOCATE( Matrix % CILUValues )
152

153
     IF ( ASSOCIATED(Matrix % CMassValues) )  DEALLOCATE( Matrix % CMassValues )
154
     IF ( ASSOCIATED(Matrix % CDampValues) )  DEALLOCATE( Matrix % CDampValues )
155

156
     IF ( ALLOCATED( Matrix % GRows ) )      DEALLOCATE( Matrix % GRows )
157
     IF ( ALLOCATED( Matrix % RowOwner ) )   DEALLOCATE( Matrix % RowOwner )
158
     IF ( ASSOCIATED( Matrix % GOrder) )      DEALLOCATE( Matrix % GOrder )
159

160
     CALL FreeMatrix( Matrix % CircuitMatrix )
161
     CALL FreeMatrix( Matrix % AddMatrix )
162
     CALL FreeMatrix( Matrix % EMatrix )
163
     CALL FreeMatrix( Matrix % ConstraintMatrix )
164
     CALL FreeMatrix( Matrix % CollectionMatrix )
165

166
     IF  (ALLOCATED(Matrix % Dvalues) ) DEALLOCATE(Matrix % DValues)
167
     IF  (ALLOCATED(Matrix % ConstrainedDOF) ) DEALLOCATE(Matrix % ConstrainedDOF)
168
     IF  (ASSOCIATED(Matrix % EPerm) ) DEALLOCATE(Matrix % EPerm)
169
     IF  (ASSOCIATED(Matrix % BulkResidual) ) DEALLOCATE(Matrix % BulKResidual)
170
     IF  (ASSOCIATED(Matrix % DiagScaling) ) DEALLOCATE(Matrix % DiagScaling)
171
     IF  (ASSOCIATED(Matrix % Tvalues) ) DEALLOCATE(Matrix % Tvalues)
172
     IF  (ASSOCIATED(Matrix % BulkMassValues) ) DEALLOCATE(Matrix % BulkMassValues)
173
     IF  (ASSOCIATED(Matrix % BulkDampValues) ) DEALLOCATE(Matrix % BulkDampValues)
174
     IF  (ASSOCIATED(Matrix % HaloValues) ) DEALLOCATE(Matrix % HaloValues)
175
     IF  (ASSOCIATED(Matrix % HaloMassValues) ) DEALLOCATE(Matrix % HaloMassValues)
176

177
     IF(ASSOCIATED(Matrix % ParallelInfo)) THEN
178
       IF(ASSOCIATED(Matrix % ParallelInfo % GlobalDOFs)) DEALLOCATE(Matrix % ParallelInfo % GlobalDOFs)
179
       IF(ASSOCIATED(Matrix % ParallelInfo % GInterface)) DEALLOCATE(Matrix % ParallelInfo % GInterface)
180
       IF(ASSOCIATED(Matrix % ParallelInfo % Gorder)) DEALLOCATE(Matrix % ParallelInfo % GOrder)
181
 
182
       IF(ASSOCIATED(Matrix % ParallelInfo % NeighbourList)) THEN
183
         DO i=1,SIZE(Matrix % ParallelInfo % NeighbourList)
184
           IF (ASSOCIATED(Matrix % ParallelInfo % NeighbourList(i) % Neighbours)) &
185
             DEALLOCATE(Matrix % ParallelInfo % NeighbourList(i) % Neighbours)
186
         END DO
187
         DEALLOCATE(Matrix % ParallelInfo % NeighbourList)
188
       END IF
189

190
       IF(ASSOCIATED(Matrix % ParallelInfo % FaceNeighbourList)) THEN
191
         DO i=1,SIZE(Matrix % ParallelInfo % FaceNeighbourList)
192
           IF (ASSOCIATED(Matrix % ParallelInfo % FaceNeighbourList(i) % Neighbours)) &
193
             DEALLOCATE(Matrix % ParallelInfo % FaceNeighbourList(i) % Neighbours)
194
         END DO
195
         DEALLOCATE(Matrix % ParallelInfo % FaceNeighbourList)
196
       END IF
197
       IF(ASSOCIATED(Matrix % ParallelInfo % FaceInterface)) DEALLOCATE(Matrix % ParallelInfo % FaceInterface)
198

199
       IF(ASSOCIATED(Matrix % ParallelInfo % EdgeNeighbourList)) THEN
200
         DO i=1,SIZE(Matrix % ParallelInfo % EdgeNeighbourList)
201
           IF (ASSOCIATED(Matrix % ParallelInfo % EdgeNeighbourList(i) % Neighbours)) &
202
             DEALLOCATE(Matrix % ParallelInfo % EdgeNeighbourList(i) % Neighbours)
203
         END DO
204
         DEALLOCATE(Matrix % ParallelInfo % EdgeNeighbourList)
205
       END IF
206
       IF(ASSOCIATED(Matrix % ParallelInfo % EdgeInterface)) DEALLOCATE(Matrix % ParallelInfo % EdgeInterface)
207
       DEALLOCATE(Matrix % ParallelInfo)
208
     END IF
209
         
210

211
     p=>Matrix % ParMatrix
212
     IF(ASSOCIATED(p)) THEN
213
       s => Matrix % ParMatrix % SplittedMatrix
214

215
       s % InsideMatrix % EMatrix => NULL()
216
       CALL FreeMatrix(s % InsideMatrix)
217

218
       DO i=1,SIZE(s % IfMatrix)
219
         m => s % IfMatrix(i)
220
         IF(ASSOCIATED(m)) THEN
221
           IF(ALLOCATED(m % Rows)) DEALLOCATE(m % Rows)
222
           IF(ALLOCATED(m % Cols)) DEALLOCATE(m % Cols)
223
           IF(ALLOCATED(m % Diag)) DEALLOCATE(m % Diag)
224
           IF(ALLOCATED(m % Grows)) DEALLOCATE(m % Grows)
225
           IF(ALLOCATED(m % Values)) DEALLOCATE(m % Values)
226
           IF(ALLOCATED(m % RowOwner)) DEALLOCATE(m % RowOwner)
227
           IF(ALLOCATED(m % ILUValues)) DEALLOCATE(m % ILUValues)
228
           IF(ALLOCATED(m % MassValues)) DEALLOCATE(m % MassValues)
229
           IF(ALLOCATED(m % DampValues)) DEALLOCATE(m % DampValues)
230
           IF(ALLOCATED(m % PrecValues)) DEALLOCATE(m % PrecValues)
231
         END IF
232
       END DO
233
       DEALLOCATE(s % IfMatrix)
234

235
       DO i=1,SIZE(s % NbsIfMatrix)
236
         m => s % NbsIfMatrix(i)
237
         IF(ASSOCIATED(m)) THEN
238
           IF(ALLOCATED(m % Rows)) DEALLOCATE(m % Rows)
239
           IF(ALLOCATED(m % Cols)) DEALLOCATE(m % Cols)
240
           IF(ALLOCATED(m % Diag)) DEALLOCATE(m % Diag)
241
           IF(ALLOCATED(m % Grows)) DEALLOCATE(m % Grows)
242
           IF(ALLOCATED(m % Values)) DEALLOCATE(m % Values)
243
           IF(ALLOCATED(m % RowOwner)) DEALLOCATE(m % RowOwner)
244
           IF(ALLOCATED(m % ILUValues)) DEALLOCATE(m % ILUValues)
245
           IF(ALLOCATED(m % MassValues)) DEALLOCATE(m % MassValues)
246
           IF(ALLOCATED(m % DampValues)) DEALLOCATE(m % DampValues)
247
           IF(ALLOCATED(m % PrecValues)) DEALLOCATE(m % PrecValues)
248
         END IF
249
       END DO
250
       DEALLOCATE(s % NbsIfMatrix)
251

252
       IF(ASSOCIATED(s % VecIndices)) THEN
253
         DO i=1,SIZE(s % VecIndices)
254
           IF(ASSOCIATED(s % Vecindices(i) % RevInd)) DEALLOCATE(s % VecIndices(i) % RevInd)
255
         END DO
256
         DEALLOCATE(s % VecIndices)
257
       END IF
258

259
       IF(ASSOCIATED(s % IfVecs)) THEN
260
         DO i=1,SIZE(s % IfVecs)
261
           IF(ASSOCIATED(s % IfVecs(i) % IfVec)) DEALLOCATE(s % IfVecs(i) % IfVec)
262
         END DO
263
         DEALLOCATE(s % ifVecs)
264
       END IF
265

266
       IF(ASSOCIATED(s % IfORows)) THEN
267
         DO i=1,SIZE(s % IfORows)
268
           IF(ASSOCIATED(s % IfORows(i) % IfVec)) DEALLOCATE(s % IfORows(i) % IfVec)
269
         END DO
270
         DEALLOCATE(s % ifORows)
271
       END IF
272

273
       IF(ASSOCIATED(s % IfLCols)) THEN
274
         DO i=1,SIZE(s % IfLCols)
275
           IF(ASSOCIATED(s % IfLCols(i) % IfVec)) DEALLOCATE(s % IfLCols(i) % IfVec)
276
         END DO
277
         DEALLOCATE(s % ifLCols)
278
       END IF
279

280
       IF(ASSOCIATED(s % ResBuf)) THEN
281
         DO i=1,SIZE(s % ResBuf)
282
           IF(ALLOCATED(s % ResBuf(i) % ResVal)) DEALLOCATE(s % ResBuf(i) % ResVal)
283
           IF(ALLOCATED(s % ResBuf(i) % ResInd)) DEALLOCATE(s % ResBuf(i) % ResInd)
284
         END DO
285
         DEALLOCATE(s % ResBuf)
286
       END IF
287

288
       IF(ASSOCIATED(s % RHS)) THEN
289
         DO i=1,SIZE(s % RHS)
290
           IF(ASSOCIATED(s % RHS(i) % RHSVec)) DEALLOCATE(s % RHS(i) % RHSVec)
291
           IF(ASSOCIATED(s % RHS(i) % RHSInd)) DEALLOCATE(s % RHS(i) % RHSInd)
292
         END DO
293
         DEALLOCATE(s % RHS)
294
       END IF
295

296
       IF (ASSOCIATED(s % Work)) DEALLOCATE(s % Work)
297
       IF (ASSOCIATED(s % TmpXVec)) DEALLOCATE(s % TmpXVec)
298
       IF (ASSOCIATED(s % TmpRVec)) DEALLOCATE(s % TmpRVec)
299

300
       IF(ASSOCIATED(s % GlueTable)) THEN
301
         IF(ASSOCIATED(s % GlueTable % Rows)) DEALLOCATE(s % GlueTable % Rows)
302
         IF(ASSOCIATED(s % GlueTable % Cols)) DEALLOCATE(s % GlueTable % Cols)
303
         IF(ASSOCIATED(s % GlueTable % Inds)) DEALLOCATE(s % GlueTable % Inds)
304
         IF(ASSOCIATED(s % GlueTable % RowOwner)) DEALLOCATE(s % GlueTable % RowOwner)
305
         DEALLOCATE(s % GlueTable)
306
       END IF
307
       DEALLOCATE(s)
308

309
       IF(ASSOCIATED(p % ParEnv % Active)) THEN
310
         DEALLOCATE(p % ParEnv % Active)
311
       END IF
312

313
       IF(ASSOCIATED(p % ParEnv % Isneighbour)) THEN
314
         DEALLOCATE(p % ParEnv % Isneighbour)
315
       END IF
316

317
       DEALLOCATE(p)
318
     END IF
319

320
#ifdef HAVE_HYPRE
321
     IF (Matrix % Hypre /= 0) THEN
322
       i = ListGetInteger( CurrentModel % Simulation,'Max Output Level',Found ) 
323
       IF(ParEnv % MyPe /= 0) i = 0 
324
       CALL SolveHypre4(Matrix % Hypre, i )
325
       Matrix % Hypre = 0
326
     END IF
327
#endif
328

329
#ifdef HAVE_TRILINOS
330
     IF (Matrix % Trilinos /= 0) THEN
331
       CALL SolveTrilinos4(Matrix % Trilinos)
332
     END IF
333
#endif
334
     DEALLOCATE( Matrix )
335
!------------------------------------------------------------------------------
336
   END SUBROUTINE FreeMatrix
337
!------------------------------------------------------------------------------
338

339

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

380
    CALL Info(Caller,'Creating list matrix',Level=14)
381

382
    GB = .FALSE.
383
    IF ( PRESENT(GlobalBubbles) ) GB = GlobalBubbles
384

385
    IF( DgSolver ) THEN
386
      DB = ListGetLogical( Solver % Values,'DG Reduced Basis',Found) 
387
    ELSE
388
      DB = .FALSE.
389
    END IF
390

391
    ! When this has been checked properly the old can be removed
392
    PSA = ListGetLogical( Solver % Values,'PSA',Found ) 
393
    IF(.NOT. Found) PSA = .TRUE.
394
    
395
    List => List_AllocateMatrix(LocalNodes)
396

397
    NDOFs = Mesh % MaxNDOFs
398
    BDOFs = Mesh % MaxBDOFs
399
    EDOFs = Mesh % MaxEdgeDOFs
400
    FDOFs = Mesh % MaxFaceDOFs
401
    IF( PRESENT( NodalDofsOnly ) ) THEN
402
      IF( NodalDofsOnly ) THEN
403
        EDOFS = 0
404
        FDOFS = 0
405
      END IF
406
    END IF
407

408
    IF( PRESENT( ProjectorDofs ) ) THEN
409
      DoProjectors = ProjectorDofs 
410
    ELSE
411
      DoProjectors = .TRUE.
412
    END IF
413

414
    IF( EDOFS > 0 .AND. .NOT. ASSOCIATED(Mesh % Edges) ) THEN
415
      CALL Warn(Caller,'Edge dofs requested but not edges exist in mesh!')
416
      EDOFS = 0
417
    END IF
418

419
    IF( FDOFS > 0 .AND. .NOT. ASSOCIATED(Mesh % Faces) ) THEN
420
      CALL Warn(Caller,'Face dofs requested but not faces exist in mesh!')
421
      FDOFS = 0
422
    END IF
423

424
    IndexSize = 128
425
    ALLOCATE( Indexes(IndexSize), STAT=istat )
426
    IF( istat /= 0 ) THEN
427
      CALL Fatal(Caller,'Allocation error for Indexes')
428
    END IF
429

430
    ! Create sparse matrix for the Discontinuous Galerkin solver 
431
    ! Using either reduced or full basis
432
    !-------------------------------------------------------------------
433
    FoundDG = .FALSE.
434
    IF ( DGSolver ) THEN    
435
      ! Create the sparse matrix for the Discontinuous Bodies solver
436
      !-------------------------------------------------------------------
437
      IF ( DB ) THEN
438
        DGIndirect = ListGetLogical( Solver % Values,'DG Indirect Connections',Found )
439

440
        IF( DGIndirect ) THEN
441
          CALL Info(Caller,'Creating also indirect connections!',Level=12)
442
          ALLOCATE( IndirectPairs( LocalNodes ), STAT=istat )
443
          IF ( istat /= 0 ) CALL Fatal(Caller,'Memory allocation error for IndirectPairs work.')
444
          IndirectPairs = 0
445
        END IF
446
                        
447

448
        DO t=1,Mesh % NumberOfBulkElements
449
          n = 0
450
          Elm => Mesh % Elements(t)
451
          IF ( .NOT. CheckElementEquation(Model,Elm,Equation) ) CYCLE
452

453
          FoundDG = FoundDG .OR. Elm % DGDOFs > 0
454
          DO j=1,Elm % DGDOFs
455
            n = n + 1
456
            Indexes(n) = Elm % DGIndexes(j)
457
          END DO
458

459
          IF( Elm % DGDofs /= Elm % TYPE % NumberOfNodes ) THEN
460
            CALL Fatal(Caller,'Mismatch in sizes in reduced basis DG!')
461
          END IF
462

463
          IF( PSA ) THEN
464
            CALL PackSortAdd(n,Indexes,Reorder)
465
          ELSE
466
            DO i=1,n
467
              k1 = Reorder(Indexes(i))
468
              IF ( k1 <= 0 ) CYCLE
469
              DO j=1,n
470
                k2 = Reorder(Indexes(j))
471
                IF ( k2 <= 0 ) CYCLE
472
                Lptr => List_GetMatrixIndex( List,k1,k2 )
473
              END DO
474
            END DO
475
          END IF
476
        END DO
477

478
        IF( Mesh % NumberOfFaces == 0 ) THEN
479
          DO t=1,Mesh % NumberOfEdges
480
            n = 0
481
            Elm => Mesh % Edges(t)
482
            IF(.NOT. ASSOCIATED( Elm % BoundaryInfo ) ) CYCLE
483

484
            Left => Elm % BoundaryInfo % Left
485
            IF(.NOT. ASSOCIATED( Left ) ) CYCLE
486
            IF ( .NOT. CheckElementEquation(Model,Left,Equation) ) CYCLE
487

488
            Right => Elm % BoundaryInfo %  Right
489
            IF(.NOT. ASSOCIATED( Right ) ) CYCLE
490
            IF ( .NOT. CheckElementEquation(Model,Right,Equation) ) CYCLE
491

492
            IF( Left % BodyId == Right % BodyId ) CYCLE
493

494
            IF( DGIndirect ) THEN
495
              DO i=1,Left % DGDOFs
496
                k1 = ReOrder( Left % DgIndexes(i) )
497
                DO j=1,Right % DGDOFs
498
                  IF( Left % NodeIndexes(i) == Right % NodeIndexes(j) ) THEN
499
                    k2 = ReOrder( Right % DgIndexes(j) )
500
                    IF( k1 /= k2 ) THEN
501
                      IndirectPairs( k1 ) = k2
502
                      EXIT
503
                    END IF
504
                  END IF
505
                END DO
506
              END DO
507
            END IF
508

509
            FoundDG = FoundDG .OR. Left % DGDOFs > 0
510
            DO j=1,Left % DGDOFs
511
              n = n + 1
512
              Indexes(n) = Left % DGIndexes(j)
513
            END DO
514

515
            FoundDG = FoundDG .OR. Right % DGDOFs > 0
516
            DO j=1,Right % DGDOFs
517
              n = n + 1
518
              Indexes(n) = Right % DGIndexes(j)
519
            END DO
520

521
            IF( PSA ) THEN
522
              CALL PackSortAdd(n,Indexes,Reorder)
523
            ELSE
524
              DO i=1,n
525
                k1 = Reorder(Indexes(i))
526
                IF ( k1 <= 0 ) CYCLE
527
                DO j=1,n
528
                  k2 = Reorder(Indexes(j))
529
                  IF ( k2 <= 0 ) CYCLE
530
                  Lptr => List_GetMatrixIndex( List,k1,k2 )
531
                END DO
532
              END DO
533
            END IF
534
          END DO
535
        END IF
536

537
        
538
        DO t=1,Mesh % NumberOfFaces
539
          n = 0
540

541
          Elm => Mesh % Faces(t)
542
          IF(.NOT. ASSOCIATED( Elm % BoundaryInfo ) ) CYCLE
543
          
544
          Left => Elm  % BoundaryInfo % Left
545
          IF(.NOT. ASSOCIATED( Left ) ) CYCLE
546
          IF ( .NOT. CheckElementEquation(Model,Left,Equation) ) CYCLE
547

548
          Right => Elm % BoundaryInfo %  Right
549
          IF(.NOT. ASSOCIATED( Right ) ) CYCLE
550
          IF ( .NOT. CheckElementEquation(Model,Right,Equation) ) CYCLE
551

552
          IF( Left % BodyId == Right % BodyId ) CYCLE
553

554
          IF( DGIndirect ) THEN
555
            DO i=1,Left % DGDOFs
556
              k1 = ReOrder( Left % DgIndexes(i) )
557
              DO j=1,Right % DGDOFs
558
                IF( Left % NodeIndexes(i) == Right % NodeIndexes(j) ) THEN
559
                  k2 = ReOrder( Right % DgIndexes(j) )
560
                  IF( k1 /= k2 ) THEN
561
                    IF( IndirectPairs(k1) > 0 .AND. IndirectPairs(k1) /= k2 ) THEN
562
                      PRINT *,'Problematic node:',k1,IndirectPairs(k1),k2
563
                    END IF
564
                    IndirectPairs( k1 ) = k2
565
                    EXIT
566
                  END IF
567
                END IF
568
              END DO
569
            END DO
570
          END IF            
571

572
          FoundDG = FoundDG .OR. Left % DGDOFs > 0
573
          DO j=1,Left % DGDOFs
574
            n = n + 1
575
            Indexes(n) = Left % DGIndexes(j)
576
          END DO
577

578
          FoundDG = FoundDG .OR. Right % DGDOFs > 0
579
          DO j=1,Right % DGDOFs
580
            n = n + 1
581
            Indexes(n) = Right % DGIndexes(j)
582
          END DO
583

584
          IF( PSA ) THEN
585
            CALL PackSortAdd(n,Indexes,Reorder)
586
          ELSE
587
            DO i=1,n
588
              k1 = Reorder(Indexes(i))
589
              IF ( k1 <= 0 ) CYCLE
590
              DO j=1,n
591
                k2 = Reorder(Indexes(j))
592
                IF ( k2 <= 0 ) CYCLE
593
                Lptr => List_GetMatrixIndex( List,k1,k2 )
594
              END DO
595
            END DO
596
          END IF
597
        END DO
598

599
        IF( DGIndirect ) THEN
600
          DO k1 = 1, LocalNodes
601
            k2 = IndirectPairs(k1)
602
            IF( k2 == 0 ) CYCLE 
603
            !PRINT *,'Exchange structure between rows:',k1,k2
604
            CALL List_ExchangeRowStructure( List, k1, k2 ) 
605
          END DO
606
        END IF
607
          
608
        
609
      ELSE
610
        ! Classical DG solver
611
        !-------------------------------------        
612
        DO t=1,Mesh % NumberOfEdges
613
          n = 0
614
          Elm => Mesh % Edges(t) % BoundaryInfo % Left
615
          IF ( ASSOCIATED( Elm ) ) THEN
616
            IF ( CheckElementEquation(Model,Elm,Equation) ) THEN
617
              FoundDG = FoundDG .OR. Elm % DGDOFs > 0
618
              DO j=1,Elm % DGDOFs
619
                n = n + 1
620
                Indexes(n) = Elm % DGIndexes(j)
621
              END DO
622
            END IF
623
          END IF
624

625
          Elm => Mesh % Edges(t) % BoundaryInfo %  Right
626
          IF ( ASSOCIATED( Elm ) ) THEN
627
            IF ( CheckElementEquation(Model,Elm,Equation) ) THEN
628
              FoundDG = FoundDG .OR. Elm % DGDOFs > 0
629
              DO j=1,Elm % DGDOFs
630
                n = n + 1
631
                Indexes(n) = Elm % DGIndexes(j)
632
              END DO
633
            END IF
634
          END IF
635

636
          IF( PSA ) THEN
637
            CALL PackSortAdd(n,Indexes,Reorder)
638
          ELSE
639
            DO i=1,n
640
              k1 = Reorder(Indexes(i))
641
              IF ( k1 <= 0 ) CYCLE
642
              DO j=1,n
643
                k2 = Reorder(Indexes(j))
644
                IF ( k2 <= 0 ) CYCLE
645
                Lptr => List_GetMatrixIndex( List,k1,k2 )
646
              END DO
647
            END DO
648
          END IF          
649
        END DO
650
        
651
        DO t=1,Mesh % NumberOfFaces
652
          n = 0
653
          Elm => Mesh % Faces(t) % BoundaryInfo % Left
654
          IF ( ASSOCIATED( Elm ) ) THEN
655
            IF ( CheckElementEquation(Model,Elm,Equation) ) THEN
656
              FoundDG = FoundDG .OR. Elm % DGDOFs > 0
657
              DO j=1,Elm % DGDOFs
658
                n = n + 1
659
                Indexes(n) = Elm % DGIndexes(j)
660
              END DO
661
            END IF
662
          END IF
663

664
          Elm => Mesh % Faces(t) % BoundaryInfo %  Right
665
          IF ( ASSOCIATED( Elm ) ) THEN
666
            IF ( CheckElementEquation(Model,Elm,Equation) ) THEN
667
              FoundDG = FoundDG .OR. Elm % DGDOFs > 0
668
              DO j=1,Elm % DGDOFs
669
                n = n + 1
670
                Indexes(n) = Elm % DGIndexes(j)
671
              END DO
672
            END IF
673
          END IF
674

675
          IF( PSA ) THEN
676
            CALL PackSortAdd(n,Indexes,Reorder)
677
          ELSE
678
            DO i=1,n
679
              k1 = Reorder(Indexes(i))
680
              IF ( k1 <= 0 ) CYCLE
681
              DO j=1,n
682
                k2 = Reorder(Indexes(j))
683
                IF ( k2 <= 0 ) CYCLE
684
                Lptr => List_GetMatrixIndex( List,k1,k2 )
685
              END DO
686
            END DO
687
          END IF            
688
        END DO
689
      END IF
690
    END IF ! DGSolver
691

692
    
693
    !     Diffuse gray radiation condition:
694
    !     ---------------------------------
695
    Radiation = ListGetLogical( Solver % Values, 'Radiation Solver', Found )
696
    IF( Radiation ) THEN        
697
      Radiation = .FALSE.
698
      DO i=1,Model % NumberOfBCs
699
        RadiationFlag = ListGetString( Model % BCs(i) % Values, 'Radiation', GotIt )
700
        IF (GotIt) THEN
701
          IF ( RadiationFlag == 'diffuse gray' ) THEN
702
            Radiation = .TRUE.
703
            EXIT
704
          END IF
705
        END IF
706
      END DO
707
    END IF
708
    
709
    IF ( Radiation ) THEN
710
      BLOCK
711
        INTEGER, ALLOCATABLE :: Inds(:),ElemInds(:),ElemInds2(:)
712
        INTEGER :: cnt, maxnodes
713

714
        n = Mesh % MaxElementNodes
715
        ALLOCATE( ElemInds(n), ElemInds2(n), STAT=istat )
716
        IF ( istat /= 0 ) CALL Fatal(Caller,'Memory allocation error for ElemInds.')
717
        
718
        ElemInds = 0
719
        ElemInds2 = 0
720

721
        ! Check the maximum number of nodes in boundary element
722
        maxnodes = 0
723
        DO i = Mesh % NumberOfBulkElements+1, &
724
            Mesh % NumberOfBulkElements + Mesh % NumberOfBoundaryElements
725
          
726
          Element => Mesh % Elements(i)
727

728
          IF( .NOT. ASSOCIATED( Element % BoundaryInfo ) ) CYCLE
729
          IF ( .NOT. ASSOCIATED(Element % BoundaryInfo % RadiationFactors) ) CYCLE
730
          
731
          n = Element % Type % NumberOfNodes
732
          maxnodes = MAX( n, maxnodes )
733
        END DO
734
        
735
        
736
        CALL Info(Caller,'Adding radiation matrix',Level=14)      
737
        DO i = Mesh % NumberOfBulkElements+1, &
738
            Mesh % NumberOfBulkElements + Mesh % NumberOfBoundaryElements
739
          
740
          Element => Mesh % Elements(i)
741

742
          IF( .NOT. ASSOCIATED( Element % BoundaryInfo ) ) CYCLE
743
          IF ( .NOT. ASSOCIATED(Element % BoundaryInfo % RadiationFactors) ) CYCLE
744

745
          NumberOfFactors = Element % BoundaryInfo % &
746
              RadiationFactors % NumberOfImplicitFactors
747
          IF(NumberOfFactors == 0 ) CYCLE
748
          
749
          n = Element % Type % NumberOfNodes
750
          
751
          IF( DgSolver ) THEN
752
            CALL DgRadiationIndexes(Element,n,ElemInds,.TRUE.)
753
          END IF
754
            
755
          DO j=1,Element % TYPE % NumberOfNodes
756

757
            IF( DgSolver ) THEN
758
              k1 = Reorder(ElemInds(j))
759
            ELSE
760
              k1 = Reorder(Element % NodeIndexes(j))
761
            END IF
762
                        
763
            IF (.NOT.ALLOCATED(inds)) THEN
764
              ALLOCATE(inds(maxnodes*NumberOfFactors),STAT=istat)
765
            ELSE IF(SIZE(inds)<maxnodes*NumberOfFactors) THEN
766
              DEALLOCATE(inds)
767
              ALLOCATE(inds(maxnodes*NumberOfFactors),STAT=istat)
768
            END IF
769
            IF ( istat /= 0 ) CALL Fatal(Caller,'Memory allocation error fo inds.')
770

771
            cnt = 0
772
            DO n=1,NumberOfFactors
773

774
              Elm => Mesh % Elements( Element % BoundaryInfo % &
775
                  RadiationFactors % Elements(n) )
776

777
              IF( DgSolver ) THEN
778
                CALL DgRadiationIndexes(Elm,Elm % TYPE % NumberOfNodes,ElemInds2,.TRUE.)
779
              END IF
780
              
781
              DO k=1,Elm % TYPE % NumberOfNodes
782
                cnt = cnt + 1
783

784
                IF( DGSolver ) THEN
785
                  inds(cnt) = Reorder(ElemInds2(k))
786
                ELSE
787
                  inds(cnt) = Reorder(Elm % NodeIndexes(k))
788
                END IF
789
              END DO
790
            END DO
791
            
792
            CALL Sort(cnt,inds)
793

794
            CALL List_AddMatrixIndexes(List,k1,cnt,inds)
795
          END DO
796
        END DO
797
      END BLOCK
798
      CALL Info(Caller,'Done Adding radiation matrix',Level=14)
799
    END IF
800

801
    
802
    ! If this is not a DG solver then create permutation considering 
803
    ! nodal, edge, face and bubble dofs. 
804
    !-------------------------------------------------------------------
805
    IF ( .NOT. FoundDG ) THEN
806
      t = 1
807
      DO WHILE( t<=Mesh % NumberOfBulkElements+Mesh % NumberOFBoundaryElements )
808
         Element => Mesh % Elements(t)
809

810
         IF ( PRESENT(Equation) ) THEN
811
           DO WHILE( t<=Mesh % NumberOfBulkElements+Mesh % NumberOfBoundaryElements )
812
             Element => Mesh % Elements(t)
813
             IF ( CheckElementEquation(Model,Element,Equation) ) EXIT
814
             t = t + 1
815
           END DO
816
           IF ( t > Mesh % NumberOfBulkElements+Mesh % NumberOfBoundaryElements ) EXIT
817
         END IF
818

819
         n = Element % NDOFs 
820
         IF( EDOFS > 0 ) n = n + Element % TYPE % NumberOfEdges * EDOFs 
821
         IF( FDOFS > 0 ) n = n + Element % TYPE % NumberOfFaces * FDOFs
822
         IF ( GB ) n = n + Element % BDOFs
823

824
         IF ( n > IndexSize ) THEN
825
            IndexSize = n
826
            IF ( ALLOCATED( Indexes ) ) DEALLOCATE( Indexes )
827
            ALLOCATE( Indexes(n), STAT=istat )
828
            IF( istat /= 0 ) CALL Fatal(Caller,'Allocation error for Indexes of size: '//I2S(n))
829
         END IF
830

831
         n = 0
832
         DOFsPerNode = Element % NDOFs / Element % TYPE % NumberOfNodes
833
         IF (DOFsPerNode > 0) THEN
834
           DO i=1,Element % TYPE % NumberOfNodes
835
             DO j=1,DOFsPerNode
836
               n = n + 1
837
               Indexes(n) = NDOFs * (Element % NodeIndexes(i)-1) + j
838
             END DO
839
           END DO
840
         END IF
841
!         DO i=1,Element % NDOFs
842
!            n = n + 1
843
!            Indexes(n) = Element % NodeIndexes(i)
844
!         END DO
845

846
         IF ( EDOFs > 0 ) THEN
847
            IF ( ASSOCIATED(Element % EdgeIndexes) ) THEN
848
              DO j=1,Element % TYPE % NumberOFEdges
849
                 DO i=1, Mesh % Edges(Element % EdgeIndexes(j)) % BDOFs
850
                   n = n + 1
851
                   Indexes(n) = EDOFs * (Element % EdgeIndexes(j)-1) + i &
852
                                + NDOFs * Mesh % NumberOfNodes
853
                END DO
854
             END DO
855

856
             IF ( GB ) THEN
857
               Edge1 => Mesh % Edges(Element % EdgeIndexes(1))
858
               IF(Element % Type % ElementCode==Edge1 % Type % ElementCode) THEN
859
                 DO i=1, Element % BDOFs
860
                   n = n + 1
861
                   Indexes(n) = EDOFs*(Element % EdgeIndexes(1)-1) + i + &
862
                       NDOFs * Mesh % NumberOfNodes
863
                 END DO
864
               END IF
865
             END IF
866
           END IF
867
         END IF
868

869
         IF ( FDOFS > 0 ) THEN
870
           IF ( ASSOCIATED(Element % FaceIndexes) ) THEN
871
             DO j=1,Element % TYPE % NumberOFFaces
872
               DO i=1, Mesh % Faces(Element % FaceIndexes(j)) % BDOFs
873
                 n = n + 1
874
                 Indexes(n) = FDOFs*(Element % FaceIndexes(j)-1) + i + &
875
                     NDOFs * Mesh % NumberOfNodes + EDOFs*Mesh % NumberOfEdges
876
               END DO
877
             END DO
878

879
             IF ( GB ) THEN
880
               Face1 => Mesh % Faces(Element % FaceIndexes(1))
881
               IF(Element % Type % ElementCode==Face1 % Type % ElementCode) THEN
882
                 DO i=1, Element % BDOFs
883
                   n = n + 1
884
                   Indexes(n) = FDOFs*(Element % FaceIndexes(1)-1) + i + &
885
                       NDOFs * Mesh % NumberOfNodes + EDOFs*Mesh % NumberOfEdges
886
                 END DO
887
               END IF
888
             END IF
889
           END IF
890
         END IF
891

892
         IF ( GB .AND. ASSOCIATED(Element % BubbleIndexes) ) THEN
893
            DO i=1,Element % BDOFs
894
              n = n + 1
895
              Indexes(n) = FDOFs*Mesh % NumberOfFaces + &
896
                   NDOFs * Mesh % NumberOfNodes + EDOFs*Mesh % NumberOfEdges + &
897
                        Element % BubbleIndexes(i)
898
            END DO
899
         END IF
900

901
         IF( PSA ) THEN
902
           CALL PackSortAdd(n,Indexes,Reorder)
903
         ELSE
904
           DO i=1,n
905
             k1 = Reorder(Indexes(i))
906
             IF ( k1 <= 0 ) CYCLE
907
             DO j=1,n
908
               k2 =  Reorder(Indexes(j))
909
               IF ( k2 <= 0 ) CYCLE
910
               Lptr => List_GetMatrixIndex( List,k1,k2 )
911
             END DO
912
           END DO
913
         END IF
914
         t = t + 1
915
      END DO
916

917
      IF ( ALLOCATED( Indexes ) ) DEALLOCATE( Indexes )
918

919
      
920
      DO i=Mesh % NumberOfBulkElements+1, Mesh % NumberOfBulkElements+ &
921
                     Mesh % NumberOfBoundaryElements
922
        IF ( Mesh % Elements(i) % TYPE % ElementCode <  102 .OR. &
923
             Mesh % Elements(i) % TYPE % ElementCode >= 200 ) CYCLE
924

925
        k1 = Reorder( Mesh % Elements(i) % NodeIndexes(1) )
926
        IF ( k1 > 0 ) THEN
927
          DO k=1,Mesh % Elements(i) % TYPE % NumberOFNodes
928
            k2 = Reorder( Mesh % Elements(i) % NodeIndexes(k) )
929
            IF ( k2 > 0 ) THEN
930
              Lptr => List_GetMatrixIndex( List,k1,k2 )
931
              Lptr => List_GetMatrixIndex( List,k2,k1 )
932
            END IF
933
          END DO
934
        END IF
935

936
        ! This is a connection element, make a matrix connection for that 
937
        IF ( Mesh % Elements(i) % TYPE % ElementCode == 102 ) THEN
938
          k2 = Reorder( Mesh % Elements(i) % NodeIndexes(2) )
939
          IF ( k2 > 0 ) Lptr => List_GetMatrixIndex( List,k2,k2 )
940
        END IF
941
      END DO
942

943

944
      ! Add connection from projectors. These are only needed if the projector is treated
945
      ! implicely. For explicit projectors or when using Lagrange coefficients the 
946
      ! connections are not needed. 
947
      !----------------------------------------------------------------------------------
948
      IF( DoProjectors ) THEN
949
        DO This=1,Model % NumberOfBCs
950
          Projector => Model % BCs(This) % PMatrix
951
          IF ( .NOT. ASSOCIATED(Projector) ) CYCLE
952

953
          IF( ListGetLogical( Model % BCs(This) % Values,&
954
              'Periodic BC Explicit',Found)) CYCLE
955
          IF( ListGetLogical( Model % BCs(This) % Values,&
956
              'Periodic BC Use Lagrange Coefficient',Found)) CYCLE
957

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

960
          DO i=1,Projector % NumberOfRows
961
            k = Reorder( Projector % InvPerm(i) )
962
            IF ( k > 0 ) THEN
963
              DO l=Projector % Rows(i),Projector % Rows(i+1)-1
964
                IF ( Projector % Cols(l) <= 0 ) CYCLE
965
                m = Reorder( Projector % Cols(l) )
966
                IF ( m > 0 ) THEN
967
                  Lptr => List_GetMatrixIndex( List,k,m )
968
                  Lptr => List_GetMatrixIndex( List,m,k ) ! keep structure symm.
969
                  CList => List(k) % Head
970
                  DO WHILE( ASSOCIATED( CList ) )
971
                    Lptr => List_GetMatrixIndex( List,m,CList % Index )
972
                    Lptr => List_GetMatrixIndex( List,CList % Index,m ) ! keep structure symm.
973
                    CList => CList % Next
974
                  END DO
975
                END IF
976
              END DO
977
            END IF
978
          END DO
979
        END DO
980
      END IF ! DoProjectors
981

982
    END IF
983
    
984
    Model % TotalMatrixElements = 0
985
    DO i=1,LocalNodes
986
      j = List(i) % Degree
987
      Model % TotalMatrixElements = Model % TotalMatrixElements + j
988
    END DO
989
    
990
      
991
    DoNonZeros = .TRUE.
992
    IF( PRESENT( CalcNonZeros ) ) DoNonZeros = CalcNonZeros
993
    
994
    IF( DoNonZeros ) THEN
995
      ALLOCATE( InvPerm(LocalNodes), STAT=istat )
996
      IF ( istat /= 0 ) CALL Fatal(Caller,'Memory allocation error for InvPerm.')
997
      
998
      InvPerm = 0
999
      k = 0
1000
      DO i=1,SIZE(Reorder)
1001
        IF (Reorder(i)>0) THEN
1002
          k = k + 1
1003
          InvPerm(Reorder(i)) = k
1004
        END IF
1005
      END DO
1006
      
1007
      Model % Rownonzeros = 0
1008
      DO i=1,LocalNodes
1009
        j = List(i) % Degree
1010
        Model % RowNonzeros(InvPerm(i)) = j
1011
      END DO      
1012
      DEALLOCATE( InvPerm ) 
1013
    END IF
1014

1015
  CONTAINS
1016

1017
    ! Takes elemental indexes (for nodes, edges, faces etc.) and using the initial permutation
1018
    ! pack the nonzeros and then sort them and finally add them to the list matrix structure.
1019
    ! Adding the whole row is computationally advantageous to adding just one entry at a time.
1020
    !-----------------------------------------------------------------------------------------
1021
    SUBROUTINE PackSortAdd(n,Ind,Perm)
1022
      INTEGER :: n
1023
      INTEGER :: Ind(:),Perm(:)
1024
      INTEGER :: i,j,m
1025
      
1026
      m = 0
1027
      DO i=1,n
1028
        j = Perm(Ind(i))
1029
        IF(j==0) CYCLE
1030
        m = m+1
1031
        Ind(m) = j
1032
      END DO
1033

1034
      CALL Sort(m,Ind)
1035

1036
      DO i=1,m
1037
        CALL List_AddMatrixIndexes(List,Ind(i),m,Ind)
1038
      END DO
1039

1040
    END SUBROUTINE PackSortAdd
1041
        
1042

1043

1044
#if 0    
1045
    SUBROUTINE DgRadiationIndexes(Element,n,ElemInds)
1046

1047
      TYPE(Element_t), POINTER :: Element
1048
      INTEGER :: n
1049
      INTEGER :: ElemInds(:)
1050

1051
      TYPE(Element_t), POINTER :: Left, Right, Parent
1052
      INTEGER :: i,i1,n1,mat_id
1053
      LOGICAL :: Found
1054
           
1055
      Left => Element % BoundaryInfo % Left
1056
      Right => Element % BoundaryInfo % Right
1057
      
1058
      IF(ASSOCIATED(Left) .AND. ASSOCIATED(Right)) THEN
1059
        DO i=1,2
1060
          IF(i==1) THEN
1061
            Parent => Left
1062
          ELSE
1063
            Parent => Right
1064
          END IF
1065
          
1066
          mat_id = ListGetInteger( CurrentModel % Bodies(Parent % BodyId) % Values, &
1067
              'Material', Found, minv=1,maxv=CurrentModel % NumberOfMaterials )
1068
          IF( ListCheckPresent( CurrentModel % Materials(mat_id) % Values,'Emissivity') ) EXIT
1069

1070
          IF( i==2) THEN
1071
            CALL Fatal(Caller,'DG boundary parents should have emissivity!')
1072
          END IF          
1073
        END DO                
1074
      ELSE IF(ASSOCIATED(Left)) THEN
1075
        Parent => Left
1076
      ELSE IF(ASSOCIATED(Right)) THEN
1077
        Parent => Right
1078
      ELSE
1079
        CALL Fatal(Caller,'DG boundary should have parents!')
1080
      END IF
1081

1082
      n1 = Parent % TYPE % NumberOfNodes
1083
      DO i=1,n
1084
        DO i1 = 1, n1
1085
          IF( Parent % NodeIndexes( i1 ) == Element % NodeIndexes(i) ) THEN
1086
            ElemInds(i) = Parent % DGIndexes( i1 )
1087
            EXIT
1088
          END IF
1089
        END DO
1090
      END DO
1091
      
1092
    END SUBROUTINE DgRadiationIndexes
1093
#endif
1094
    
1095
!------------------------------------------------------------------------------
1096
  END SUBROUTINE MakeListMatrix
1097
!------------------------------------------------------------------------------
1098

1099

1100
  ! Pick thet correct indexes for radition when using discontinuous Galerkin.
1101
  ! In internal BCs we expect to find 'emissivity' given on either side.
1102
  !--------------------------------------------------------------------------
1103
  SUBROUTINE DgRadiationIndexes(Element,n,ElemInds,DiffuseGray)
1104

1105
    TYPE(Element_t), POINTER :: Element
1106
    INTEGER :: n
1107
    INTEGER :: ElemInds(:)
1108
    LOGICAL :: DiffuseGray
1109

1110
    TYPE(Element_t), POINTER :: Left, Right, Parent
1111
    INTEGER :: i,i1,n1,mat_id
1112
    LOGICAL :: Found
1113

1114
    Left => Element % BoundaryInfo % Left
1115
    Right => Element % BoundaryInfo % Right
1116

1117
    IF(ASSOCIATED(Left) .AND. ASSOCIATED(Right)) THEN
1118

1119
      IF( DiffuseGray ) THEN
1120
        DO i=1,2
1121
          IF(i==1) THEN
1122
            Parent => Left
1123
          ELSE
1124
            Parent => Right
1125
          END IF
1126

1127
          mat_id = ListGetInteger( CurrentModel % Bodies(Parent % BodyId) % Values, &
1128
              'Material', Found, minv=1,maxv=CurrentModel % NumberOfMaterials )
1129
          IF(.NOT. Found ) THEN
1130
            CALL Fatal('DGRadiationIndexes','Body '//I2S(Parent % BodyId)//' has no Material associated!')
1131
          END IF
1132
          IF( ListCheckPresent( CurrentModel % Materials(mat_id) % Values,'Emissivity') ) EXIT
1133

1134
          IF( i==2) THEN
1135
            CALL Fatal('DGRadiationIndexes','DG boundary parents should have emissivity!')
1136
          END IF
1137
        END DO
1138
      ELSE
1139
        IF( Left % BodyId > Right % BodyId ) THEN
1140
          Parent => Left
1141
        ELSE
1142
          Parent => Right
1143
        END IF
1144
      END IF
1145
    ELSE IF(ASSOCIATED(Left)) THEN
1146
      Parent => Left
1147
    ELSE IF(ASSOCIATED(Right)) THEN
1148
      Parent => Right
1149
    ELSE
1150
      CALL Fatal('DGRadiationIndexes','DG boundary should have parents!')
1151
    END IF
1152

1153
    n1 = Parent % TYPE % NumberOfNodes
1154
    DO i=1,n
1155
      DO i1 = 1, n1
1156
        IF( Parent % NodeIndexes( i1 ) == Element % NodeIndexes(i) ) THEN
1157
          ElemInds(i) = Parent % DGIndexes( i1 )
1158
          EXIT
1159
        END IF
1160
      END DO
1161
    END DO
1162

1163
  END SUBROUTINE DgRadiationIndexes
1164

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

1208
    CALL Info(Caller,'Creating list matrix',Level=14)
1209

1210
    GB = .FALSE.
1211
    IF ( PRESENT(GlobalBubbles) ) GB = GlobalBubbles
1212

1213
    ! No colouring equals a single colour
1214
    NumberOfMeshColours = 1
1215
    IF (ASSOCIATED(Solver % ColourIndexList)) THEN
1216
       NumberOfMeshColours = Solver % ColourIndexList % n
1217
       ! If boundary mesh exists, it has been coloured as well
1218
       IF (ASSOCIATED(Solver % BoundaryColourIndexList)) &
1219
            NumberOfMeshColours = NumberOfMeshColours + Solver % BoundaryColourIndexList % n
1220
    END IF
1221

1222
    nthr=1
1223
    !$ nthr = omp_get_max_threads()
1224
    NeedLocking = (NumberOfMeshColours == 1) .AND. (nthr > 1)
1225
    CALL ListMatrixArray_Allocate(List, LocalNodes, Atomic=NeedLocking)
1226
    
1227
    NDOFs = Mesh % MaxNDOFs
1228
    BDOFs = Mesh % MaxBDOFs
1229
    EDOFs = Mesh % MaxEdgeDOFs
1230
    FDOFs = Mesh % MaxFaceDOFs
1231
    IF( PRESENT( NodalDofsOnly ) ) THEN
1232
      IF( NodalDofsOnly ) THEN
1233
        EDOFS = 0
1234
        FDOFS = 0
1235
      END IF
1236
    END IF
1237

1238
    IF( PRESENT( ProjectorDofs ) ) THEN
1239
      DoProjectors = ProjectorDofs 
1240
    ELSE
1241
      DoProjectors = .TRUE.
1242
    END IF
1243

1244
    IF( EDOFS > 0 .AND. .NOT. ASSOCIATED(Mesh % Edges) ) THEN
1245
      CALL Warn(Caller,'Edge dofs requested but not edges exist in mesh!')
1246
      EDOFS = 0
1247
    END IF
1248

1249
    IF( FDOFS > 0 .AND. .NOT. ASSOCIATED(Mesh % Faces) ) THEN
1250
      CALL Warn(Caller,'Face dofs requested but not faces exist in mesh!')
1251
      FDOFS = 0
1252
    END IF
1253

1254
    ! Create the permutation for the Discontinuous Galerkin solver 
1255
    !-------------------------------------------------------------------
1256
    FoundDG = .FALSE.
1257
    IF ( DGSolver ) THEN
1258
       IndexSize = 128
1259
       ALLOCATE( Indexes(IndexSize), STAT=istat )
1260
       IF( istat /= 0 ) CALL Fatal(Caller,'Allocation error for Indexes')
1261

1262
       ! TODO: Add multithreading
1263
       DO t=1,Mesh % NumberOfEdges
1264
         n = 0
1265
         Elm => Mesh % Edges(t) % BoundaryInfo % Left
1266
         IF ( ASSOCIATED( Elm ) ) THEN
1267
             IF ( CheckElementEquation(Model,Elm,Equation) ) THEN
1268
                FoundDG = FoundDG .OR. Elm % DGDOFs > 0
1269
                DO j=1,Elm % DGDOFs
1270
                   n = n + 1
1271
                   Indexes(n) = Elm % DGIndexes(j)
1272
                END DO
1273
             END IF
1274
         END IF
1275

1276
         Elm => Mesh % Edges(t) % BoundaryInfo %  Right
1277
         IF ( ASSOCIATED( Elm ) ) THEN
1278
             IF ( CheckElementEquation(Model,Elm,Equation) ) THEN
1279
                FoundDG = FoundDG .OR. Elm % DGDOFs > 0
1280
                DO j=1,Elm % DGDOFs
1281
                   n = n + 1
1282
                   Indexes(n) = Elm % DGIndexes(j)
1283
                END DO
1284
             END IF
1285
         END IF
1286

1287
         DO i=1,n
1288
            k1 = Reorder(Indexes(i))
1289
            IF ( k1 <= 0 ) CYCLE
1290
            DO j=1,n
1291
              k2 = Reorder(Indexes(j))
1292
              IF ( k2 <= 0 ) CYCLE
1293
              
1294
              CALL ListMatrixArray_AddEntry(List, k1, k2)
1295
            END DO
1296
         END DO
1297
      END DO
1298
      DO t=1,Mesh % NumberOfFaces
1299
         n = 0
1300
         Elm => Mesh % Faces(t) % BoundaryInfo % Left
1301
         IF ( ASSOCIATED( Elm ) ) THEN
1302
             IF ( CheckElementEquation(Model,Elm,Equation) ) THEN
1303
                FoundDG = FoundDG .OR. Elm % DGDOFs > 0
1304
                DO j=1,Elm % DGDOFs
1305
                   n = n + 1
1306
                   Indexes(n) = Elm % DGIndexes(j)
1307
                END DO
1308
             END IF
1309
         END IF
1310

1311
         Elm => Mesh % Faces(t) % BoundaryInfo %  Right
1312
         IF ( ASSOCIATED( Elm ) ) THEN
1313
             IF ( CheckElementEquation(Model,Elm,Equation) ) THEN
1314
                FoundDG = FoundDG .OR. Elm % DGDOFs > 0
1315
                DO j=1,Elm % DGDOFs
1316
                   n = n + 1
1317
                   Indexes(n) = Elm % DGIndexes(j)
1318
                END DO
1319
             END IF
1320
         END IF
1321

1322
         DO i=1,n
1323
            k1 = Reorder(Indexes(i))
1324
            IF ( k1 <= 0 ) CYCLE
1325
            DO j=1,n
1326
              k2 = Reorder(Indexes(j))
1327
              IF ( k2 <= 0 ) CYCLE
1328
              
1329
              CALL ListMatrixArray_AddEntry(List, k1, k2)
1330
            END DO
1331
         END DO
1332
      END DO
1333
      DEALLOCATE(Indexes)
1334
    END IF
1335

1336
    ! If this is not a GD solver then create permutation considering 
1337
    ! nodal, edge, face and bubble dofs. 
1338
    !-------------------------------------------------------------------
1339
    IF (.NOT. FoundDG) THEN
1340
      
1341
      !$OMP PARALLEL &
1342
      !$OMP SHARED(LocalNodes, List, Equation, NDOFs, EDOFS, FDOFS, GB, &
1343
      !$OMP        Reorder, Model, Mesh, NumberOfMeshColours, CurrentColourStart, NeedLocking, &
1344
      !$OMP        CurrentColourEnd, CurrentColourList, ElementsList, Solver, BoundaryColour) &
1345
      !$OMP PRIVATE(Element, Indexes, istat, IndexSize, IndexReord, &
1346
      !$OMP         IPerm, n, i, j, nReord, k1, k2, Lptr, DOFsPerNode, &
1347
      !$OMP         CurrentColour, Edge1, Face1) &
1348
      !$OMP DEFAULT(NONE)
1349

1350
      IndexSize = 0
1351

1352
      ! Loop over mesh colours
1353
      DO CurrentColour=1,NumberOfMeshColours
1354

1355
         !$OMP SINGLE
1356
         ! Test if matrix is actually coloured or not
1357
         IF (NumberOfMeshColours == 1) THEN
1358
           CurrentColourList => NULL()
1359
           ElementsList => Mesh % Elements(1:Mesh % NumberOfBulkElements+Mesh % NumberOFBoundaryElements)
1360
           CurrentColourStart = 1
1361
           CurrentColourEnd = Mesh % NumberOfBulkElements+Mesh % NumberOFBoundaryElements
1362
         ELSE IF (CurrentColour <= Solver % ColourIndexList % n) THEN
1363
           CALL Info(Caller,'ListMatrix add colour: '//I2S(CurrentColour),Level=10)
1364
           CurrentColourList => Solver % ColourIndexList
1365
           ElementsList => Mesh % Elements(1:Mesh % NumberOfBulkElements)
1366
           CurrentColourStart = CurrentColourList % Ptr(CurrentColour)
1367
           CurrentColourEnd = CurrentColourList % Ptr(CurrentColour+1)-1
1368
         ELSE
1369
           BoundaryColour = CurrentColour-Solver % ColourIndexList % n
1370
           CALL Info(Caller,'ListMatrix add boundary colour: '//I2S(BoundaryColour),Level=10)
1371

1372
           CurrentColourList => Solver % BoundaryColourIndexList
1373
           ! Boundary elements are stored after bulk elements in Mesh
1374
           ElementsList => Mesh % Elements(Mesh % NumberOfBulkElements+1:&
1375
                 Mesh % NumberOfBulkElements + Mesh % NumberOfBoundaryElements)
1376
           CurrentColourStart = CurrentColourList % Ptr(BoundaryColour)
1377
           CurrentColourEnd = CurrentColourList % Ptr(BoundaryColour+1)-1
1378
         END IF
1379
         !$OMP END SINGLE
1380
         
1381
         !$OMP DO
1382
         DO t=CurrentColourStart,CurrentColourEnd
1383
            IF (ASSOCIATED(CurrentColourList)) THEN
1384
               Element => ElementsList(CurrentColourList % ind(t))
1385
            ELSE
1386
               Element => ElementsList(t)
1387
            END IF
1388

1389

1390
            IF ( PRESENT(Equation) ) THEN
1391
               IF ( .NOT. CheckElementEquation(Model,Element,Equation) ) CYCLE
1392
            END IF
1393

1394
            n = Element % NDOFs 
1395
            IF( EDOFS > 0 ) n = n + Element % TYPE % NumberOfEdges * EDOFs 
1396
            IF( FDOFS > 0 ) n = n + Element % TYPE % NumberOfFaces * FDOFs
1397
            IF ( GB ) n = n + Element % BDOFs
1398
            
1399
            ! Reallocate index permutation if needed
1400
            IF ( n > IndexSize ) THEN
1401
               IF (ALLOCATED(Indexes)) DEALLOCATE(Indexes, IndexReord, IPerm)
1402
               
1403
               IndexSize = MAX(MAX(128, IndexSize*2), n)
1404
               ALLOCATE(Indexes(IndexSize), &
1405
                   IndexReord(IndexSize), &
1406
                   IPerm(IndexSize), STAT=istat )
1407
               IF( istat /= 0 ) CALL Fatal(Caller,'Allocation error for Indexes of size: '//I2S(n))
1408
            END IF
1409
            
1410
            n = 0
1411
            DOFsPerNode = Element % NDOFs / Element % TYPE % NumberOfNodes
1412
            IF (DOFsPerNode > 0) THEN
1413
              DO i=1,Element % TYPE % NumberOfNodes
1414
                DO j=1,DOFsPerNode
1415
                  n = n + 1
1416
                  Indexes(n) = NDOFs * (Element % NodeIndexes(i)-1) + j
1417
                END DO
1418
              END DO
1419
            END IF
1420
!            DO i=1,Element % NDOFs
1421
!               n = n + 1
1422
!               Indexes(n) = Element % NodeIndexes(i)
1423
!            END DO
1424
            
1425
            IF ( EDOFs > 0 ) THEN
1426
               IF ( ASSOCIATED(Element % EdgeIndexes) ) THEN
1427
                  DO j=1,Element % TYPE % NumberOFEdges
1428
                     DO i=1, Mesh % Edges(Element % EdgeIndexes(j)) % BDOFs
1429
                        n = n + 1
1430
                        Indexes(n) = EDOFs * (Element % EdgeIndexes(j)-1) + i &
1431
                             + NDOFs * Mesh % NumberOfNodes
1432
                     END DO
1433
                  END DO
1434
                  IF ( GB .AND. ASSOCIATED(Element % BoundaryInfo)) THEN
1435
                    Edge1 => Mesh % Edges(Element % EdgeIndexes(1))
1436
                    IF(Element % Type % ElementCode==Edge1 % Type % ElementCode) THEN
1437
                      DO i=1, Element % BDOFs
1438
                        n = n + 1
1439
                        Indexes(n) = EDOFs*(Element % EdgeIndexes(1)-1) + i + &
1440
                            NDOFs * Mesh % NumberOfNodes
1441
                      END DO
1442
                    END IF
1443
                  END IF
1444
               END IF
1445
            END IF
1446
            
1447
            IF ( FDOFS > 0 ) THEN
1448
               IF ( ASSOCIATED(Element % FaceIndexes) ) THEN
1449
                  DO j=1,Element % TYPE % NumberOFFaces
1450
                     DO i=1, Mesh % Faces(Element % FaceIndexes(j)) % BDOFs
1451
                        n = n + 1
1452
                        Indexes(n) = FDOFs*(Element % FaceIndexes(j)-1) + i + &
1453
                             NDOFs * Mesh % NumberOfNodes + EDOFs*Mesh % NumberOfEdges
1454
                     END DO
1455
                  END DO
1456
               END IF
1457

1458
               IF ( GB.AND. ASSOCIATED(Element % BoundaryInfo) ) THEN
1459
                 Face1 => Mesh % Faces(Element % FaceIndexes(1))
1460
                 IF(Element % Type % ElementCode==Face1 % Type % ElementCode) THEN
1461
                   DO i=1, Element % BDOFs
1462
                     n = n + 1
1463
                     Indexes(n) = FDOFs*(Element % FaceIndexes(1)-1) + i + &
1464
                         NDOFs * Mesh % NumberOfNodes + EDOFs*Mesh % NumberOfEdges
1465
                   END DO
1466
                 END IF
1467
               END IF
1468
            END IF
1469
            
1470
            IF ( GB .AND. ASSOCIATED(Element % BubbleIndexes) ) THEN
1471
               DO i=1,Element % BDOFs
1472
                  n = n + 1
1473
                  Indexes(n) = FDOFs*Mesh % NumberOfFaces + &
1474
                       NDOFs * Mesh % NumberOfNodes + EDOFs*Mesh % NumberOfEdges + &
1475
                       Element % BubbleIndexes(i)
1476
               END DO
1477
            END IF
1478

1479
            nReord=0
1480
            DO i=1,n
1481
               IF (Reorder(Indexes(i)) > 0) THEN
1482
                  nReord=nReord+1
1483
                  IndexReord(nReord)=Reorder(Indexes(i))
1484
               END IF
1485
            END DO
1486
            CALL InsertionSort(nReord, IndexReord, IPerm)
1487

1488
            DO i=1,nReord
1489
               k1 = IndexReord(IPerm(i))
1490
             
1491
               ! Bulk add all sorted nonzero indices to a matrix row in one go
1492
               CALL ListMatrixArray_AddEntries(List, k1, nReord, IndexReord, IPerm, NeedLocking)
1493
            END DO
1494
         END DO
1495
         !$OMP END DO
1496
      END DO ! Loop over mesh colours
1497

1498
      IF (ALLOCATED(Indexes)) DEALLOCATE(Indexes, IndexReord, IPerm)
1499
      !$OMP END PARALLEL
1500

1501
!
1502
!     Diffuse gray radiation condition:
1503
!     ---------------------------------
1504
      Radiation = ListGetLogical( Solver % Values, 'Radiation Solver', Found )
1505
      IF( Radiation ) THEN        
1506
        Radiation = .FALSE.
1507
        DO i=1,Model % NumberOfBCs
1508
          RadiationFlag = ListGetString( Model % BCs(i) % Values, 'Radiation', GotIt )
1509
          IF (GotIt) THEN
1510
            IF ( RadiationFlag == 'diffuse gray' ) THEN
1511
              Radiation = .TRUE.
1512
              EXIT
1513
            END IF
1514
          END IF
1515
        END DO
1516
      END IF
1517
      
1518
      IF ( Radiation ) THEN
1519
        CALL Info(Caller,'Adding radiation matrix entries',Level=14)
1520
        DO i = Mesh % NumberOfBulkElements+1, &
1521
          Mesh % NumberOfBulkElements + Mesh % NumberOfBoundaryElements
1522

1523
          Element => Mesh % Elements(i)
1524
          IF ( ASSOCIATED(Element % BoundaryInfo % RadiationFactors) ) THEN
1525
             DO j=1,Element % TYPE % NumberOfNodes
1526
                k1 = Reorder(Element % NodeIndexes(j))
1527

1528
                NumberOfFactors = Element % BoundaryInfo % &
1529
                  RadiationFactors % NumberOfImplicitFactors
1530

1531
                DO n=1,NumberOfFactors
1532

1533
                  Elm => Mesh % Elements( Element % BoundaryInfo % &
1534
                              RadiationFactors % Elements(n) )
1535

1536
                  DO k=1,Elm % TYPE % NumberOfNodes
1537
                     k2 = Reorder( Elm % NodeIndexes(k) )
1538
                     CALL ListMatrixArray_AddEntry(List, k1, k2)
1539
                  END DO
1540
                END DO
1541
             END DO
1542
          END IF
1543
        END DO
1544
        CALL Info(Caller,'Done Adding radiation matrix',Level=14)
1545
      END IF
1546

1547
      ! TODO: Add multithreading
1548
      DO i=Mesh % NumberOfBulkElements+1, Mesh % NumberOfBulkElements+ &
1549
                     Mesh % NumberOfBoundaryElements
1550
        IF ( Mesh % Elements(i) % TYPE % ElementCode <  102 .OR. &
1551
             Mesh % Elements(i) % TYPE % ElementCode >= 200 ) CYCLE
1552

1553
        k1 = Reorder( Mesh % Elements(i) % NodeIndexes(1) )
1554
        IF ( k1 > 0 ) THEN
1555
          DO k=1,Mesh % Elements(i) % TYPE % NumberOFNodes
1556
            k2 = Reorder( Mesh % Elements(i) % NodeIndexes(k) )
1557
            IF ( k2 > 0 ) THEN
1558
              CALL ListMatrixArray_AddEntry(List, k1, k2)
1559
              CALL ListMatrixArray_AddEntry(List, k2, k1)
1560
            END IF
1561
          END DO
1562
        END IF
1563

1564
        ! This is a connection element, make a matrix connection for that 
1565
        IF ( Mesh % Elements(i) % TYPE % ElementCode == 102 ) THEN
1566
          k2 = Reorder( Mesh % Elements(i) % NodeIndexes(2) )
1567
          IF ( k2 > 0 ) CALL ListMatrixArray_AddEntry(List, k2, k2)
1568
        END IF
1569
      END DO
1570

1571

1572
      ! Add connection from projectors. These are only needed if the projector is treated
1573
      ! implicely. For explicit projectors or when using Lagrange coefficients the 
1574
      ! connections are not needed. 
1575
      !----------------------------------------------------------------------------------
1576
      IF( DoProjectors ) THEN
1577
        DO This=1,Model % NumberOfBCs
1578
          Projector => Model % BCs(This) % PMatrix
1579
          IF ( .NOT. ASSOCIATED(Projector) ) CYCLE
1580

1581
          IF( ListGetLogical( Model % BCs(This) % Values,&
1582
              'Periodic BC Explicit',Found)) CYCLE
1583
          IF( ListGetLogical( Model % BCs(This) % Values,&
1584
              'Periodic BC Use Lagrange Coefficient',Found)) CYCLE
1585

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

1588
          ! TODO: Add multithreading
1589
          DO i=1,Projector % NumberOfRows
1590
            k = Reorder( Projector % InvPerm(i) )
1591
            IF ( k > 0 ) THEN
1592
              DO l=Projector % Rows(i),Projector % Rows(i+1)-1
1593
                IF ( Projector % Cols(l) <= 0 ) CYCLE
1594
                m = Reorder( Projector % Cols(l) )
1595
                IF ( m > 0 ) THEN
1596
                  CALL ListMatrixArray_AddEntry(List, k, m)
1597
                  CALL ListMatrixArray_AddEntry(List, m, k)
1598
                  CList => List % Rows(k) % Head
1599
                  DO WHILE( ASSOCIATED( CList ) )
1600
                    CALL ListMatrixArray_AddEntry(List, m, CList % Index)
1601
                    CALL ListMatrixArray_AddEntry(List, CList % Index, m)
1602
                    CList => CList % Next
1603
                  END DO
1604
                END IF
1605
              END DO
1606
            END IF
1607
          END DO
1608
        END DO
1609
      END IF ! DoProjectors
1610

1611
    END IF ! (.NOT. FoundDG)
1612

1613

1614

1615
    Model % TotalMatrixElements = 0
1616
    DO i=1,LocalNodes
1617
      j = List % Rows(i) % Degree
1618
      Model % TotalMatrixElements = Model % TotalMatrixElements + j
1619
    END DO
1620
    
1621
      
1622
    DoNonZeros = .TRUE.
1623
    IF( PRESENT( CalcNonZeros ) ) DoNonZeros = CalcNonZeros
1624
    
1625
    IF( DoNonZeros ) THEN
1626
      ALLOCATE( InvPerm(LocalNodes) )
1627
      InvPerm = 0
1628
      k = 0
1629
      DO i=1,SIZE(Reorder)
1630
        IF (Reorder(i)>0) THEN
1631
          k = k + 1
1632
          InvPerm(Reorder(i)) = k
1633
        END IF
1634
      END DO
1635
      
1636
      Model % Rownonzeros = 0
1637
      DO i=1,LocalNodes
1638
        j = List % Rows(i) % Degree
1639
        Model % RowNonzeros(InvPerm(i)) = j
1640
      END DO      
1641
      DEALLOCATE( InvPerm ) 
1642
    END IF
1643

1644
!------------------------------------------------------------------------------
1645
  END SUBROUTINE MakeListMatrixArray
1646
!------------------------------------------------------------------------------
1647

1648
  
1649

1650
!------------------------------------------------------------------------------
1651
!>    Initialize a CRS format matrix to the effect that it will be ready to
1652
!>    accept values when CRS_GlueLocalMatrix is called (build up the index
1653
!>    tables of a CRS format matrix)....
1654
!------------------------------------------------------------------------------
1655
  SUBROUTINE InitializeMatrix( Matrix, n, List, DOFs, Reorder, InvInitialReorder )
1656
!------------------------------------------------------------------------------
1657
    INTEGER :: DOFs, n
1658
    TYPE(Matrix_t),POINTER :: Matrix
1659
    TYPE(ListMatrix_t) :: List(:)
1660
    INTEGER, OPTIONAL :: Reorder(:), InvInitialReorder(:)
1661
!------------------------------------------------------------------------------
1662
    TYPE(ListMatrixEntry_t), POINTER :: Clist
1663
    INTEGER :: i,j,k,l,m,k1,k2
1664
    INTEGER, POINTER CONTIG :: Rows(:), Cols(:)
1665
    LOGICAL :: DoReorder 
1666

1667
    
1668
    Rows => Matrix % Rows
1669
    Cols => Matrix % Cols
1670
    
1671
    DoReorder = .FALSE.
1672
    IF( PRESENT( Reorder ) ) THEN
1673
      IF( .NOT. PRESENT( InvInitialReorder ) ) THEN
1674
        CALL Fatal('InitializeMatrix','Need both old and new numbering!')
1675
      END IF
1676
      DoReorder = .TRUE.
1677
    END IF
1678
    
1679
    
1680
    IF( DoReorder ) THEN
1681
      ! In case of reordering we need to compute the row offset in advance
1682
      Rows(1) = 1
1683
      DO i=1,n       
1684
        DO l=1,DOFs
1685
          j = Reorder( InvInitialReorder(i) )
1686
          k1 = DOFs * (j-1) + l
1687
          Rows(k1+1) = Dofs * List(i) % Degree
1688
        END DO
1689
      END DO
1690
      DO i=1,Dofs*n
1691
        Rows(i+1) = Rows(i) + Rows(i+1)
1692
      END DO
1693
      
1694
      !$OMP PARALLEL DO SHARED(Rows, Cols, List, n, DOFs, Reorder, InvInitialReorder) &
1695
      !$OMP PRIVATE(CList, l, j, k1, k2, k, m) &
1696
      !$OMP DEFAULT(NONE)
1697
      DO i=1,n
1698
        DO l=1,DOFs
1699
          CList => List(i) % Head
1700
          j = Reorder( InvInitialReorder(i) )
1701
          k1 = DOFs * (j-1) + l
1702
          k2 = Rows(k1)-1
1703
          
1704
          DO WHILE( ASSOCIATED( CList ) )
1705
            k = Reorder( InvInitialReorder(Clist % INDEX))
1706
            k = DOFs*(k-1)
1707
            DO m=k+1,k+DOFs
1708
              k2 = k2+1
1709
              Cols(k2) = m
1710
            END DO
1711
            CList => Clist % Next
1712
          END DO
1713
        END DO
1714
      END DO
1715
      !$OMP END PARALLEL DO      
1716
    ELSE
1717
      Rows(1) = 1
1718
      DO i=1,n       
1719
        DO l=1,DOFs
1720
          k1 = DOFs * (i-1) + l
1721
          Rows(k1+1) = Dofs * List(i) % Degree
1722
        END DO
1723
      END DO
1724
      DO i=1,Dofs*n
1725
        Rows(i+1) = Rows(i) + Rows(i+1)
1726
      END DO
1727

1728
      ! If there is no renumbering then the reordering is one-to-one mapping
1729
      !$OMP PARALLEL DO SHARED(Rows, Cols, List, n, DOFs) &
1730
      !$OMP PRIVATE(CList, l, j, k1, k2, k, m) &
1731
      !$OMP DEFAULT(NONE)
1732

1733
      DO i=1,n
1734
        DO l=1,DOFs
1735
          CList => List(i) % Head
1736
          j = i
1737
          k1 = DOFs * (j-1) + l
1738
          k2 = Rows(k1)-1
1739
          
1740
          DO WHILE( ASSOCIATED( CList ) )
1741
            k = Clist % index 
1742
            k = DOFs*(k-1)
1743
            DO m=k+1,k+DOFs
1744
              k2 = k2+1
1745
              Cols(k2) = m
1746
            END DO
1747
            CList => Clist % Next
1748
          END DO
1749
        END DO
1750
      END DO
1751
      !$OMP END PARALLEL DO      
1752
    END IF
1753
    
1754
    IF ( Matrix % FORMAT == MATRIX_CRS ) CALL CRS_SortMatrix( Matrix )
1755
!------------------------------------------------------------------------------
1756
  END SUBROUTINE InitializeMatrix
1757
!------------------------------------------------------------------------------
1758

1759

1760
!------------------------------------------------------------------------------
1761
   FUNCTION CreateMatrix( Model, Solver, Mesh, Perm, DOFs, MatrixFormat, &
1762
          OptimizeBW, Equation, DGSolver, GlobalBubbles, &
1763
          NodalDofsOnly, ProjectorDofs, ThreadedStartup, &
1764
          UseGivenPerm ) RESULT(Matrix)
1765
!------------------------------------------------------------------------------
1766
     IMPLICIT NONE
1767
     TYPE(Model_t) :: Model
1768
     TYPE(Mesh_t) :: Mesh
1769
     TYPE(Solver_t), TARGET :: Solver
1770
     INTEGER :: DOFs, MatrixFormat
1771
     INTEGER, TARGET :: Perm(:)
1772
     LOGICAL :: OptimizeBW
1773
     LOGICAL, OPTIONAL :: DGSolver, GlobalBubbles
1774
     LOGICAL, OPTIONAL :: NodalDofsOnly, ProjectorDofs
1775
     LOGICAL, OPTIONAL :: ThreadedStartup
1776
     LOGICAL, OPTIONAL :: USeGivenPerm
1777

1778
     CHARACTER(LEN=*), OPTIONAL :: Equation
1779

1780
     TYPE(Matrix_t),POINTER :: Matrix
1781
!------------------------------------------------------------------------------
1782
     TYPE(ListMatrix_t), POINTER :: ListMatrix(:)
1783
     TYPE(ListMatrixArray_t) :: ListMatrixArray
1784
     TYPE(Matrix_t), POINTER :: A
1785
     TYPE(Element_t), POINTER :: Element
1786
     TYPE(ListMatrixEntry_t), POINTER :: CList
1787
     CHARACTER(:), ALLOCATABLE :: Eq,str
1788
     LOGICAL :: GotIt, DG, GB, UseOptimized, Found, UseGiven
1789
     INTEGER i,j,k,l,k1,t,n, p,m, minEdgeDOFs, maxEdgeDOFs, &
1790
           minFaceDOFs, maxFaceDOFs, BDOFs, cols, istat, &
1791
           NDOFs
1792
     INTEGER, POINTER :: Ivals(:)
1793
     INTEGER, ALLOCATABLE, SAVE :: InvInitialReorder(:)
1794
     INTEGER :: nthr
1795
     LOGICAL :: UseThreads
1796
     LOGICAL, ALLOCATABLE :: ConstrainedNode(:)
1797
     CHARACTER(*), PARAMETER :: Caller = 'CreateMatrix'
1798
   
1799
!------------------------------------------------------------------------------
1800

1801
     NULLIFY( Matrix )
1802

1803
     DG = .FALSE.
1804
     IF ( PRESENT(DGSolver) )  DG = DGSolver
1805

1806
     GB = .FALSE.
1807
     IF ( PRESENT(GlobalBubbles) ) GB = GlobalBubbles
1808

1809
     UseGiven = .FALSE.
1810
     IF( PRESENT( UseGivenPerm ) ) THEN
1811
       IF( UseGivenPerm ) THEN
1812
         CALL Info(Caller,'Using given Perm table to create the matrix',Level=6)
1813
         UseGiven = .TRUE.
1814
         OptimizeBW = .FALSE.
1815
       END IF
1816
     END IF
1817
       
1818
     IF( OptimizeBW ) THEN
1819
       IF( ListGetLogical( Solver % Values,'DG Reduced Basis',Found ) ) THEN
1820
         CALL Info(Caller,'Suppressing bandwidth optimization for discontinuous bodies',Level=8)
1821
         OptimizeBW = .FALSE.
1822
       END IF
1823
       IF( ListGetLogical( Solver % Values,'Apply Conforming BCs',Found ) ) THEN
1824
         CALL Info(Caller,'Suppressing bandwidth optimization for conforming bcs',Level=8)
1825
         OptimizeBW = .FALSE.
1826
       END IF
1827
     END IF
1828

1829
     
1830
     UseThreads = .FALSE.
1831
     nthr = 1
1832
     IF ( PRESENT(ThreadedStartup) ) THEN
1833
#ifdef _OPENMP
1834
       IF (ThreadedStartup) THEN
1835
          UseThreads = .TRUE.
1836
          nthr = omp_get_max_threads()
1837
       END IF
1838
#endif
1839
     END IF
1840

1841
     minEdgeDOFs = HUGE(minEdgeDOFs)
1842
     maxEdgeDOFs = 0
1843
     minFaceDOFs = HUGE(minFaceDOFs)
1844
     maxFaceDOFs = 0
1845
     BDOFs = 0
1846
     NDOFs = 0
1847

1848
     !$OMP PARALLEL SHARED(Mesh, GB) PRIVATE(j, Element) &
1849
     !$OMP          REDUCTION(min:minEdgeDOFs) REDUCTION(max:maxEdgeDOFs) &
1850
     !$OMP          REDUCTION(min:minFaceDOFs) REDUCTION(max:maxFaceDOFs) &
1851
     !$OMP          REDUCTION(max:BDOFs) &
1852
     !$OMP          DEFAULT(NONE) NUM_THREADS(nthr)
1853
     
1854
     !$OMP DO 
1855
     DO i=1,Mesh % NumberOfEdges
1856
        minEdgeDOFs = MIN( minEdgeDOFs, Mesh % Edges(i) % BDOFs )
1857
        maxEdgeDOFs = MAX( maxEdgeDOFs, Mesh % Edges(i) % BDOFs )
1858
     END DO
1859
     !$OMP END DO NOWAIT
1860
     !$OMP DO
1861
     DO i=1,Mesh % NumberOfFaces
1862
        minFaceDOFs = MIN( minFaceDOFs, Mesh % Faces(i) % BDOFs )
1863
        maxFaceDOFs = MAX( maxFaceDOFs, Mesh % Faces(i) % BDOFs )
1864
     END DO
1865
     !$OMP END DO NOWAIT
1866

1867
     IF ( GB ) THEN
1868
       DO i=1,Mesh % NumberOfBoundaryElements
1869
          j = i + Mesh % NumberOfBulkElements
1870
          Element => Mesh % Elements(j)
1871
          IF(Element % Type % ElementCode >= 300) THEN
1872
            minFaceDOFs = MIN( minFaceDOFs, Element % BDOFs )
1873
            maxFaceDOFs = MAX( maxFaceDOFs, Element % BDOFs )
1874
          ELSE
1875
            minEdgeDOFs = MIN( minEdgeDOFs, Element % BDOFs )
1876
            maxEdgeDOFs = MAX( maxEdgeDOFs, Element % BDOFs )
1877
          END IF
1878
       END DO
1879
     END IF
1880

1881
     !$OMP DO
1882
     DO i=1,Mesh % NumberOfBulkElements
1883
        BDOFs = MAX( BDOFs, Mesh % Elements(i) % BDOFs )
1884
     END DO
1885
     !$OMP END DO NOWAIT
1886
     !$OMP END PARALLEL
1887
     DO i=1,Mesh % NumberOfBulkElements
1888
        NDOFs = MAX( NDOFs, Mesh % Elements(i) % NDOFs )
1889
     END DO
1890
     
1891
     Mesh % MaxEdgeDOFs = maxEdgeDOFs
1892
     IF(minEdgeDOFs <= maxEdgeDOFs ) THEN
1893
       Mesh % MinEdgeDOFs = minEdgeDOFs
1894
     ELSE
1895
       Mesh % MinEdgeDOFs = maxEdgeDOFs
1896
     END IF
1897

1898
     Mesh % MaxFaceDOFs = maxFaceDOFs
1899
     IF(minFaceDOFs <= maxFaceDOFs ) THEN
1900
       Mesh % MinFaceDOFs = minFaceDOFs
1901
     ELSE
1902
       Mesh % MinFaceDOFs = maxFaceDOFs
1903
     END IF
1904
     
1905
     Mesh % MaxBDOFs = BDOFs
1906
     
1907
     IF (PRESENT( Equation)) THEN
1908
       n = LEN(Equation)
1909
       ALLOCATE(CHARACTER(n)::Eq)
1910
       n=StringToLowerCase(Eq,Equation)
1911
     ELSE
1912
       Eq = ' '
1913
     END IF
1914

1915
     
1916
     IF( UseGiven ) THEN
1917
       k = MAXVAL( Perm ) 
1918
       ALLOCATE( InvInitialReorder(k), STAT=istat )
1919
       IF( istat /= 0 ) THEN
1920
         CALL Fatal(Caller,'Allocation error for InvInitialReorder of size: '//I2S(k))
1921
       END IF
1922
       InvInitialReorder = 0
1923
       DO i=1,SIZE(Perm)
1924
         IF (Perm(i)>0) InvInitialReorder(Perm(i)) = i
1925
       END DO
1926
       GOTO 10
1927
     END IF
1928

1929
       
1930
     Perm = 0
1931
     IF ( PRESENT(Equation) ) THEN
1932
       CALL Info(Caller,'Creating initial permutation',Level=14)
1933
       k = InitialPermutation( Perm,Model,Solver,Mesh,Eq,DG,GB )
1934
       IF ( k <= 0 ) THEN
1935
         IF(ALLOCATED(InvInitialReorder)) DEALLOCATE(InvInitialReorder)
1936
         RETURN
1937
       END IF
1938
     ELSE
1939
       k = SIZE( Perm )
1940
     END IF
1941
     
1942
     IF ( k == SIZE(Perm) ) THEN
1943
       IF(PRESENT(NodalDofsOnly)) THEN
1944
         IF(NodalDofsOnly) k=Mesh % NumberOfNodes
1945
       END IF       
1946
       DO i=1,k 
1947
         Perm(i) = i
1948
       END DO
1949
     END IF
1950

1951
     IF( ParEnv % PEs > 1 .AND. &
1952
         ListGetLogical( Solver % Values,'Skip Pure Halo Nodes',Found ) ) THEN
1953
       CALL Info(Caller,'Skipping pure halo nodes',Level=14)
1954
       j = 0
1955
       DO i=1,Mesh % NumberOfNodes 
1956
         ! These are pure halo nodes that need not be communicated. They are created only 
1957
         ! for sufficient geometric information on the boundaries.
1958
         IF( .NOT. ANY( ParEnv % Mype == Mesh % ParallelInfo % NeighbourList(i) % Neighbours ) ) THEN
1959
           Perm(i) = 0
1960
         ELSE IF( Perm(i) > 0 ) THEN
1961
           j = j + 1
1962
           Perm(i) = j
1963
         END IF
1964
       END DO
1965
       PRINT *,'Eliminating '//I2S(k-j)//' halo nodes out of '&
1966
           //I2S(k)//' in partition '//I2S(ParEnv % MyPe)
1967
       k = j
1968
     END IF
1969

1970
     
1971
     IF( OptimizeBW ) THEN
1972
       CALL Info(Caller,'Creating inverse of initial order of size: '//I2S(k),Level=14)
1973
       ALLOCATE( InvInitialReorder(k), STAT=istat )
1974
       IF( istat /= 0 ) THEN
1975
         CALL Fatal(Caller,'Allocation error for InvInitialReorder of size: '//I2S(k))
1976
       END IF
1977

1978
       ! We need to keep the initial numbering only in case we optimize the bandwidth!
1979
       InvInitialReorder = 0
1980
       DO i=1,SIZE(Perm)
1981
         IF (Perm(i)>0) InvInitialReorder(Perm(i)) = i
1982
       END DO
1983
     END IF
1984
     
1985
     UseOptimized = ListGetLogical( Solver % Values, &
1986
         'Optimize Bandwidth Use Always', GotIt )
1987
          
1988
10   Matrix => NULL()
1989

1990
     ! check if matrix structures really need to be created:
1991
     ! -----------------------------------------------------
1992
     IF ( ListGetLogical( Solver % Values, 'No matrix',GotIt)) THEN
1993
       IF(ALLOCATED(InvInitialReorder)) DEALLOCATE(InvInitialReorder)
1994
       RETURN
1995
     END IF
1996

1997
     !------------------------------------------------------------------------------
1998
     IF (UseThreads) THEN
1999
       CALL Info(Caller,'Creating threaded list matrix array for equation',Level=14)
2000
       IF ( PRESENT(Equation) ) THEN
2001
         CALL MakeListMatrixArray( Model, Solver, Mesh, ListMatrixArray, Perm, k, Eq, DG, GB,&
2002
               NodalDofsOnly, ProjectorDofs, CalcNonZeros = .FALSE. )
2003
         n = OptimizeBandwidth( ListMatrixArray % Rows, Perm, InvInitialReorder, &
2004
               k, OptimizeBW, UseOptimized, Eq )
2005
       ELSE
2006
         CALL MakeListMatrixArray( Model, Solver, Mesh, ListMatrixArray, Perm, k, &
2007
               DGSolver=DG, GlobalBubbles=GB, NodalDofsOnly=NodalDofsOnly, &
2008
               ProjectorDofs=ProjectorDofs, CalcNonZeros = .FALSE. )
2009
         n = OptimizeBandwidth( ListMatrixArray % Rows, Perm, InvInitialReorder, &
2010
               k, OptimizeBW,UseOptimized, ' ' )
2011
       END IF
2012
       
2013
       !------------------------------------------------------------------------------
2014
       ! Initialize the matrix. Multithreading only supports CRS.
2015
       !------------------------------------------------------------------------------
2016
       CALL Info(Caller,'Initializing list matrix array for equation',Level=14)
2017
       IF ( MatrixFormat == MATRIX_CRS) THEN
2018
         Matrix => CRS_CreateMatrix( DOFs*k, Model % TotalMatrixElements, Ndeg=DOFs, &
2019
             Reorder=Perm, AllocValues=.TRUE., SetRows = .FALSE.)
2020
         Matrix % FORMAT = MatrixFormat
2021
         IF( OptimizeBW ) THEN
2022
           CALL InitializeMatrix( Matrix, k, ListMatrixArray % Rows, &
2023
               DOFs, Perm, InvInitialReorder )
2024
         ELSE
2025
           CALL InitializeMatrix( Matrix, k, ListMatrixArray % Rows, DOFs )          
2026
         END IF
2027
       ELSE
2028
         CALL Fatal(Caller,'Multithreaded startup only supports CRS matrix format')
2029
       END IF
2030
       
2031
       CALL Info(Caller,'Sparse atrix created',Level=14)
2032

2033
       CALL ListMatrixArray_Free( ListMatrixArray )       
2034
     ELSE
2035
       NULLIFY( ListMatrix )
2036
       CALL Info(Caller,'Creating list matrix for equation: '//TRIM(Eq),Level=14)
2037

2038
       IF ( PRESENT(Equation) ) THEN
2039
         CALL MakeListMatrix( Model, Solver, Mesh, ListMatrix, Perm, k, Eq, DG, GB,&
2040
               NodalDofsOnly, ProjectorDofs, CalcNonZeros = .FALSE.)
2041
         n = OptimizeBandwidth( ListMatrix, Perm, InvInitialReorder, &
2042
               k, OptimizeBW, UseOptimized, Eq )
2043
       ELSE
2044
         CALL MakeListMatrix( Model, Solver, Mesh, ListMatrix, Perm, k, &
2045
               DGSolver=DG, GlobalBubbles=GB, NodalDofsOnly=NodalDofsOnly, &
2046
               ProjectorDofs=ProjectorDofs, CalcNonZeros = .FALSE.)
2047
         n = OptimizeBandwidth( ListMatrix, Perm, InvInitialReorder, &
2048
               k, OptimizeBW,UseOptimized, ' ' )
2049
       END IF
2050
       
2051
       !------------------------------------------------------------------------------
2052
       ! Initialize the matrix. 
2053
       !------------------------------------------------------------------------------
2054
       CALL Info(Caller,'Initializing list matrix for equation',Level=14)
2055
       SELECT CASE( MatrixFormat )
2056
       CASE( MATRIX_CRS )
2057
         
2058
         Matrix => CRS_CreateMatrix( DOFs*k, Model % TotalMatrixElements, Ndeg=DOFs, &
2059
             Reorder=Perm, AllocValues=.TRUE., SetRows = .FALSE.)
2060
         Matrix % FORMAT = MatrixFormat
2061
         IF( OptimizeBW ) THEN
2062
           CALL InitializeMatrix( Matrix, k, ListMatrix, &
2063
               DOFs, Perm, InvInitialReorder )
2064
         ELSE
2065
           CALL InitializeMatrix( Matrix, k, ListMatrix, DOFs )
2066
         END IF
2067
      CASE( MATRIX_BAND )
2068
        Matrix => Band_CreateMatrix( DOFs*k, DOFs*n,.FALSE.,.TRUE. )
2069
         
2070
       CASE( MATRIX_SBAND )
2071
         Matrix => Band_CreateMatrix( DOFs*k, DOFs*n,.TRUE.,.TRUE. )
2072
       END SELECT
2073
       CALL Info(Caller,'Sparse matrix created',Level=14)
2074

2075
       CALL List_FreeMatrix( k, ListMatrix )
2076
     END IF
2077
     
2078
     NULLIFY( Matrix % MassValues, Matrix % DampValues, Matrix % Force, Matrix % RHS_im )
2079

2080
     IF (ListGetLogical(Solver % Values, 'Allocate Preconditioning Matrix', GotIt)) THEN
2081
       CALL Info(Caller, 'Allocating for separate preconditioning matrix!', Level=20)
2082
       ALLOCATE(Matrix % PrecValues(SIZE(Matrix % Values)) )
2083
       Matrix % PrecValues = 0.0d0
2084
     END IF
2085

2086
     !------------------------------------------------------------------------------
2087
     Matrix % Solver => Solver
2088
     Matrix % DGMatrix = DG
2089
     Matrix % Subband = DOFs * n
2090
     Matrix % COMPLEX = .FALSE.
2091
     Matrix % FORMAT  = MatrixFormat
2092
!------------------------------------------------------------------------------
2093

2094
     n = ListGetInteger( Solver % Values, 'Constraint DOFs', Found )
2095
     IF ( n>0 ) THEN       
2096
       Matrix % ConstraintMatrix => AllocateMatrix()
2097
       A => Matrix % ConstraintMatrix
2098
       A % NumberOfRows = n
2099
       ALLOCATE( A % Rows(n+1), A % Diag(n), A % RHS(n), &
2100
           ConstrainedNode(Mesh % NumberOfNodes), STAT=istat )
2101
       IF( istat /= 0 ) THEN
2102
         CALL Fatal(Caller,'Allocation error for CRS matrix topology: '//I2S(n))
2103
       END IF
2104

2105
       DO i=1,n
2106
         A % RHS(i:i) = ListGetConstReal( Solver % Values,  &
2107
           'Constraint DOF ' // i2s(i) // ' Value' )
2108
       END DO
2109

2110
       Cols = 0
2111
       A % Rows(1) = 1
2112
       DO i=1,n
2113
         str = 'Constraint DOF '//i2s(i)//' Body' 
2114
         ivals => ListGetIntegerArray( Solver % Values, str, Found )
2115
         IF ( ASSOCIATED(ivals) ) THEN
2116
           ConstrainedNode = .FALSE.
2117
           DO k=1,Solver % Mesh % NumberOfBulkElements
2118
             Element => Solver % Mesh % Elements(k)
2119
             IF ( ALL(ivals /= Element % Bodyid) ) CYCLE
2120
             IF ( ALL( Perm(Element % NodeIndexes) > 0 ) ) &
2121
               ConstrainedNode(Element % NodeIndexes) = .TRUE.
2122
           END DO
2123
           Cols = Cols+DOFs*COUNT(ConstrainedNode)
2124
         END IF
2125

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

2142
       ALLOCATE( A % Cols(cols), A % Values(cols), STAT=istat )
2143
       IF( istat /= 0 ) THEN
2144
         CALL Fatal(Caller,'Allocation error for CRS cols and values: '//I2S(cols))
2145
       END IF
2146
       A % Cols = 0
2147
       A % Values = 0
2148

2149
       DO i=1,n
2150
         str = 'Constraint DOF ' // i2s(i) // ' Body' 
2151
         ivals => ListGetIntegerArray( Solver % Values, str, Found )
2152
         IF ( ASSOCIATED(ivals) ) THEN
2153
           DO k=1,Solver % Mesh % NumberOfBulkElements
2154
             Element => Solver % Mesh % Elements(k)
2155
             IF ( ALL(ivals /= Element % Bodyid) ) CYCLE
2156

2157
             IF ( ALL( Perm(Element % NodeIndexes) > 0 ) ) THEN
2158
               DO p=1,Element % TYPE % NumberOfNodes
2159
                 l = Perm(Element % NodeIndexes(p))
2160
                 DO m=1,DOFs 
2161
                   k1 = DOFs*(l-1)+m
2162
                   CALL CRS_MakeMatrixIndex( A,i,k1 )
2163
                 END DO
2164
               END DO
2165
             END IF
2166
           END DO
2167
         END IF
2168

2169
         str = 'Constraint DOF ' // i2s(i) // ' BC'
2170
         ivals => ListGetIntegerArray( Solver % Values, str, Found )
2171
         IF ( ASSOCIATED(ivals) ) THEN
2172
           DO k=Solver % Mesh % NumberOfBulkElements+1, &
2173
                Solver % Mesh % NumberOfBulkElements+Solver % Mesh % NumberOfBoundaryElements
2174
             Element => Solver % Mesh % Elements(k)
2175
             IF ( ALL(ivals /= Element % Boundaryinfo % Constraint) ) CYCLE
2176
             IF ( ALL(Perm(Element % NodeIndexes) > 0) ) THEN
2177
               DO p=1,Element % TYPE % NumberOfNodes
2178
                 l = Perm(Element % NodeIndexes(p))
2179
                 DO m=1,DOFs 
2180
                   k1 = DOFs*(l-1)+m
2181
                   CALL CRS_MakeMatrixIndex( A,i,k1 )
2182
                 END DO
2183
               END DO
2184
             END IF
2185
           END DO
2186
         END IF
2187
       END DO
2188
       CALL CRS_SortMatrix(A)
2189
     END IF
2190

2191
!     DEALLOCATE( Model % RowNonZeros )
2192
     IF( ALLOCATED(InvInitialReorder) ) DEALLOCATE( InvInitialReorder )
2193
!------------------------------------------------------------------------------
2194
   END FUNCTION CreateMatrix
2195
!------------------------------------------------------------------------------
2196

2197

2198
!------------------------------------------------------------------------------
2199
   FUNCTION CreateOdeMatrix( Model, Solver, Dofs ) RESULT(Matrix)
2200
!------------------------------------------------------------------------------
2201
     TYPE(Model_t) :: Model
2202
     TYPE(Solver_t), TARGET :: Solver
2203
     INTEGER :: DOFs
2204
     TYPE(Matrix_t), POINTER :: Matrix
2205
!------------------------------------------------------------------------------
2206
     LOGICAL :: Found
2207
     INTEGER i,j,k
2208
!------------------------------------------------------------------------------
2209

2210
     Matrix => NULL()
2211

2212
     IF ( ListGetLogical( Solver % Values, 'No matrix',Found)) RETURN
2213
     
2214
     ! Create a list matrix that allows for unspecified entries in the matrix 
2215
     ! structure to be introduced.
2216
     Matrix => AllocateMatrix()
2217
     Matrix % FORMAT = MATRIX_LIST
2218
     
2219
     ! Initialize matrix indices
2220
     DO i = 1, Dofs
2221
       DO j = 1, Dofs
2222
          CALL List_AddMatrixIndex(Matrix % ListMatrix, i, j) 
2223
       END DO
2224
     END DO
2225

2226
     CALL List_ToCRSMatrix(Matrix)
2227
     CALL CRS_SortMatrix(Matrix,.TRUE.)
2228
     
2229
     CALL Info('CreateOdeMatrix','Number of rows in ode matrix: '//&
2230
         I2S(Matrix % NumberOfRows), Level=9)
2231
     CALL Info('CreateOdeMatrix','Number of entries in ode matrix: '//&
2232
         I2S(SIZE(Matrix % Cols) ), Level=9)
2233
     
2234
     Matrix % Solver => Solver
2235
     Matrix % DGMatrix = .FALSE.
2236
     Matrix % Subband = DOFs
2237
     Matrix % COMPLEX = .FALSE.
2238
     ! Matrix % FORMAT  = MatrixFormat
2239

2240
!------------------------------------------------------------------------------
2241
   END FUNCTION CreateOdeMatrix
2242
!------------------------------------------------------------------------------
2243

2244

2245
!------------------------------------------------------------------------------
2246
   FUNCTION CreateDiagMatrix( Model, Solver, Dofs, TimeOrder ) RESULT(Matrix)
2247
!------------------------------------------------------------------------------
2248
     TYPE(Model_t) :: Model
2249
     TYPE(Solver_t), TARGET :: Solver
2250
     INTEGER :: DOFs
2251
     TYPE(Matrix_t), POINTER :: Matrix
2252
     INTEGER, OPTIONAL :: TimeOrder
2253
!------------------------------------------------------------------------------
2254
     LOGICAL :: Found
2255
     INTEGER i,j,k
2256
!------------------------------------------------------------------------------
2257

2258
     Matrix => NULL()
2259

2260
     !IF ( ListGetLogical( Solver % Values, 'No matrix',Found)) RETURN
2261
     
2262
     ! Create a list matrix that allows for unspecified entries in the matrix 
2263
     ! structure to be introduced.
2264
     Matrix => AllocateMatrix()
2265
     Matrix % FORMAT = MATRIX_LIST
2266
     
2267
     ! Initialize matrix indices
2268
     DO i = 1, Dofs
2269
       CALL List_AddMatrixIndex(Matrix % ListMatrix, i, i) 
2270
     END DO
2271

2272
     CALL List_ToCRSMatrix(Matrix)
2273
     CALL CRS_SortMatrix(Matrix,.TRUE.)
2274
     
2275
     CALL Info('CreateOdeMatrix','Number of rows in diag matrix: '//&
2276
         I2S(Matrix % NumberOfRows), Level=9)
2277

2278
     IF( PRESENT( TimeOrder ) ) THEN
2279
       IF( TimeOrder >= 1 ) THEN
2280
         ALLOCATE( Matrix % MassValues( SIZE( Matrix % Values ) ) )
2281
         Matrix % MassValues = 0.0_dp
2282
       END IF
2283
       IF( TimeOrder >= 2 ) THEN
2284
         ALLOCATE( Matrix % DampValues( SIZE( Matrix % Values ) ) )
2285
         Matrix % DampValues = 0.0_dp
2286
       END IF
2287
     END IF
2288
     
2289
     Matrix % Solver => Solver
2290
     Matrix % DGMatrix = .FALSE.
2291
     Matrix % Subband = 1
2292
     Matrix % COMPLEX = .FALSE.
2293
     ! Matrix % FORMAT  = MatrixFormat
2294

2295
!------------------------------------------------------------------------------
2296
   END FUNCTION CreateDiagMatrix
2297
!------------------------------------------------------------------------------
2298

2299

2300

2301
!------------------------------------------------------------------------------
2302
  SUBROUTINE RotateMatrix( Matrix,Vector,n,DIM,DOFs,NodeIndexes,  &
2303
                   Normals,Tangent1,Tangent2 )
2304
!------------------------------------------------------------------------------
2305

2306
    REAL(KIND=dp) :: Matrix(:,:),Vector(:)
2307
    REAL(KIND=dp), POINTER :: Normals(:,:), Tangent1(:,:),Tangent2(:,:)
2308
    INTEGER :: n,DIM,DOFs,NodeIndexes(:)
2309
!------------------------------------------------------------------------------
2310

2311
    INTEGER :: i,j,k,l
2312
    REAL(KIND=dp) :: s,R(n*DOFs,n*DOFs),Q(n*DOFs,n*DOFs),N1(3),T1(3),T2(3)
2313
!------------------------------------------------------------------------------
2314

2315
    DO i=1,MIN(n,SIZE(NodeIndexes))
2316
      IF ( NodeIndexes(i)<=0 .OR. NodeIndexes(i)>SIZE(Normals,1) ) CYCLE
2317

2318
      R = 0.0d0
2319
      DO j=1,n*DOFs
2320
        R(j,j) = 1.0d0
2321
      END DO
2322

2323
      N1 = Normals( NodeIndexes(i),: )
2324

2325
      SELECT CASE(DIM)
2326
      CASE (2)
2327
        R(DOFs*(i-1)+1,DOFs*(i-1)+1) =  N1(1)
2328
        R(DOFs*(i-1)+1,DOFs*(i-1)+2) =  N1(2)
2329

2330
        R(DOFs*(i-1)+2,DOFs*(i-1)+1) = -N1(2)
2331
        R(DOFs*(i-1)+2,DOFs*(i-1)+2) =  N1(1)
2332
      CASE (3)
2333
        T1 = Tangent1( NodeIndexes(i),: )
2334
        T2 = Tangent2( NodeIndexes(i),: )
2335

2336
        R(DOFs*(i-1)+1,DOFs*(i-1)+1) = N1(1)
2337
        R(DOFs*(i-1)+1,DOFs*(i-1)+2) = N1(2)
2338
        R(DOFs*(i-1)+1,DOFs*(i-1)+3) = N1(3)
2339

2340
        R(DOFs*(i-1)+2,DOFs*(i-1)+1) = T1(1)
2341
        R(DOFs*(i-1)+2,DOFs*(i-1)+2) = T1(2)
2342
        R(DOFs*(i-1)+2,DOFs*(i-1)+3) = T1(3)
2343

2344
        R(DOFs*(i-1)+3,DOFs*(i-1)+1) = T2(1)
2345
        R(DOFs*(i-1)+3,DOFs*(i-1)+2) = T2(2)
2346
        R(DOFs*(i-1)+3,DOFs*(i-1)+3) = T2(3)
2347
      END SELECT
2348

2349
      DO j=1,n*DOFs
2350
        DO k=1,n*DOFs
2351
          s = 0.0D0
2352
          DO l=1,n*DOFs
2353
            s = s + R(j,l) * Matrix(l,k)
2354
          END DO
2355
          Q(j,k) = s
2356
        END DO
2357
      END DO
2358

2359
      DO j=1,n*DOFs
2360
        DO k=1,n*DOFs
2361
          s = 0.0D0
2362
          DO l=1,n*DOFs
2363
            s = s + Q(j,l) * R(k,l)
2364
          END DO
2365
          Matrix(j,k) = s
2366
        END DO
2367
      END DO
2368

2369
      DO j=1,n*DOFs
2370
        s = 0.0D0
2371
        DO k=1,n*DOFs
2372
          s = s + R(j,k) * Vector(k)
2373
        END DO
2374
        Q(j,1) = s
2375
      END DO
2376
      Vector(1:n*DOFs) = Q(:,1)
2377
    END DO
2378
!------------------------------------------------------------------------------
2379
  END SUBROUTINE RotateMatrix
2380
!------------------------------------------------------------------------------
2381

2382

2383

2384
!------------------------------------------------------------------------------
2385
!> Given the normal return the tangent directions. The 
2386
!> First tangent direction will always be on the xy-plane if
2387
!> also the normal is in the xy-plane.
2388
!------------------------------------------------------------------------------
2389
  SUBROUTINE TangentDirections( Normal,Tangent1,Tangent2 )
2390
!------------------------------------------------------------------------------
2391
   REAL(KIND=dp) :: Normal(3),Tangent1(3)
2392
   REAL(KIND=dp), OPTIONAL :: Tangent2(3)
2393
!------------------------------------------------------------------------------
2394
   REAL(KIND=dp) :: n1,n2,n3
2395
!------------------------------------------------------------------------------
2396
   n1 = ABS(Normal(1))
2397
   n2 = ABS(Normal(2))
2398
   n3 = ABS(Normal(3))
2399

2400
   IF( PRESENT( Tangent2 ) ) THEN   
2401
     IF ( n1 <= n3 .AND. n2 <= n3 ) THEN
2402
       Tangent1(1) =  0.0_dp
2403
       Tangent1(2) = -Normal(3)
2404
       Tangent1(3) =  Normal(2)
2405
     ELSE
2406
       Tangent1(1) = -Normal(2)
2407
       Tangent1(2) =  Normal(1)
2408
       Tangent1(3) =  0.0_dp
2409
     END IF
2410

2411
     Tangent1 = Tangent1 / SQRT(SUM(Tangent1**2))
2412
     Tangent2(1) = Normal(2)*Tangent1(3) - Normal(3)*Tangent1(2)
2413
     Tangent2(2) = Normal(3)*Tangent1(1) - Normal(1)*Tangent1(3)
2414
     Tangent2(3) = Normal(1)*Tangent1(2) - Normal(2)*Tangent1(1)
2415
     Tangent2 = Tangent2 / SQRT(SUM(Tangent2**2))
2416
   ELSE
2417
     ! This is a 2D tangent only
2418
     Tangent1(1) = Normal(2)
2419
     Tangent1(2) = -Normal(1)
2420
     Tangent1(3) = 0.0_dp
2421
     Tangent1 = Tangent1 / SQRT(SUM(Tangent1**2))
2422
   END IF
2423
     
2424
!------------------------------------------------------------------------------
2425
 END SUBROUTINE TangentDirections
2426
!------------------------------------------------------------------------------
2427

2428

2429
!------------------------------------------------------------------------------
2430
!> Given one or two tangents return the normal direction.
2431
!> It is assumed that if one tangent is given the normal is in xy-plane otherwise
2432
!> in 3D. Note that the sign of normal vector is not unique.
2433
!------------------------------------------------------------------------------
2434
  FUNCTION NormalDirection( Tangent1,Tangent2 ) RESULT ( Normal ) 
2435
!------------------------------------------------------------------------------
2436
   REAL(KIND=dp) :: Normal(3),Tangent1(3)
2437
   REAL(KIND=dp), OPTIONAL :: Tangent2(3)
2438
!------------------------------------------------------------------------------
2439
   REAL(KIND=dp) :: t1(3),t2(3)
2440
!------------------------------------------------------------------------------
2441

2442
   IF(PRESENT(Tangent2)) THEN
2443
     Normal = CrossProduct(Tangent1,Tangent2)
2444
   ELSE
2445
     Normal(1) = Tangent1(2)
2446
     Normal(2) = -Tangent1(1)
2447
   END IF
2448
   Normal = Normal/SQRT(SUM(Normal**2))       
2449
!------------------------------------------------------------------------------
2450
 END FUNCTION NormalDirection
2451
!------------------------------------------------------------------------------
2452

2453

2454
!------------------------------------------------------------------------------
2455
!>    Integrates a user-defined function over the specified bulk elements.
2456
!> \deprecated Is this used for something?
2457
!------------------------------------------------------------------------------
2458
   FUNCTION VolumeIntegrate( Model, ElementList, IntegrandFunctionName ) &
2459
       RESULT(Integral)
2460
!------------------------------------------------------------------------------
2461
   TYPE(Model_t) :: Model                     !< All model information (mesh, materials, BCs, etc...)
2462
   INTEGER, DIMENSION(:) :: ElementList       !< List of elements that belong to the integration volume
2463
   CHARACTER(LEN=*) :: IntegrandFunctionName  !< Name the function has in the .sif file
2464
   REAL(KIND=dp) :: Integral                  !< The value of the volume integral
2465
!------------------------------------------------------------------------------
2466
     INTEGER :: n
2467
     TYPE(Element_t), POINTER :: CurrentElement
2468
     INTEGER, POINTER :: NodeIndexes(:)
2469
     TYPE(Nodes_t)   :: ElementNodes
2470

2471
     TYPE(GaussIntegrationPoints_t), TARGET :: IntegStuff
2472
     REAL(KIND=dp), DIMENSION(:), POINTER :: &
2473
         U_Integ,V_Integ,W_Integ,S_Integ
2474

2475
     REAL(KIND=dp), DIMENSION(Model % MaxElementNodes) :: IntegrandFunction
2476
     REAL(KIND=dp) :: s,ug,vg,wg
2477
     REAL(KIND=dp) :: Basis(Model % MaxElementNodes)
2478
     REAL(KIND=dp) :: dBasisdx(Model % MaxElementNodes,3),SqrtElementMetric
2479
     REAL(KIND=dp) :: IntegrandAtGPt, dV
2480
     INTEGER :: N_Integ, t, tg, i, istat
2481
     LOGICAL :: stat
2482

2483
! Need MaxElementNodes only in allocation
2484
     n = Model % MaxElementNodes
2485
     ALLOCATE( ElementNodes % x( n ),   &
2486
               ElementNodes % y( n ),   &
2487
               ElementNodes % z( n ), STAT=istat )
2488
     IF( istat /= 0 ) THEN
2489
       CALL Fatal('VolumeIntegrate','Allocation error for ElementNodes')
2490
     END IF
2491

2492
     Integral = 0.0d0
2493

2494
! Loop over all elements in the list
2495
     DO i=1,SIZE(ElementList)
2496

2497
       t = ElementList(i)
2498

2499
       IF ( t < 1 .OR. t > Model % NumberOfBulkElements ) THEN
2500
! do something
2501
       END IF
2502

2503
       CurrentElement => Model % Elements(t)
2504
       n = CurrentElement % TYPE % NumberOfNodes
2505
       NodeIndexes => CurrentElement % NodeIndexes
2506

2507
!------------------------------------------------------------------------------
2508
! Get element nodal coordinates
2509
!------------------------------------------------------------------------------
2510
       ElementNodes % x(1:n) = Model % Nodes % x(NodeIndexes)
2511
       ElementNodes % y(1:n) = Model % Nodes % y(NodeIndexes)
2512
       ElementNodes % z(1:n) = Model % Nodes % z(NodeIndexes)
2513

2514
! Read this from Simulation block in the .sif file
2515
       IntegrandFunction(1:n) = ListGetReal( Model % Simulation, &
2516
           IntegrandFunctionName, n, NodeIndexes )
2517

2518
!------------------------------------------------------------------------------
2519
!    Gauss integration stuff
2520
!------------------------------------------------------------------------------
2521
       IntegStuff = GaussPoints( CurrentElement )
2522
       U_Integ => IntegStuff % u
2523
       V_Integ => IntegStuff % v
2524
       W_Integ => IntegStuff % w
2525
       S_Integ => IntegStuff % s
2526
       N_Integ =  IntegStuff % n
2527

2528
!------------------------------------------------------------------------------
2529
! Loop over Gauss integration points
2530
!------------------------------------------------------------------------------
2531
       DO tg=1,N_Integ
2532

2533
         ug = U_Integ(tg)
2534
         vg = V_Integ(tg)
2535
         wg = W_Integ(tg)
2536

2537
!------------------------------------------------------------------------------
2538
! Need SqrtElementMetric and Basis at the integration point
2539
!------------------------------------------------------------------------------
2540
         stat = ElementInfo( CurrentElement,ElementNodes,ug,vg,wg, &
2541
             SqrtElementMetric,Basis,dBasisdx )
2542

2543
         s = SqrtElementMetric * S_Integ(tg)
2544

2545
! Calculate the function to be integrated at the Gauss point
2546
         IntegrandAtGPt = SUM( IntegrandFunction(1:n) * Basis )
2547

2548
! Use general coordinate system for dV
2549
         dV = CoordinateSqrtMetric( SUM( ElementNodes % x(1:n) * Basis), &
2550
             SUM( ElementNodes % y(1:n) * Basis), &
2551
             SUM( ElementNodes % z(1:n) * Basis) )
2552

2553
         Integral = Integral + s*IntegrandAtGPt*dV
2554

2555
       END DO! of the Gauss integration points
2556

2557
     END DO! of the bulk elements
2558

2559
     DEALLOCATE( ElementNodes % x, &
2560
         ElementNodes % y, &
2561
         ElementNodes % z )
2562
!------------------------------------------------------------------------------
2563
   END FUNCTION VolumeIntegrate
2564
!------------------------------------------------------------------------------
2565

2566

2567

2568
!------------------------------------------------------------------------------
2569
!>    Integrates the normal component of a user-defined vector function
2570
!>    over the specified boundary elements. 
2571
!>  \deprecated Is this used for something?
2572
!------------------------------------------------------------------------------
2573
   FUNCTION FluxIntegrate( Model, ElementList, IntegrandFunctionName ) &
2574
!------------------------------------------------------------------------------
2575
       RESULT(Integral)
2576
!------------------------------------------------------------------------------
2577
   TYPE(Model_t) :: Model                     !< All model information (mesh, materials, BCs, etc...)
2578
   INTEGER, DIMENSION(:) :: ElementList       !< List of elements that belong to the integration boundary
2579
   CHARACTER(LEN=*) :: IntegrandFunctionName  !< Name the function has in the .sif file
2580
   REAL(KIND=dp) :: Integral                  !< The value of the boundary integral
2581
!------------------------------------------------------------------------------
2582
     INTEGER :: n
2583
     TYPE(Element_t), POINTER :: CurrentElement
2584
     INTEGER, POINTER :: NodeIndexes(:)
2585
     TYPE(Nodes_t)    :: ElementNodes
2586

2587
     TYPE(GaussIntegrationPoints_t), TARGET :: IntegStuff
2588
     REAL(KIND=dp), DIMENSION(:), POINTER :: &
2589
         U_Integ,V_Integ,W_Integ,S_Integ
2590

2591
     REAL(KIND=dp), DIMENSION(Model % MaxElementNodes,3) :: IntegrandFunction
2592
!     REAL(KIND=dp), POINTER :: IntegrandFunction(:,:)
2593
     CHARACTER(LEN=2) :: Component
2594
     CHARACTER(:), ALLOCATABLE :: IntegrandFunctionComponent
2595
     REAL(KIND=dp) :: s,ug,vg,wg
2596
     REAL(KIND=dp) :: Basis(Model % MaxElementNodes)
2597
     REAL(KIND=dp) :: dBasisdx(Model % MaxElementNodes,3),SqrtElementMetric
2598
     REAL(KIND=dp) :: Normal(3)
2599
     REAL(KIND=dp) :: IntegrandAtGPt(3), FluxAtGPt, dS
2600
     INTEGER :: N_Integ, t, tg, i, j, DIM, istat
2601
     LOGICAL :: stat
2602

2603
     DIM = CoordinateSystemDimension()
2604

2605
! Need MaxElementNodes only in allocation
2606
     n = Model % MaxElementNodes
2607
     ALLOCATE( ElementNodes % x( n ),   &
2608
               ElementNodes % y( n ),   &
2609
               ElementNodes % z( n ), STAT=istat )
2610
     IF( istat /= 0 ) THEN
2611
       CALL Fatal('FluxIntegrate','Allocation error for ElementNodes')
2612
     END IF
2613

2614
     Integral = 0.0d0
2615

2616
! Loop over all elements in the list
2617
     DO i=1,SIZE(ElementList)
2618

2619
       t = ElementList(i)
2620

2621
       IF ( t < 1 .OR. t > Model % NumberOfBulkElements ) THEN
2622
! do something
2623
       END IF
2624

2625
       CurrentElement => Model % Elements(t)
2626
       n = CurrentElement % TYPE % NumberOfNodes
2627
       NodeIndexes => CurrentElement % NodeIndexes
2628

2629
       Model % CurrentElement => Model % Elements(t)
2630
!------------------------------------------------------------------------------
2631
! Get element nodal coordinates
2632
!------------------------------------------------------------------------------
2633
       ElementNodes % x(1:n) = Model % Nodes % x(NodeIndexes)
2634
       ElementNodes % y(1:n) = Model % Nodes % y(NodeIndexes)
2635
       ElementNodes % z(1:n) = Model % Nodes % z(NodeIndexes)
2636

2637
! Read the integrand from Simulation block in the .sif file
2638
! It is assumed to be a contravariant vector, but
2639
! ListGetRealArray doesn t exist, so we READ it component by component
2640
! naming them with suffixes " 1" etc.
2641
       DO j=1,DIM
2642
         IntegrandFunctionComponent = TRIM(IntegrandFunctionName)//' '//I2S(j)
2643
         IntegrandFunction(1:n,j) = ListGetReal( Model % Simulation, &
2644
             IntegrandFunctionComponent, n, NodeIndexes )
2645
       END DO
2646

2647
!------------------------------------------------------------------------------
2648
!    Gauss integration stuff
2649
!------------------------------------------------------------------------------
2650
       IntegStuff = GaussPoints( CurrentElement )
2651
       U_Integ => IntegStuff % u
2652
       V_Integ => IntegStuff % v
2653
       W_Integ => IntegStuff % w
2654
       S_Integ => IntegStuff % s
2655
       N_Integ =  IntegStuff % n
2656

2657
!------------------------------------------------------------------------------
2658
! Loop over Gauss integration points
2659
!------------------------------------------------------------------------------
2660
       DO tg=1,N_Integ
2661

2662
         ug = U_Integ(tg)
2663
         vg = V_Integ(tg)
2664
         wg = W_Integ(tg)
2665

2666
!------------------------------------------------------------------------------
2667
! Need SqrtElementMetric and Basis at the integration point
2668
!------------------------------------------------------------------------------
2669
         stat = ElementInfo( CurrentElement,ElementNodes,ug,vg,wg, &
2670
             SqrtElementMetric,Basis,dBasisdx )
2671

2672
! If we want to allow covariant integrand vectors given
2673
! we would also need Metric, and read it as follows...
2674
!     IF ( CurrentCoordinateSystem() /= Cartesian ) THEN
2675
!       x = SUM( nodes % x(1:n) * Basis )
2676
!       y = SUM( nodes % y(1:n) * Basis )
2677
!       z = SUM( nodes % z(1:n) * Basis )
2678
!     END IF
2679
!
2680
!     CALL CoordinateSystemInfo( Metric,SqrtMetric,Symb,dSymb,x,y,z )
2681
!     s = SqrtMetric * SqrtElementMetric * S_Integ(t)
2682
! And also FluxAtGpt...
2683
! Now we won t need Metric, so get SqrtMetric later...
2684

2685
         s = SqrtElementMetric * S_Integ(tg)
2686

2687
! Get normal to the element boundary at the integration point
2688
! N.B. NormalVector returns covariant normal vector
2689
         Normal = NormalVector( CurrentElement,ElementNodes,ug,vg,.TRUE. )
2690

2691
! Calculate the contravariant vector function to be integrated
2692
! at the Gauss point
2693
         DO j=1,DIM
2694
           IntegrandAtGPt(j) = SUM( IntegrandFunction(1:n,j) * Basis )
2695
         END DO
2696
! Calculate the normal component of the vector function
2697
         FluxAtGPt = SUM( IntegrandAtGPt * Normal )
2698

2699
! Use general coordinate system for dS
2700
! Would be included in s by SqrtMetric
2701
         dS = CoordinateSqrtMetric( SUM( ElementNodes % x(1:n) * Basis), &
2702
             SUM( ElementNodes % y(1:n) * Basis), &
2703
             SUM( ElementNodes % z(1:n) * Basis) )
2704

2705
         Integral = Integral + s*FluxAtGPt*dS
2706

2707
       END DO! of the Gauss integration points
2708

2709
     END DO! of the boundary elements
2710

2711
     DEALLOCATE( ElementNodes % x, &
2712
         ElementNodes % y, &
2713
         ElementNodes % z )
2714
!------------------------------------------------------------------------------
2715
   END FUNCTION FluxIntegrate
2716
!------------------------------------------------------------------------------
2717

2718

2719

2720
!------------------------------------------------------------------------------
2721
!>    Integrates A user-defined vector function 
2722
!>    over the specified boundary elements.
2723
!> \deprecated Is this used for something?
2724
!------------------------------------------------------------------------------
2725
   FUNCTION SurfaceIntegrate( Model, ElementList, IntegrandFunctionName ) &
2726
       RESULT(Integral)
2727
!------------------------------------------------------------------------------
2728
   TYPE(Model_t) :: Model                     !< All model information (mesh, materials, BCs, etc...)
2729
   INTEGER, DIMENSION(:) :: ElementList       !< List of elements that belong to the integration boundary
2730
   CHARACTER(LEN=*) :: IntegrandFunctionName  !< Name the function has in the .sif file
2731
   REAL(KIND=dp), DIMENSION(3) :: Integral    !< The vector value of the integral
2732
!------------------------------------------------------------------------------
2733
     INTEGER :: n
2734
     TYPE(Element_t), POINTER :: CurrentElement
2735
     INTEGER, POINTER :: NodeIndexes(:)
2736
     TYPE(Nodes_t)    :: ElementNodes
2737

2738
     TYPE(GaussIntegrationPoints_t), TARGET :: IntegStuff
2739
     REAL(KIND=dp), DIMENSION(:), POINTER :: &
2740
         U_Integ,V_Integ,W_Integ,S_Integ
2741

2742
     REAL(KIND=dp), DIMENSION(Model % MaxElementNodes,3) :: IntegrandFunction
2743
!     REAL(KIND=dp), POINTER :: IntegrandFunction(:,:)
2744
     CHARACTER(LEN=2) :: Component
2745
     CHARACTER(:), ALLOCATABLE :: IntegrandFunctionComponent
2746
     REAL(KIND=dp) :: s,ug,vg,wg
2747
     REAL(KIND=dp) :: Basis(Model % MaxElementNodes)
2748
     REAL(KIND=dp) :: dBasisdx(Model % MaxElementNodes,3),SqrtElementMetric
2749
!     REAL(KIND=dp) :: Normal(3)
2750
     REAL(KIND=dp) :: IntegrandAtGPt(3), dS
2751
     INTEGER :: N_Integ, t, tg, i, j, DIM, istat
2752
     LOGICAL :: stat
2753

2754
     DIM = CoordinateSystemDimension()
2755

2756
! Need MaxElementNodes only in allocation
2757
     n = Model % MaxElementNodes
2758
     ALLOCATE( ElementNodes % x( n ),   &
2759
               ElementNodes % y( n ),   &
2760
               ElementNodes % z( n ), STAT=istat )
2761
     IF( istat /= 0 ) THEN
2762
       CALL Fatal('SurfaceIntegrate','Allocation error for ElementNodes')
2763
     END IF
2764

2765
     Integral = 0.0d0
2766

2767
! Loop over all elements in the list
2768
     DO i=1,SIZE(ElementList)
2769

2770
       t = ElementList(i)
2771

2772
       IF ( t < 1 .OR. t > Model % NumberOfBulkElements ) THEN
2773
! do something
2774
       END IF
2775

2776
       CurrentElement => Model % Elements(t)
2777
       Model % CurrentElement => CurrentElement
2778
       n = CurrentElement % TYPE % NumberOfNodes
2779
       NodeIndexes => CurrentElement % NodeIndexes
2780

2781
!------------------------------------------------------------------------------
2782
! Get element nodal coordinates
2783
!------------------------------------------------------------------------------
2784
       ElementNodes % x(1:n) = Model % Nodes % x(NodeIndexes)
2785
       ElementNodes % y(1:n) = Model % Nodes % y(NodeIndexes)
2786
       ElementNodes % z(1:n) = Model % Nodes % z(NodeIndexes)
2787

2788
! Read the integrand from Simulation block in the .sif file
2789
! It is assumed to be a contravariant vector, but
2790
! ListGetRealArray doesn t exist, so we READ it component by component
2791
! naming them with suffixes " 1" etc.
2792
       DO j=1,DIM
2793
         IntegrandFunctionComponent = TRIM(IntegrandFunctionName)//' '//I2S(j)
2794
         IntegrandFunction(1:n,j) = ListGetReal( Model % Simulation, &
2795
             IntegrandFunctionComponent, n, NodeIndexes )
2796
       END DO
2797

2798
!------------------------------------------------------------------------------
2799
!    Gauss integration stuff
2800
!------------------------------------------------------------------------------
2801
       IntegStuff = GaussPoints( CurrentElement )
2802
       U_Integ => IntegStuff % u
2803
       V_Integ => IntegStuff % v
2804
       W_Integ => IntegStuff % w
2805
       S_Integ => IntegStuff % s
2806
       N_Integ =  IntegStuff % n
2807

2808
!------------------------------------------------------------------------------
2809
! Loop over Gauss integration points
2810
!------------------------------------------------------------------------------
2811
       DO tg=1,N_Integ
2812

2813
         ug = U_Integ(tg)
2814
         vg = V_Integ(tg)
2815
         wg = W_Integ(tg)
2816

2817
!------------------------------------------------------------------------------
2818
! Need SqrtElementMetric and Basis at the integration point
2819
!------------------------------------------------------------------------------
2820
         stat = ElementInfo( CurrentElement,ElementNodes,ug,vg,wg, &
2821
             SqrtElementMetric,Basis,dBasisdx )
2822

2823
! If we want to allow covariant integrand vectors given
2824
! we would also need Metric, and read it as follows...
2825
!     IF ( CurrentCoordinateSystem() /= Cartesian ) THEN
2826
!       x = SUM( nodes % x(1:n) * Basis )
2827
!       y = SUM( nodes % y(1:n) * Basis )
2828
!       z = SUM( nodes % z(1:n) * Basis )
2829
!     END IF
2830
!
2831
!     CALL CoordinateSystemInfo( Metric,SqrtMetric,Symb,dSymb,x,y,z )
2832
!     s = SqrtMetric * SqrtElementMetric * S_Integ(t)
2833
! Now we won t need Metric, so get SqrtMetric later...
2834

2835
         s = SqrtElementMetric * S_Integ(tg)
2836

2837
! If you need normal directly at the integration point
2838
!         Normal = NormalVector( CurrentElement,ElementNodes,ug,vg,.TRUE. )
2839

2840
! Calculate the contravariant vector function to be integrated
2841
! at the Gauss point
2842
         DO j=1,DIM
2843
           IntegrandAtGPt(j) = SUM( IntegrandFunction(1:n,j) * Basis )
2844
         END DO
2845

2846
! Use general coordinate system for dS
2847
! Would be included in s by SqrtMetric
2848
         dS = CoordinateSqrtMetric( SUM( ElementNodes % x(1:n) * Basis), &
2849
             SUM( ElementNodes % y(1:n) * Basis), &
2850
             SUM( ElementNodes % z(1:n) * Basis) )
2851

2852
         DO j=1,DIM
2853
           Integral(j) = Integral(j) + s*IntegrandAtGPt(j)*dS
2854
         END DO
2855

2856
       END DO! of the Gauss integration points
2857

2858
     END DO! of the boundary elements
2859

2860
     DEALLOCATE( ElementNodes % x, &
2861
         ElementNodes % y, &
2862
         ElementNodes % z )
2863
!------------------------------------------------------------------------------
2864
   END FUNCTION SurfaceIntegrate
2865
!------------------------------------------------------------------------------
2866

2867

2868

2869
!------------------------------------------------------------------------------
2870
!>    Integrates the normal component of a user-defined vector function
2871
!>    over a specified line element.
2872
!> \deprecated Is this used for something?
2873
!------------------------------------------------------------------------------
2874
   FUNCTION LineIntegrate( Model, LineElement, LineElementNodes, &
2875
       IntegrandFunctionName, QuadrantTreeExists, RootQuadrant ) &
2876
       RESULT(Integral)
2877
!------------------------------------------------------------------------------
2878
!
2879
!  ARGUMENTS:
2880
!
2881
!  TYPE(Model_t), POINTER :: Model,  
2882
!     INPUT: All model information (mesh, materials, BCs, etc...)
2883
!
2884
!  TYPE(Element_t) :: LineElement
2885
!     INPUT: Line element that belongs to the line of integration
2886
!
2887
!  REAL(KIND=dp), DIMENSION(LineElement % Type % NumberOfNodes,3) ::
2888
!     LineElementNodes
2889
!     INPUT: List of nodal point coordinates
2890
!
2891
!  CHARACTER(LEN=*) :: IntegrandFunctionName
2892
!     INPUT: Name the function has in the .sif file or somewhere else
2893
!
2894
!  LOGICAL :: QuadrantTreeExists
2895
!     INPUT: QuadrantTree has been built, use it in element search
2896
!
2897
!  TYPE(Quadrant_t), POINTER :: RootQuadrant
2898
!     OUTPUT: Quadrant tree structure root
2899
!
2900
!  FUNCTION RETURN VALUE:
2901
!    REAL(KIND=dp) :: Integral
2902
!     The value of the flux integral
2903
!      
2904
!------------------------------------------------------------------------------
2905
   TYPE(Model_t) :: Model
2906
! TARGET only for CurrentElement
2907
   TYPE(Element_t), TARGET :: LineElement
2908
   REAL(KIND=dp), DIMENSION(:,:), TARGET CONTIG :: LineElementNodes
2909
   CHARACTER(LEN=*) :: IntegrandFunctionName
2910
   REAL(KIND=dp) :: Integral
2911
   LOGICAL :: QuadrantTreeExists
2912
   TYPE(Quadrant_t), POINTER :: RootQuadrant
2913
!------------------------------------------------------------------------------
2914
     INTEGER :: n
2915
! Only one element at a time, need CurrentElement only for NormalVector!
2916
     TYPE(Element_t), POINTER :: CurrentElement
2917
! LineElement nodes don t belong to global node structure
2918
! Need the structure for the function calls
2919
     TYPE(Nodes_t) :: ElementNodes
2920

2921
     TYPE(GaussIntegrationPoints_t), TARGET :: IntegStuff
2922
     REAL(KIND=dp), DIMENSION(:), POINTER :: &
2923
         U_Integ,V_Integ,W_Integ,S_Integ
2924

2925
! IntegrandFunction at the bulk element nodal points
2926
     REAL(KIND=dp), DIMENSION(Model % MaxElementNodes,3) :: IntegrandAtNodes
2927
! IntegrandFunction at the Gauss points
2928
     REAL(KIND=dp), DIMENSION(LineElement % TYPE % GaussPoints,3) :: IntegrandFunction
2929
     CHARACTER(LEN=2) :: Component
2930
     CHARACTER(:), ALLOCATABLE :: IntegrandFunctionComponent
2931
     REAL(KIND=dp) :: s,ug,vg,wg
2932
     REAL(KIND=dp) :: Basis(LineElement % TYPE % NumberOfNodes)
2933
     REAL(KIND=dp) :: dBasisdx(LineElement % TYPE % NumberOfNodes,3),SqrtElementMetric
2934
     REAL(KIND=dp) :: Normal(3)
2935
! IntegrandFunction already at GPts
2936
     REAL(KIND=dp) :: FluxAtGPt, dS
2937
     INTEGER :: N_Integ, t, tg, i, j, DIM
2938
     LOGICAL :: stat
2939
! Search for the bulk element each Gauss point belongs to
2940
     TYPE(Element_t), POINTER :: BulkElement
2941
     TYPE(Nodes_t) :: BulkElementNodes
2942
     INTEGER, POINTER :: NodeIndexes(:)
2943
     REAL(KIND=dp), DIMENSION(3) :: LocalCoordsInBulkElem
2944
     REAL(KIND=dp), DIMENSION(3) :: Point
2945
     INTEGER :: nBulk, maxlevel=10, k, Quadrant
2946
     TYPE(Quadrant_t), POINTER :: LeafQuadrant
2947

2948
     DIM = CoordinateSystemDimension()
2949

2950
     n = LineElement % TYPE % NumberOfNodes
2951

2952
     Integral = 0.0d0
2953

2954
!------------------------------------------------------------------------------
2955
! Get element nodal coordinates
2956
!------------------------------------------------------------------------------
2957
! Move from LineElementNodes to Nodes_t structure
2958
     ElementNodes % x => LineElementNodes(1:n,1)
2959
     ElementNodes % y => LineElementNodes(1:n,2)
2960
     ElementNodes % z => LineElementNodes(1:n,3)
2961

2962
!------------------------------------------------------------------------------
2963
!    Gauss integration stuff
2964
!------------------------------------------------------------------------------
2965
     IntegStuff = GaussPoints( LineElement )
2966
     U_Integ => IntegStuff % u
2967
     V_Integ => IntegStuff % v
2968
     W_Integ => IntegStuff % w
2969
     S_Integ => IntegStuff % s
2970
     N_Integ =  IntegStuff % n
2971

2972
!------------------------------------------------------------------------------
2973
! Loop over Gauss integration points
2974
!------------------------------------------------------------------------------
2975
     DO tg=1,N_Integ
2976

2977
       ug = U_Integ(tg)
2978
       vg = V_Integ(tg)
2979
       wg = W_Integ(tg)
2980

2981
!------------------------------------------------------------------------------
2982
! Need SqrtElementMetric and Basis at the integration point
2983
!------------------------------------------------------------------------------
2984
       stat = ElementInfo( LineElement,ElementNodes,ug,vg,wg, &
2985
           SqrtElementMetric,Basis,dBasisdx )
2986

2987
! If we want to allow covariant integrand vectors given... (see FluxIntegrate)
2988

2989
       s = SqrtElementMetric * S_Integ(tg)
2990

2991
! Find in which bulk element the Gauss point belongs to
2992
       IF ( QuadrantTreeExists ) THEN
2993
! Find the last existing quadrant that the point belongs to
2994
         Point = [ SUM( ElementNodes % x(1:n) * Basis), &
2995
             SUM( ElementNodes % y(1:n) * Basis), &
2996
             SUM( ElementNodes % z(1:n) * Basis) ]
2997
!         PRINT*,'Point:', Point
2998
         CALL FindLeafElements(Point, DIM, RootQuadrant, LeafQuadrant)         
2999
!         PRINT*,'Elems in LeafQuadrant',LeafQuadrant % NElemsInQuadrant
3000
! Go through the bulk elements in the last ChildQuadrant only
3001
         nBulk = Model % MaxElementNodes
3002
         ALLOCATE( BulkElementNodes % x( nBulk ),   &
3003
             BulkElementNodes % y( nBulk ),   &
3004
             BulkElementNodes % z( nBulk ) )
3005
!         PRINT*,'Elements:', LeafQuadrant % Elements
3006
         DO k=1, LeafQuadrant % NElemsInQuadrant
3007
           BulkElement => Model % Elements( &
3008
               LeafQuadrant % Elements(k) )
3009
           nBulk = BulkElement % TYPE % NumberOfNodes
3010
           NodeIndexes => BulkElement % NodeIndexes
3011
           BulkElementNodes % x(1:nBulk) = Model % Nodes % x(NodeIndexes)
3012
           BulkElementNodes % y(1:nBulk) = Model % Nodes % y(NodeIndexes)
3013
           BulkElementNodes % z(1:nBulk) = Model % Nodes % z(NodeIndexes)
3014
           IF ( PointInElement( BulkElement,BulkElementNodes, &
3015
               [ SUM( ElementNodes % x(1:n) * Basis), &
3016
               SUM( ElementNodes % y(1:n) * Basis), &
3017
               SUM( ElementNodes % z(1:n) * Basis) ], &
3018
               LocalCoordsInBulkElem) ) EXIT
3019
         END DO
3020
!         PRINT*,'Point in Element: ', LeafQuadrant % Elements(k)
3021
        ELSE
3022
! Go through all BulkElements
3023
! Need MaxElementNodes only in allocation
3024
         nBulk = Model % MaxElementNodes
3025
         ALLOCATE( BulkElementNodes % x( nBulk ),   &
3026
             BulkElementNodes % y( nBulk ),   &
3027
             BulkElementNodes % z( nBulk ) )
3028
         DO k=1,Model % NumberOfBulkElements
3029
           BulkElement => Model % Elements(k)
3030
           nBulk = BulkElement % TYPE % NumberOfNodes
3031
           NodeIndexes => BulkElement % NodeIndexes
3032
           BulkElementNodes % x(1:nBulk) = Model % Nodes % x(NodeIndexes)
3033
           BulkElementNodes % y(1:nBulk) = Model % Nodes % y(NodeIndexes)
3034
           BulkElementNodes % z(1:nBulk) = Model % Nodes % z(NodeIndexes)
3035
           IF ( PointInElement(BulkElement,BulkElementNodes, &
3036
               [ SUM(ElementNodes % x(1:n) * Basis), &
3037
                 SUM(ElementNodes % y(1:n) * Basis), &
3038
                 SUM(ElementNodes % z(1:n) * Basis) ], &
3039
               LocalCoordsInBulkElem) ) EXIT
3040
         END DO
3041
!         PRINT*,'Point in Element: ', k
3042
       END IF
3043
! Calculate value of the function in the bulk element
3044
! Read the integrand from Simulation block in the .sif file
3045
! It is assumed to be a contravariant vector, but
3046
! ListGetRealArray doesn t exist, so we read it component by component
3047
! naming them with suffixes " 1" etc.
3048

3049
      DO j=1,DIM
3050
        WRITE (Component, '(" ",I1.1)') j
3051
        IntegrandFunctionComponent = IntegrandFunctionName(1: &
3052
            LEN_TRIM(IntegrandFunctionName))
3053
        IntegrandFunctionComponent(LEN_TRIM(IntegrandFunctionName)+1: &
3054
            LEN_TRIM(IntegrandFunctionName)+2) = Component
3055
        IntegrandAtNodes(1:nBulk,j) = ListGetReal( Model % Simulation, &
3056
            IntegrandFunctionComponent(1:LEN_TRIM(IntegrandFunctionComponent)), &
3057
            nBulk, NodeIndexes )
3058
      END DO
3059

3060
      DO j=1,DIM
3061
        IntegrandFunction(tg,j) = InterpolateInElement( BulkElement, &
3062
            IntegrandAtNodes(1:nBulk,j),LocalCoordsInBulkElem(1), &
3063
            LocalCoordsInBulkElem(2),LocalCoordsInBulkElem(3) )
3064
      END DO
3065

3066
      DEALLOCATE( BulkElementNodes % x, &
3067
          BulkElementNodes % y, &
3068
          BulkElementNodes % z )
3069

3070
! Get normal to the element boundary at the integration point
3071
! N.B. NormalVector returns covariant normal vector
3072
! NormalVector defined weirdly, doesn t accept LineElement as an argument,
3073
! but wants a pointer
3074
         CurrentElement => LineElement
3075
         Normal = NormalVector( CurrentElement,ElementNodes,ug,vg,.FALSE. )
3076
! Might be consistently in wrong direction, since no check
3077
!         Normal = NormalVector( CurrentElement,ElementNodes,ug,vg,.TRUE. )
3078

3079
! Contravariant vector function to be integrated is already
3080
! at the Gauss point
3081

3082
! Calculate the normal component of the vector function
3083
         FluxAtGPt = SUM( IntegrandFunction(tg,1:DIM) * Normal(1:DIM) )
3084

3085
! Use general coordinate system for dS
3086
! Would be included in s by SqrtMetric
3087
         dS = CoordinateSqrtMetric( SUM( ElementNodes % x(1:n) * Basis), &
3088
             SUM( ElementNodes % y(1:n) * Basis), &
3089
             SUM( ElementNodes % z(1:n) * Basis) )
3090

3091
         Integral = Integral + s*FluxAtGPt*dS
3092

3093
       END DO! of the Gauss integration points
3094

3095
!------------------------------------------------------------------------------
3096
   END FUNCTION LineIntegrate
3097
!------------------------------------------------------------------------------
3098

3099

3100
!------------------------------------------------------------------------------
3101
   FUNCTION ElementArea( Mesh,Element,N ) RESULT(A)
3102
!------------------------------------------------------------------------------
3103
     TYPE(Mesh_t), POINTER :: Mesh
3104
     INTEGER :: N
3105
     TYPE(Element_t) :: Element
3106
!------------------------------------------------------------------------------
3107

3108
     REAL(KIND=dp), TARGET :: NX(N),NY(N),NZ(N)
3109

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

3112
     TYPE(GaussIntegrationPoints_t), TARGET :: IntegStuff
3113
     INTEGER :: N_Integ,t
3114

3115
     REAL(KIND=dp) :: Metric(3,3),Symb(3,3,3),dSymb(3,3,3,3), &
3116
              SqrtMetric,SqrtElementMetric
3117

3118
     TYPE(Nodes_t) :: Nodes
3119

3120
     LOGICAL :: stat
3121

3122
     REAL(KIND=dp) :: Basis(n)
3123
     REAL(KIND=dp) :: dBasisdx(n,3)
3124

3125
     REAL(KIND=dp), DIMENSION(:), POINTER :: U_Integ,V_Integ,W_Integ,S_Integ
3126
!------------------------------------------------------------------------------
3127
 
3128
#if 0
3129
     IF ( ( CurrentCoordinateSystem() == AxisSymmetric .OR. &
3130
            CurrentCoordinateSystem() == CylindricSymmetric ) .AND. Element % TYPE % ELementCode / 100 == 2 ) THEN
3131
       R1 = Mesh % Nodes % x(Element % NodeIndexes(1))
3132
       R2 = Mesh % Nodes % x(Element % NodeIndexes(2))
3133

3134
       Z1 = Mesh % Nodes % y(Element % NodeIndexes(1))
3135
       Z2 = Mesh % Nodes % y(Element % NodeIndexes(2))
3136

3137
       A = PI*ABS(R1+R2)*SQRT((Z1-Z2)*(Z1-Z2)+(R1-R2)*(R1-R2))
3138
     ELSE 
3139
#endif
3140
       Nodes % x => NX
3141
       Nodes % y => NY
3142
       Nodes % z => NZ
3143

3144
       Nodes % x = Mesh % Nodes % x(Element % NodeIndexes)
3145
       Nodes % y = Mesh % Nodes % y(Element % NodeIndexes)
3146
       Nodes % z = Mesh % Nodes % z(Element % NodeIndexes)
3147

3148
       IntegStuff = GaussPoints( element )
3149
       U_Integ => IntegStuff % u
3150
       V_Integ => IntegStuff % v
3151
       W_Integ => IntegStuff % w
3152
       S_Integ => IntegStuff % s
3153
       N_Integ  = IntegStuff % n
3154
!
3155
!------------------------------------------------------------------------------
3156
!   Now we start integrating
3157
!------------------------------------------------------------------------------
3158
!
3159
       A = 0.0
3160
       DO t=1,N_Integ
3161
!
3162
!        Integration stuff
3163
!
3164
         u = U_Integ(t)
3165
         v = V_Integ(t)
3166
         w = W_Integ(t)
3167
!
3168
!------------------------------------------------------------------------------
3169
!        Basis function values & derivatives at the integration point
3170
!------------------------------------------------------------------------------
3171
         stat = ElementInfo( Element,Nodes,u,v,w,SqrtElementMetric, &
3172
                    Basis,dBasisdx )
3173
!------------------------------------------------------------------------------
3174
!        Coordinatesystem dependent info
3175
!------------------------------------------------------------------------------
3176
         IF ( CurrentCoordinateSystem() /= Cartesian ) THEN
3177
           X = SUM( Nodes % x(1:n)*Basis )
3178
           Y = SUM( Nodes % y(1:n)*Basis )
3179
           Z = SUM( Nodes % z(1:n)*Basis )
3180

3181
           SqrtMetric = CoordinateSqrtMetric( x,y,z )
3182
           A =  A + SqrtMetric * SqrtElementMetric * S_Integ(t)
3183
         ELSE
3184
           A =  A + SqrtElementMetric * S_Integ(t)
3185
         END IF
3186
       END DO
3187
#if 0
3188
     END IF
3189
#endif
3190

3191
   END FUNCTION ElementArea
3192
!------------------------------------------------------------------------------
3193

3194

3195
   !------------------------------------------------------------------------------
3196
   !> If element has two of the same indexes regard the element as degenerate.
3197
   !------------------------------------------------------------------------------
3198
   FUNCTION DegenerateElement( Element ) RESULT ( Stat ) 
3199
     TYPE(Element_t), POINTER :: Element
3200
     LOGICAL Stat
3201

3202
     INTEGER :: i,n
3203
     INTEGER, POINTER :: Indexes(:)
3204
     
3205
     Stat = .FALSE.
3206

3207
     n = Element % TYPE % NumberOfNodes
3208
     Indexes => Element % NodeIndexes
3209
     
3210
     DO i = 1, n
3211
       IF( ANY( Indexes(i+1:n) == Indexes(i) ) ) THEN
3212
         Stat = .TRUE.           
3213
         EXIT
3214
       END IF
3215
     END DO
3216
     
3217
   END FUNCTION DegenerateElement
3218

3219
   !------------------------------------------------------------------------------
3220
   !> Return the aspect ratio of an element 
3221
   !------------------------------------------------------------------------------
3222
   FUNCTION ElementAspectRatio(Model, Element ) RESULT ( AspectRatio ) 
3223
     IMPLICIT NONE
3224
     TYPE(Model_t) Model
3225
     TYPE(Element_t) :: Element
3226
     REAL(KIND=dp) :: AspectRatio
3227
     REAL(KIND=dp) :: CharLen(2)
3228

3229
     CharLen = ElementCharacteristicLengths(Model, Element)
3230
     IF (CharLen(1) <= 0) THEN
3231
       AspectRatio = HUGE(AspectRatio)
3232
     ELSE
3233
       AspectRatio = CharLen(2)/CharLen(1)
3234
     END IF
3235
   END FUNCTION ElementAspectRatio
3236

3237
   !------------------------------------------------------------------------------
3238
   !> Return the characteristic lengths of an element 
3239
   !------------------------------------------------------------------------------
3240
   FUNCTION ElementCharacteristicLengths(Model, Element ) RESULT ( Charlengths ) 
3241
     IMPLICIT NONE
3242
     TYPE(Model_t) :: Model
3243
     TYPE(Element_t) :: Element
3244
     REAL(KIND=dp) :: Charlengths(2)
3245
     REAL(KIND=dp) :: Dist
3246

3247
     TYPE(Nodes_t) :: en
3248
     INTEGER :: i,j,n
3249

3250
     INTEGER :: istat
3251
     
3252
     n = Element % TYPE % NumberOfNodes
3253

3254
     ALLOCATE( en % x( n ),   &
3255
               en % y( n ),   &
3256
               en % z( n ), STAT=istat )
3257

3258
     IF( istat /= 0 ) THEN
3259
       CALL Fatal('ElementCharacteristicLengths','Allocation error for ElementNodes')
3260
     END IF
3261

3262
     en % x(1:n) = Model % Nodes % x(Element % NodeIndexes)
3263
     en % y(1:n) = Model % Nodes % y(Element % NodeIndexes)
3264
     en % z(1:n) = Model % Nodes % z(Element % NodeIndexes)
3265
     
3266
     Charlengths = 0._dp
3267
     DO i = 1, n
3268
       DO j = 1, n
3269
         IF (i /= j) THEN
3270
           Dist = SQRT((en % x(i)-en % x(j))**2. + (en % y(i)-en % y(j))**2. + (en % z(i)-en % z(j))**2.)
3271
           IF (Dist < Charlengths(1)) THEN
3272
             Charlengths(1) = Dist
3273
           ELSE IF (Dist > Charlengths(2)) THEN
3274
             Charlengths(2) = Dist
3275
           END IF
3276
         END IF
3277
       END DO
3278
     END DO
3279
     
3280
   END FUNCTION ElementCharacteristicLengths
3281
   
3282
   !------------------------------------------------------------------------------
3283
   !> Return normal of degenerate Element 
3284
   !------------------------------------------------------------------------------
3285
   FUNCTION NormalOfDegenerateElement(Model, Element ) RESULT ( Normal ) 
3286
     IMPLICIT NONE
3287
     TYPE(Model_t) :: Model
3288
     TYPE(Element_t) :: Element
3289
     REAL(KIND=dp) :: a(3), b(3), c(3), Normal(3)
3290

3291
     TYPE(Nodes_t) :: en
3292
     INTEGER :: i,n
3293

3294
     INTEGER :: istat
3295
     
3296
     n = Element % TYPE % NumberOfNodes
3297

3298
     ALLOCATE( en % x( n ),   &
3299
               en % y( n ),   &
3300
               en % z( n ), STAT=istat )
3301

3302
     IF( istat /= 0 ) THEN
3303
       CALL Fatal('NormalOfDegenerateElement','Allocation error for ElementNodes')
3304
     END IF
3305

3306
     en % x(1:n) = Model % Nodes % x(Element % NodeIndexes)
3307
     en % y(1:n) = Model % Nodes % y(Element % NodeIndexes)
3308
     en % z(1:n) = Model % Nodes % z(Element % NodeIndexes)
3309

3310
     a = (/ en % x(1), en % y(1), en % z(1) /)
3311
     b = (/ en % x(2), en % y(2), en % z(2) /)
3312
     c = (/ en % x(n), en % y(n), en % z(n) /)
3313

3314
     Normal = crossproduct(a-b, a-c)
3315

3316
     Normal = Normal / SQRT(SUM(c**2))
3317
     
3318
   END FUNCTION NormalOfDegenerateElement
3319

3320

3321
!------------------------------------------------------------------------------
3322
   FUNCTION FindBoundaryEdgeIndex(Mesh,Boundary,nedge) RESULT(n)
3323
!------------------------------------------------------------------------------
3324
     IMPLICIT NONE
3325
     INTEGER :: n,nedge
3326
     TYPE(Mesh_t), POINTER :: Mesh
3327
     TYPE(Element_t) :: Boundary
3328
!------------------------------------------------------------------------------
3329
     INTEGER :: i,j,k,jb1,jb2,je1,je2
3330
     TYPE(Element_t), POINTER :: Parent, Edge, Face
3331
!------------------------------------------------------------------------------
3332
     n = 0
3333
     SELECT CASE(Boundary % TYPE % ElementCode / 100 )
3334
     CASE(1)
3335
       RETURN
3336
     CASE(2)
3337
       IF ( nedge==1 ) THEN
3338
         Parent => Boundary % BoundaryInfo % Left
3339
         IF ( .NOT. ASSOCIATED(Parent) ) &
3340
             Parent => Boundary % BoundaryInfo % Right
3341

3342
         jb1 = Boundary % NodeIndexes(1)
3343
         jb2 = Boundary % NodeIndexes(2)
3344
         DO i=1,Parent % TYPE % NumberOfEdges
3345
           Edge => Mesh % Edges(Parent % EdgeIndexes(i))
3346
           je1 = Edge % NodeIndexes(1)
3347
           je2 = Edge % NodeIndexes(2)
3348
           IF ( jb1==je1.AND.jb2==je2 .OR. jb1==je2.AND.jb2==je1) EXIT
3349
         END DO
3350
         n = Parent % EdgeIndexes(i)
3351
       END IF
3352
     CASE(3,4)
3353
       j = FindBoundaryFaceIndex(Mesh,Boundary)
3354
       Face => Mesh % Faces(j)
3355
       IF ( nedge>0.AND.nedge<=Face % TYPE % NumberOfEdges ) &
3356
           n = Face % EdgeIndexes(nedge) 
3357
     END SELECT
3358
!------------------------------------------------------------------------------
3359
   END FUNCTION FindBoundaryEdgeIndex
3360
!------------------------------------------------------------------------------
3361

3362

3363
!------------------------------------------------------------------------------
3364
   FUNCTION FindBoundaryFaceIndex(Mesh,Boundary) RESULT(n)
3365
!------------------------------------------------------------------------------
3366
     IMPLICIT NONE
3367
     INTEGER :: n
3368
     TYPE(Element_t) :: Boundary
3369
     TYPE(Mesh_t), POINTER :: Mesh
3370
!------------------------------------------------------------------------------
3371
     INTEGER :: i,j,k,m
3372
     TYPE(Element_t), POINTER :: Parent, Face
3373
!------------------------------------------------------------------------------
3374
     Parent => Boundary % BoundaryInfo % Left
3375
     IF ( .NOT. ASSOCIATED(Parent) ) &
3376
         Parent => Boundary % BoundaryInfo % Right
3377

3378
     DO i=1,Parent % TYPE % NumberOfFaces
3379
       Face => Mesh % Faces(Parent % FaceIndexes(i))
3380
       m = 0
3381
       DO j=1,Face % TYPE % NumberOfNodes
3382
         DO k=1,Boundary % TYPE % NumberOfNodes
3383
           IF ( Face % NodeIndexes(j)==Boundary % NodeIndexes(k)) m=m+1
3384
         END DO
3385
       END DO
3386
       IF ( m==Face % TYPE % NumberOfNodes) EXIT
3387
     END DO
3388
     n = Parent % FaceIndexes(i)
3389
!------------------------------------------------------------------------------
3390
   END FUNCTION FindBoundaryFaceIndex
3391
!------------------------------------------------------------------------------
3392

3393
   
3394
!-----------------------------------------------------------------------------   
3395
!> Given basis function values at surface element find the corresponding local
3396
!> coordinate in the parent element. 
3397
!------------------------------------------------------------------------------
3398
   SUBROUTINE FindParentUVW( Element, n, Parent, np, U, V, W, Basis ) 
3399
!------------------------------------------------------------------------------
3400
     IMPLICIT NONE
3401
     TYPE( Element_t ), POINTER :: Element
3402
     TYPE( Element_t ), POINTER :: Parent
3403
     INTEGER :: n, np
3404
     REAL( KIND=dp ) :: U, V, W, Basis(:)
3405
!------------------------------------------------------------------------------
3406
    INTEGER :: i, j, nParent, check 
3407
    REAL(KIND=dp) :: NodalParentU(n), NodalParentV(n), NodalParentW(n)
3408
!------------------------------------------------------------------------------
3409

3410
    Check = 0
3411

3412
    DO i = 1,n
3413
      DO j = 1,np
3414
        IF( Element % NodeIndexes(i) == Parent % NodeIndexes(j) ) THEN
3415
          Check = Check + 1
3416
          NodalParentU(i) = Parent % Type % NodeU(j)
3417
          NodalParentV(i) = Parent % Type % NodeV(j)
3418
          NodalParentW(i) = Parent % Type % NodeW(j)
3419
        END IF
3420
      END DO
3421
    END DO
3422

3423
    IF( Check /= n ) THEN
3424
      IF(n /= Element % TYPE % NumberOfNodes ) THEN
3425
        CALL Warn('FindParentUVW','Inconsistent size for "n"!')
3426
      END IF
3427
      IF(np /= Parent % TYPE % NumberOfNodes ) THEN
3428
        CALL Warn('FindParentUVW','Inconsistent size for "np"!')
3429
      END IF
3430
      CALL Fatal('FindParentUVW','Could not find all nodes in parent!') 
3431
    END IF
3432
    
3433
    U = SUM( Basis(1:n) * NodalParentU(1:n) )
3434
    V = SUM( Basis(1:n) * NodalParentV(1:n) )
3435
    W = SUM( Basis(1:n) * NodalParentW(1:n) )
3436
!------------------------------------------------------------------------------      
3437
  END SUBROUTINE FindParentUVW
3438
!------------------------------------------------------------------------------      
3439

3440

3441

3442

3443

3444
  
3445

3446
  SUBROUTINE EvaluateVariableAtGivenPoint(No,Values,Mesh,Var,Var2,Var3,Element,LocalCoord,&
3447
      LocalBasis,LocalNode,LocalDGNode,DoGrad,DoDiv,GotEigen,GotEdge,Parent)
3448

3449
    INTEGER :: No
3450
    REAL(KIND=dp) :: Values(:)
3451
    TYPE(Mesh_t), POINTER :: Mesh
3452
    TYPE(Variable_t), POINTER :: Var
3453
    TYPE(Variable_t), POINTER, OPTIONAL :: Var2, Var3
3454
    TYPE(Element_t), POINTER, OPTIONAL :: Element
3455
    REAL(KIND=dp), OPTIONAL :: LocalCoord(3)
3456
    INTEGER, OPTIONAL :: LocalNode
3457
    INTEGER, OpTIONAL :: LocalDGNode
3458
    REAL(KIND=dp), OPTIONAL, TARGET :: LocalBasis(:)
3459
    LOGICAL, OPTIONAL :: DoGrad, DoDiv
3460
    LOGICAL, OPTIONAL :: GotEigen, GotEdge
3461
    TYPE(Element_t), POINTER, OPTIONAL :: Parent
3462
    
3463
    LOGICAL :: Found, EdgeBasis, AVBasis, DgVar, IpVar, ElemVar, DoEigen, &
3464
        PiolaVersion, PElem, NeedDerBasis, Stat, UseGivenNode, IsGrad, IsDiv, &
3465
        IsEigen
3466
    INTEGER :: i1,i2,ii,i,j,k,l,comps, n, n2, nd, np, NoEigenValues, iMode
3467
    INTEGER, TARGET :: DGIndexes(27), Indexes(100), DofIndexes(100), NodeIndex(1)
3468
    INTEGER, POINTER :: pToIndexes(:)
3469
    REAL(KIND=dp) :: u,v,w,detJ
3470
    TYPE(Nodes_t), SAVE :: Nodes
3471
    INTEGER :: lr
3472
    REAL(KIND=dp), TARGET :: Basis(100),dBasisdx(100,3),NodeBasis(1)
3473
    REAL(KIND=dp) :: WBasis(54,3),RotWBasis(54,3)
3474
    REAL(KIND=dp), ALLOCATABLE, SAVE :: fdg(:), fip(:)
3475
    REAL(KIND=dp), POINTER :: pToBasis(:)
3476
    TYPE(Variable_t), POINTER :: pVar
3477
    TYPE(Element_t), POINTER :: Element2
3478
    COMPLEX(KIND=dp), POINTER :: cValues(:)
3479

3480
    INTERFACE 
3481
      SUBROUTINE Ip2DgFieldInElement( Mesh, Parent, nip, fip, np, fdg )
3482
        USE Types
3483
        TYPE(Mesh_t), POINTER :: Mesh
3484
        TYPE(Element_t), POINTER :: Parent
3485
        INTEGER :: nip, np
3486
        REAL(KIND=dp) :: fip(:), fdg(:)
3487
      END SUBROUTINE Ip2DgFieldInElement
3488
    END INTERFACE
3489

3490
    IF(PRESENT(GotEdge)) GotEdge = .FALSE.
3491
    IF(PRESENT(GotEigen)) GotEIgen = .FALSE.
3492
        
3493
    IF(.NOT. ASSOCIATED(Var)) RETURN
3494
    IF(.NOT. ASSOCIATED(Var % Values)) RETURN
3495

3496
    ! Check that a vector field was not given by components 
3497
    comps = 1
3498
    IF(PRESENT(Var2)) THEN
3499
      IF(ASSOCIATED(Var2)) THEN
3500
        comps = 2
3501
        IF( PRESENT(Var3)) THEN
3502
          IF(ASSOCIATED(Var3)) THEN
3503
            comps = 3
3504
          END IF
3505
        END IF
3506
      END IF
3507
    END IF
3508

3509
    IF(.NOT. PRESENT(Element) ) THEN
3510
      CALL Fatal('EvaluteVariableAtGivenPoint','This routine most often really likes to have Element provided!')
3511
    END IF
3512
    
3513
    EdgeBasis = ( Var % TYPE == variable_on_edges )
3514

3515
    DGVar = ( Var % TYPE == variable_on_nodes_on_elements ) 
3516
    IpVar = ( Var % TYPE == variable_on_gauss_points )
3517
    ElemVar = ( Var % TYPE == Variable_on_elements )       
3518

3519
    ! We do not need P-elements if the value is to be found in node since
3520
    ! the higher order p-basis does not have any effect there.
3521
    pElem = .FALSE.
3522
    IF(PRESENT(LocalCoord)) THEN
3523
!    IF(.NOT. PRESENT(LocalNode)) THEN
3524
      IF(ASSOCIATED(Var % Solver)) THEN
3525
        pElem = isActivePElement(Element, Var % Solver) 
3526
      END IF
3527
    END IF
3528
      
3529
    PiolaVersion = .FALSE.
3530
    n = Element % TYPE % NumberOfNodes                  
3531

3532
    ! Elemental variable is a special case which is constant on the element
3533
    ! This does not depend on any other things. 
3534
    IF( ElemVar ) THEN
3535
      l = Element % ElementIndex 
3536
      IF( SIZE( Var % Perm ) >= l ) THEN
3537
        l = Var % Perm(l)
3538
      END IF
3539
      IF( l > 0 ) THEN
3540
        IF( Var % Dofs > 1 ) THEN
3541
          DO ii=1,Var % Dofs
3542
            Values(No+ii) = Var % Values(Var%Dofs*(l-1)+ii)
3543
          END DO
3544
        ELSE
3545
          Values(No+1) = Var % Values(l)              
3546
          IF( comps >= 2 ) Values(No+2) = Var % Values(l)
3547
          IF( comps >= 3 ) Values(No+3) = Var % Values(l)
3548
        END IF
3549
      END IF
3550
      No = No + MAX( Var % Dofs, comps )       
3551
      RETURN
3552
    END IF
3553

3554
    ! We may need to compute dBasisdx if we want derivative or divergence at the point. 
3555
    NeedDerBasis = .FALSE.
3556
    IsGrad = .FALSE.
3557
    IsDiv = .FALSE.
3558
    IF(PRESENT(DoGrad)) IsGrad = DoGrad
3559
    IF(PRESENT(DoDiv)) IsDiv = DoDiv
3560
    NeedDerBasis = IsGrad .OR. IsDiv
3561
    
3562
    pToBasis => NULL()
3563
    pToIndexes => NULL()
3564
    
3565
    IsEigen = .FALSE.
3566
    IF( PRESENT( GotEigen ) ) THEN
3567
      IsEigen = ASSOCIATED( Var % EigenValues )
3568
      IF(IsEigen) NoEigenValues = SIZE( Var % EigenValues )
3569
      IF( comps > 1 .AND. IsEigen ) THEN
3570
        CALL Warn('EvaluetVariableAtGivenPoint','Eigenmode cannot be given in components!')
3571
        IsEigen = .FALSE.
3572
      END IF
3573
      GotEigen = IsEigen
3574
    END IF
3575
    
3576
    ! Given node is the quickest way to estimate the values at nodes.
3577
    ! However, it only applies to nodal fields. 
3578
    NodeIndex(1) = 0
3579
    IF(.NOT. (EdgeBasis .OR. IpVar .OR. ElemVar .OR. NeedDerBasis .OR. pElem) ) THEN      
3580
      IF( DgVar ) THEN        
3581
        IF(PRESENT(LocalDGNode)) NodeIndex(1) = LocalDGNode
3582
        IF(NodeIndex(1) == 0 ) THEN
3583
          IF( PRESENT(LocalNode)) THEN
3584
            PToIndexes => PickDGIndexes(Element,Parent)
3585
            DO i=1, n
3586
              IF( Element % NodeIndexes(i) == LocalNode ) THEN
3587
                NodeIndex(1) = pToIndexes(i)
3588
                EXIT
3589
              END IF
3590
            END DO
3591
          END IF
3592
        END IF
3593
      ELSE IF( PRESENT(LocalNode) ) THEN
3594
        NodeIndex(1) = LocalNode
3595
      END IF
3596

3597
      IF( NodeIndex(1) > 0 ) THEN
3598
        pToIndexes => NodeIndex        
3599
        NodeBasis(1) = 1.0_dp
3600
        pToBasis => NodeBasis
3601
        n = 1
3602
        nd = 1        
3603
      ELSE IF( PRESENT( LocalBasis ) ) THEN
3604
        ! The 2nd quickest is to use existing local nodal basis functions
3605
        PtoBasis => LocalBasis
3606
        nd = n        
3607
        IF( DgVar ) THEN
3608
          PtoIndexes => PickDGIndexes(Element,Parent)
3609
        ELSE
3610
          PtoIndexes => Element % NodeIndexes 
3611
        END IF
3612
      END IF
3613
    END IF    
3614

3615
    ! We don't have basis to evaluate the the fields. Continue to find suitable basis. 
3616
    IF(ASSOCIATED(PtoIndexes) ) THEN
3617
      Element2 => Element
3618
      n2 = n
3619
    ELSE      
3620
      IF(.NOT. PRESENT(LocalCoord) ) THEN
3621
        CALL Fatal('EvaluteVariableAtGivenPoint',&
3622
            'No recipe to evaluate variable without local coordinates: '//TRIM(Var % Name))
3623
      END IF
3624

3625
      u = LocalCoord(1)
3626

3627
      IF( PRESENT(Parent)) THEN
3628
        Element2 => Parent
3629
        n2 = Element2 % TYPE % NumberOfNodes
3630
        IF(.NOT. PRESENT(LocalBasis)) THEN
3631
          CALL Fatal('EvaluateVariableAtGivenPoint','"LocalBasis" is needed when using "Parent" element!')
3632
        END IF
3633
        CALL FindParentUVW( Element, n, Parent, n2, U, V, W, LocalBasis )
3634
      ELSE
3635
        Element2 => Element
3636
        n2 = n
3637
        u = LocalCoord(1)
3638
        v = LocalCoord(2)
3639
        w = LocalCoord(3)
3640
      END IF
3641

3642
      IF( ASSOCIATED( Var % Solver ) ) THEN
3643
        nd = mGetElementDOFs( Indexes, Element2, Var % Solver)
3644
      ELSE
3645
        nd = mGetElementDOFs( Indexes, Element2 )
3646
      END IF          
3647

3648
      IF( EdgeBasis ) THEN
3649
        IF( ListGetLogical(Var % Solver % Values, 'Quadratic Approximation', Found) ) THEN
3650
          PiolaVersion = .TRUE.
3651
        ELSE
3652
          PiolaVersion = ListGetLogical(Var % Solver % Values,'Use Piola Transform', Found )   
3653
        END IF
3654
        np = n2 * Var % Solver % Def_Dofs(Element2 % type % ElementCode/100,Element2 % BodyId,1)
3655
        AVBasis = (np > 0 )
3656
        IF( PRESENT( GotEdge ) ) GotEdge = .TRUE.
3657
      END IF
3658
            
3659
      IF( EdgeBasis ) THEN
3660
        CALL CopyElementNodesFromMesh(Nodes,Mesh,n2,Element2 % NodeIndexes)
3661
        stat = ElementInfo( Element2, Nodes, u, v, w, &
3662
            detJ, Basis, dBasisdx,  EdgeBasis = WBasis, &
3663
            RotBasis = RotWBasis, USolver = Var % Solver)
3664
      ELSE          
3665
        IF( pElem ) THEN
3666
          ! The standard element of SaveLine, SaveScalars etc. is most likely standard nodal element.
3667
          ! If the user gives something else we may be trouble...
3668
          !CALL CopyElementNodesFromMesh(Nodes,Mesh,n2,Element2 % NodeIndexes)
3669
          PtoIndexes => Indexes
3670
          CALL CopyElementNodesFromMesh(Nodes,Mesh,nd,pToIndexes)
3671
          CALL ConvertToPReference(Element2 % TYPE % ElementCode,u,v,w)            
3672
          IF( NeedDerBasis ) THEN
3673
            stat = ElementInfo( Element2, Nodes, u, v, w, detJ, Basis, dBasisdx, &
3674
                USolver = Var % Solver )
3675
          ELSE
3676
            stat = ElementInfo( Element2, Nodes, u, v, w, detJ, Basis, &
3677
                USolver = Var % Solver )
3678
          END IF
3679
        ELSE
3680
          IF( NeedDerBasis ) THEN
3681
            CALL CopyElementNodesFromMesh(Nodes,Mesh,n2,Element2 % NodeIndexes)
3682
            stat = ElementInfo( Element2, Nodes, u, v, w, detJ, Basis, dBasisdx )
3683
          ELSE
3684
            CALL CopyElementNodesFromMesh(Nodes,Mesh,n2,Element2 % NodeIndexes)          
3685
            stat = ElementInfo( Element2, Nodes, u, v, w, detJ, Basis )
3686
          END IF
3687
        END IF
3688

3689
        IF( DgVar ) THEN
3690
          PtoIndexes => PickDGIndexes(Element2)
3691
        END IF
3692
      END IF
3693
      IF(.NOT. ASSOCIATED(pToIndexes)) PtoIndexes => Indexes
3694
      IF(.NOT. ASSOCIATED(pToBasis)) PtoBasis => Basis
3695
    END IF
3696
    
3697
    IF( EdgeBasis ) THEN
3698
      DofIndexes(1:nd) = Var % Perm(PToIndexes(1:nd))
3699
      IF( IsEigen ) THEN
3700
        DO iMode = 1, NoEigenValues
3701
          DO j=1,3          
3702
            No = No + 1
3703
            IF( ALL(DofIndexes(np+1:nd) > 0 ) ) THEN            
3704
              cValues => Var % EigenVectors(iMode,:)
3705
              Values(No) = SUM( WBasis(1:nd-np,j) * cValues(DofIndexes(np+1:nd)))
3706
            ELSE
3707
              Values(No) = 0.0_dp
3708
            END IF
3709
          END DO
3710
        END DO
3711
      ELSE
3712
        DO j=1,3
3713
          No = No + 1
3714
          IF( ALL(DofIndexes(np+1:nd) > 0 ) ) THEN
3715
            Values(No) = SUM( WBasis(1:nd-np,j) * Var % Values(DofIndexes(np+1:nd)))
3716
          ELSE
3717
            Values(No) = 0.0_dp
3718
          END IF
3719
        END DO
3720
      END IF
3721
      IF( AVBasis ) THEN
3722
        No = No + 1
3723
        Values(No) = 0.0_dp
3724
        IF( ALL(DofIndexes(1:np) > 0 ) ) THEN
3725
          Values(No) = SUM( Basis(1:np) * Var % Values(DofIndexes(1:np)))
3726
        END IF
3727
      END IF
3728

3729
    ELSE IF ( IpVar ) THEN
3730
      l = 0 
3731
      IF( SIZE( Var % Perm ) > Element2 % ElementIndex ) THEN 
3732
        i1 = Var % Perm(Element2 % ElementIndex)
3733
        i2 = Var % Perm(Element2 % ElementIndex+1)
3734
        l = i2-i1
3735
      END IF
3736

3737
      IF(l>0) THEN
3738
        IF( .NOT. ALLOCATED(fip)) THEN
3739
          ALLOCATE( fip(l) )
3740
        ELSE IF (SIZE(fip) < l) THEN
3741
          DEALLOCATE( fip )
3742
          ALLOCATE( fip(l) )
3743
        END IF
3744

3745
        IF(.NOT. ALLOCATED(fdg)) THEN
3746
           ALLOCATE( fdg(n) )
3747
        ELSE IF(SIZE(fdg) < n) THEN
3748
           DEALLOCATE( fdg )
3749
           ALLOCATE( fdg(n) )
3750
        END IF
3751

3752
        DO ii=1,MAX(Var % Dofs,comps)
3753
          IF( Var % Dofs > 1 ) THEN
3754
            CONTINUE
3755
          ELSE          
3756
            IF( ii == 1 ) THEN
3757
              pVar => Var
3758
            ELSE IF( ii == 2 ) THEN
3759
              pVar => Var2
3760
            ELSE IF( ii == 3 ) THEN
3761
              pVar => Var3
3762
            END IF
3763
            fip(1:l) = pVar % Values(i1+1:i2)
3764
          END IF
3765

3766
          CALL Ip2DgFieldInElement( Mesh, Element2, l, fip, n2, fdg )              
3767
          Values(No+ii) = SUM( PtoBasis(1:n) * fdg(1:n) )
3768
        END DO
3769
      END IF
3770
      
3771
      No = No + MAX(Var % Dofs, Comps )      
3772
    ELSE
3773
      IF(.NOT. ASSOCIATED(pToBasis)) THEN
3774
        CALL Fatal('EvaluteVariableAtGivenPoint',&
3775
            'pToBasis not associated for variable: '//TRIM(Var % Name))
3776
      END IF
3777
      IF(.NOT. ASSOCIATED(pToIndexes)) THEN
3778
        CALL Fatal('EvaluteVariableAtGivenPoint',&
3779
            'pToIndexes not associated for variable: '//TRIM(Var % Name))
3780
      END IF
3781

3782
      ! This maybe should be checked earlier, but better late than never perhaps ?
3783
      IF (DGVar .AND. pElem ) nd = Element2 % Type % NumberOfNodes
3784
      
3785
      IF( ASSOCIATED(Var % Perm) ) THEN
3786
        DofIndexes(1:nd) = Var % Perm(PToIndexes(1:nd))
3787
      ELSE
3788
        DofIndexes(1:nd) = PtoIndexes(1:nd)
3789
      END IF
3790

3791

3792
      
3793
        IF( IsGrad ) THEN          
3794
          IF( Var % Dofs /= 1 ) THEN
3795
            CALL Fatal('EvaluteVariableAtGivenPoint','Gradient only possible for one dof!')
3796
          END IF
3797
          IF( ALL(Dofindexes(1:nd) > 0 ) ) THEN
3798
            DO ii=1,3
3799
              Values(No+ii) = SUM( dBasisdx(1:nd,ii) * Var % Values(DofIndexes(1:nd)))
3800
            END DO
3801
          ELSE
3802
            Values(No+1:No+3) = 0._dp
3803
          END IF
3804
          No = No + 3
3805
        ELSE IF( IsEigen ) THEN          
3806
          DO iMode = 1, NoEigenValues
3807
            cValues => Var % EIgenVectors(iMode,:)            
3808
            IF( ALL(Dofindexes(1:nd) > 0 ) ) THEN
3809
              DO ii=1,Var % DOfs              
3810
                Values(No+ii) = Values(No+ii) + SUM( PtoBasis(1:nd) * &
3811
                    cValues(Var%Dofs*(DofIndexes(1:nd)-1)+ii))
3812
              END DO
3813
            ELSE
3814
              Values(No+1:No+Var % Dofs)=0._dp      
3815
            END IF
3816
          END DO
3817
          No = No + Var % Dofs
3818
        ELSE
3819
          IF( ALL(Dofindexes(1:nd) > 0 ) ) THEN
3820
            IF( Var % Dofs > 1 ) THEN
3821
              DO ii=1,Var % Dofs
3822
                Values(No+ii) = SUM( PtoBasis(1:nd) * &
3823
                    Var % Values(Var%Dofs*(DofIndexes(1:nd)-1)+ii))
3824
              END DO
3825
            ELSE
3826
              Values(No+1) = SUM(PToBasis(1:nd) * Var % Values(DofIndexes(1:nd)))            
3827
              IF( comps >= 2 ) Values(No+2) = SUM(PToBasis(1:nd) * Var2 % Values(DofIndexes(1:nd)))
3828
              IF( comps >= 3 ) Values(No+3) = SUM(PToBasis(1:nd) * Var3 % Values(DofIndexes(1:nd)))            
3829
            END IF
3830
          ELSE
3831
            Values(No+1:No+MAX( Var % Dofs, comps ))=0._dp
3832
          ENDIF
3833
          No = No + MAX( Var % Dofs, comps )
3834
        END IF
3835
      END IF
3836

3837
#if 0
3838
    ! Debugging code
3839
    IF( .NOT. EdgeBasis .AND. ASSOCIATED(Var % Perm) .AND. No > 0 .AND. ASSOCIATED( PToIndexes ) ) THEN
3840
      BLOCK
3841
        REAL(KIND=dp) :: vmin, vmax
3842
        INTEGER :: kmax
3843
        PRINT *,'VarInfo:',TRIM(Var % Name), nd, n, pToIndexes(1:nd)
3844

3845
        kmax = 1
3846
        DO i=1,nd
3847
          IF( PtoBasis(i) > pToBasis(kmax)) kmax = i
3848
        END DO
3849

3850
        IF( ASSOCIATED(pToIndexes) ) THEN
3851
          DofIndexes(1:nd) = Var % Perm(PToIndexes(1:nd))
3852
        ELSE
3853
          DofIndexes(1:nd) = PtoIndexes(1:nd)
3854
        END IF
3855
          
3856
        vmin = MINVAL(Var % Values(DofIndexes(1:nd)),DofIndexes(1:nd)>0)
3857
        vmax = MAXVAL(Var % Values(DofIndexes(1:nd)),DofIndexes(1:nd)>0)
3858

3859
        PRINT *,'PtoIndexes:',PtoIndexes(1:nd)
3860
        PRINT *,'PtoBasis',SUM(PToBasis(1:nd)),PToBasis(1:nd)
3861
        PRINT *,'ElemRange:',vmin< Values(1) .AND. vmax > Values(1), &
3862
            vmin, vmax, Var % Values(Var % Perm(PToIndexes(kmax))), Values(1:No)
3863
      END BLOCK
3864
    END IF
3865
#endif
3866
      
3867
  CONTAINS
3868

3869
    FUNCTION PickDgIndexes(Element,PParent) RESULT ( PToInds) 
3870
      TYPE(Element_t), POINTER :: Element
3871
      INTEGER, POINTER :: PtoInds(:)
3872
      TYPE(Element_t), POINTER, OPTIONAL :: PParent
3873

3874
      TYPE(Element_t), POINTER :: Parent
3875
      INTEGER :: i,j,lr
3876
            
3877
      IF( ASSOCIATED( Element % DgIndexes ) ) THEN
3878
        PToInds => Element % DGIndexes
3879
      ELSE IF( ASSOCIATED( Element % BoundaryInfo ) ) THEN
3880
        DO lr=1,2
3881
          IF(lr==1) THEN
3882
            Parent => Element % BoundaryInfo % Left
3883
          ELSE
3884
            Parent => Element % BoundaryInfo % Right
3885
          END IF
3886
          IF(.NOT. ASSOCIATED( Parent ) ) CYCLE
3887
          IF(.NOT. ASSOCIATED( Parent % DGIndexes ) ) CYCLE
3888
          IF( PRESENT(PParent)) THEN
3889
            IF(.NOT. ASSOCIATED(Parent, PParent) ) CYCLE
3890
          END IF          
3891
          IF( ALL( Var % Perm( Parent % DGIndexes ) /= 0) ) THEN                  
3892
            DO i=1,Element % TYPE % NumberOfNodes
3893
              DO j=1,Parent % TYPE % NumberOfNodes
3894
                IF( Element % NodeIndexes(i) == Parent % NodeIndexes(j) ) THEN
3895
                  DGIndexes(i) = Parent % DGIndexes(j)
3896
                  EXIT
3897
                END IF
3898
              END DO
3899
            END DO
3900
            EXIT
3901
          END IF
3902
        END DO
3903
        PtoInds => DGIndexes
3904
      END IF
3905
    END FUNCTION PickDgIndexes
3906
    
3907
  END SUBROUTINE EvaluateVariableAtGivenPoint
3908

3909

3910

3911

3912
!------------------------------------------------------------------------------
3913
   SUBROUTINE mGetBoundaryIndexesFromParent( Mesh, Element, Indexes, indSize )
3914
!------------------------------------------------------------------------------
3915
     IMPLICIT NONE
3916

3917
     ! Parameters
3918
     TYPE(Mesh_t) :: Mesh
3919
     TYPE(Element_t), POINTER :: Element
3920
     INTEGER :: indSize, Indexes(:)
3921
     
3922
     ! Variables
3923
     TYPE(Element_t), POINTER :: Edge, Face
3924
     INTEGER :: i,j,n
3925
     TYPE(Element_t), POINTER :: Parent
3926
     
3927
     ! Clear indexes
3928
     Indexes = 0
3929
     indSize = 0
3930

3931
     Parent => Element % pDefs % localParent
3932
     IF ( .NOT. ASSOCIATED(Parent) ) RETURN
3933
             
3934
     n = Element % TYPE % NumberOfNodes
3935

3936
     ! Nodal indexes
3937
     Indexes(1:n) = Element % NodeIndexes(1:n)
3938
     indSize = n
3939
     
3940
     
3941
     ! Assign rest of indexes if necessary
3942
     SELECT CASE(Parent % TYPE % DIMENSION)
3943
     CASE (1)
3944
       CONTINUE
3945

3946
     CASE (2)
3947
       IF(.NOT. ASSOCIATED(Mesh % Edges) ) RETURN
3948
         
3949
       ! Add index for each bubble dof in edge
3950
        DO i=1,Element % BDOFs
3951
           n = n+1
3952
           
3953
           IF (SIZE(Indexes) < n) THEN
3954
              CALL Warn('mGetBoundaryIndexes','Not enough space reserved for indexes')
3955
              RETURN
3956
           END IF
3957

3958
           Indexes(n) = Mesh % NumberOfNodes + &
3959
                (Parent % EdgeIndexes(Element % PDefs % localNumber)-1) * Mesh % MaxEdgeDOFs + i
3960
        END DO
3961
     
3962
        indSize = n 
3963
      CASE (3)
3964
        IF(.NOT. ( ASSOCIATED(Mesh % Faces) .AND. ASSOCIATED(Mesh % Edges) ) ) RETURN
3965
        IF(.NOT. ASSOCIATED(Element % PDefs) ) RETURN        
3966
        IF(Element % PDefs % LocalNumber == 0 ) RETURN
3967
                
3968
        ! Get boundary face
3969
        Face => Mesh % Faces( Parent % FaceIndexes(Element % PDefs % localNumber) )
3970
        
3971
        ! Add indexes of faces edges 
3972
        DO i=1, Face % TYPE % NumberOfEdges
3973
           Edge => Mesh % Edges( Face % EdgeIndexes(i) )
3974
           
3975
           ! If edge has no dofs jump to next edge
3976
           IF (Edge % BDOFs <= 0) CYCLE
3977

3978
           DO j=1,Edge % BDOFs
3979
              n = n + 1
3980
              
3981
              IF (SIZE(Indexes) < n) THEN
3982
                 CALL Warn('mGetBoundaryIndexes','Not enough space reserved for indexes')
3983
                 RETURN
3984
              END IF
3985
              
3986
              Indexes(n) = Mesh % NumberOfNodes +&
3987
                  ( Face % EdgeIndexes(i)-1)*Mesh % MaxEdgeDOFs + j
3988
           END DO
3989
        END DO
3990
               
3991
        ! Add indexes of faces bubbles
3992
        DO i=1,Face % BDOFs
3993
           n = n + 1
3994

3995
           IF (SIZE(Indexes) < n) THEN
3996
              CALL Warn('mGetBoundaryIndexes','Not enough space reserved for indexes')
3997
              RETURN
3998
           END IF
3999

4000
           Indexes(n) = Mesh % NumberOfNodes + &
4001
                Mesh % NumberOfEdges * Mesh % MaxEdgeDOFs + &
4002
                (Parent % FaceIndexes( Element % PDefs % localNumber )-1) * Mesh % MaxFaceDOFs + i
4003
        END DO        
4004

4005
        indSize = n
4006
     CASE DEFAULT
4007
        CALL Fatal('mGetBoundaryIndexes','Unsupported dimension')
4008
     END SELECT
4009
!------------------------------------------------------------------------------
4010
   END SUBROUTINE mGetBoundaryIndexesFromParent
4011
!------------------------------------------------------------------------------
4012

4013
END MODULE ElementUtils
4014

4015
!> \} ElmerLib
4016

4017
Product

Resources

Company
