CoCalc -- Adaptive.F90

GitHub Repository: ElmerCSC/elmerfem
Path: blob/devel/fem/src/Adaptive.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, Mikko Lyly
27
! *  Email:   [email protected], [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: Autumn 2000
34
! *
35
! *****************************************************************************/
36

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

40
!--------------------------------------------------------------------------------------------------------
41
!> Module for adaptive meshing routines. The adaptivity is based on solver-specific error indicators that
42
!> are used to create a field with the desired mesh density. This may be used by some mesh generators to
43
!> create a more optimal mesh. 
44
!--------------------------------------------------------------------------------------------------------
45
MODULE Adaptive
46

47
  USE GeneralUtils
48
  USE SolverUtils, ONLY : VectorValuesRange
49
  USE ElementUtils, ONLY : ElementArea
50
  USE ModelDescription
51
  USE MeshUtils, ONLY : AllocateMesh, AllocatePDefinitions, FindMeshEdges, &
52
      LoadMesh2, MeshStabParams, PrepareMesh, ReleaseMesh, &
53
      ReleaseMeshEdgeTables, ReleaseMeshFaceTables, SetActiveElementsTable, &
54
      SetCurrentMesh, TransferCoordAndTime, UpdateSolverMesh, WriteMeshToDisk, &
55
      WriteMeshToDisk2
56
  USE MeshRemeshing
57
  USE SaveUtils, ONLY : SaveGmshOutput
58
  USE DefUtils, ONLY: GetMaterial, GetReal, GetBodyForce, GetSolverParams, GetLogical, &
59
      GetNofActive, GetActiveElement, GetElementNofBDOFs, GetBoundaryElement, &
60
      ActiveBoundaryElement, GetNofBoundaryElements, GetElementFamily, GetElementNodes
61
  USE ElementDescription, ONLY: GetEdgeMap, mGetElementDOFs
62
  USE MainUtils, ONLY : AddEquationSolution
63
  
64
  IMPLICIT NONE
65

66
  
67
CONTAINS
68

69
!------------------------------------------------------------------------------
70
  SUBROUTINE RefineMesh( Model,Solver,Quant,Perm, &
71
            InsideResidual, EdgeResidual, BoundaryResidual )
72
!------------------------------------------------------------------------------
73
    IMPLICIT NONE
74

75
    TYPE( Model_t ) :: Model
76
    TYPE(Solver_t), TARGET :: Solver
77
    REAL(KIND=dp) :: Quant(:)
78
    INTEGER :: Perm(:)
79

80
    INTERFACE
81
       SUBROUTINE BoundaryResidual( Model,Edge,Mesh,Quant,Perm,Gnorm,Indicator )
82
          USE Types
83
          TYPE(Element_t), POINTER :: Edge
84
          TYPE(Model_t) :: Model
85
          TYPE(Mesh_t), POINTER :: Mesh
86
          REAL(KIND=dp) :: Quant(:), Indicator(2), Gnorm
87
          INTEGER :: Perm(:)
88
       END SUBROUTINE BoundaryResidual
89

90

91
       SUBROUTINE EdgeResidual( Model,Edge,Mesh,Quant,Perm, Indicator)
92
          USE Types
93
          TYPE(Element_t), POINTER :: Edge
94
          TYPE(Model_t) :: Model
95
          TYPE(Mesh_t), POINTER :: Mesh
96
          REAL(KIND=dp) :: Quant(:), Indicator(2)
97
          INTEGER :: Perm(:)
98
       END SUBROUTINE EdgeResidual
99

100

101
       SUBROUTINE InsideResidual( Model,Element,Mesh,Quant,Perm,Fnorm, Indicator)
102
          USE Types
103
          TYPE(Element_t), POINTER :: Element
104
          TYPE(Model_t) :: Model
105
          TYPE(Mesh_t), POINTER :: Mesh
106
          REAL(KIND=dp) :: Quant(:), Indicator(2), Fnorm
107
          INTEGER :: Perm(:)
108
       END SUBROUTINE InsideResidual
109
    END INTERFACE
110
!------------------------------------------------------------------------------
111

112
    TYPE(Mesh_t), POINTER :: RefMesh, NewMesh, Mesh, sMesh
113
    TYPE( Nodes_t ) :: Nodes
114
    TYPE(Solver_t), POINTER :: SolverPtr
115
    INTEGER, POINTER :: Permutation(:)
116
    LOGICAL, POINTER :: EdgeSplitted(:)
117
    INTEGER, POINTER :: Referenced(:)
118
    TYPE( Element_t ), POINTER :: RefElement
119
    INTEGER :: i,j,k,n,nn,MarkedElements
120
    TYPE( Variable_t ), POINTER :: Var, Var1, NewVar, RTFlux
121
    REAL(KIND=dp) :: MaxError, ErrorLimit, minH, maxH, MaxChangeFactor, &
122
        LocalIndicator,ErrorEstimate,t,TotalTime,RemeshTime,s,FinalRef,&
123
        MaxScale,AveScale,OutFrac,MaxFrac, NodalMaxError, mError, SolNorm
124

125
    LOGICAL :: BandwidthOptimize, Found, Coarsening, &
126
        MeshNumbering, DoFinalRef
127
    INTEGER :: MaxDepth, MinDepth, NLen, MeshDim
128
    CHARACTER(:), ALLOCATABLE :: Path, VarName
129
    REAL(KIND=dp), POINTER  :: Time(:), NodalError(:), PrevValues(:), &
130
         Hvalue(:), HValue1(:), PrevNodalError(:), PrevHValue(:), hConvergence(:), ptr(:), tt(:)
131
    REAL(KIND=dp), POINTER  :: ErrorIndicator(:), eRef(:), hRef(:), Work(:)
132
    LOGICAL :: AdaptiveInterp, Parallel, AdaptiveOutput, AdaptInit
133
    LOGICAL :: WithRecovery, ConvCond
134
    TYPE(ValueList_t), POINTER :: Params
135
    CHARACTER(*), PARAMETER :: Caller = 'RefineMesh'
136
    REAL(KIND=dp), POINTER :: Wrk(:,:)
137
    REAL(KIND=dp) :: CoordScale(3)
138
    
139
    SAVE DoFinalRef
140
    
141
!---------------------------------------------------------------------------------
142
!
143
!   Initialize:
144
!   -----------
145
    CALL Info( Caller, ' ', Level=5 )
146
    CALL Info( Caller, &
147
        '----------- M E S H   R E F I N E M E N T --------------', Level=5 )
148
    TotalTime = CPUTime()
149
    RemeshTime = 0.0d0
150

151
    RefMesh => Solver % Mesh
152

153
    IF( RefMesh % DiscontMesh ) THEN
154
      CALL Fatal(Caller,'Adaptive refinement not possible for discontinuous mesh!')
155
    END IF
156

157
    Params => Solver % Values
158
    SolverPtr => Solver
159
    
160
    MinDepth = ListGetInteger( Params, 'Adaptive Min Depth', Found )
161

162
    MaxDepth = ListGetInteger( Params, 'Adaptive Max Depth', Found )
163
    IF (Found) THEN
164
      IF ( Refmesh % AdaptiveDepth > MaxDepth ) THEN
165
        ! Setting this flag will override convergence for this equation on steady-state level.
166
        Solver % Mesh % AdaptiveFinished = .TRUE.
167
        CALL Info( Caller,'Max adaptive depth reached! Doing nothing!', Level = 5 )
168
        RETURN
169

170
        IF( MinDepth > MaxDepth ) THEN
171
          CALL Warn(Caller,'"Adaptive Min Depth" greater than Max!' )
172
        END IF
173
      END IF
174
    END IF
175

176
    AdaptInit = ( RefMesh % AdaptiveDepth == 0 ) 
177

178
    IF( AdaptInit ) THEN
179
      CALL Info(Caller,'Initializing stuff on the coarsest level!')
180
      DoFinalRef = .FALSE.
181
    END IF
182

183
    IF( DoFinalRef ) THEN
184
      CALL Info( Caller, 'Final refinement done. Nothing to do!', Level=6 )
185
      RefMesh % OutputActive = .TRUE.      
186
      RefMesh % Parent % OutputActive = .FALSE.      
187
      CALL Info(Caller,'Setting adaptive restart to True!',Level=12)
188
      RefMesh % AdaptiveFinished = .TRUE.
189
      RETURN
190
    ELSE
191
      RefMesh % OutputActive = .TRUE.
192
      RefMesh % AdaptiveFinished = .FALSE.
193
    END IF
194
        
195
    
196
    ! Interpolation is costly in parallel. Do it by default only in serial. 
197
    Parallel = ( ParEnv % PEs > 1 )
198
    AdaptiveInterp = ListGetLogical( Params,'Adaptive Interpolate',Found,DefValue=.TRUE. )
199
!   IF(.NOT. Found) AdaptiveInterp = Parallel 
200
    
201
    AdaptiveOutput = ListGetLogical( Params,'Adaptive Output',Found )
202
    
203
    DO i=1,RefMesh % NumberOfBulkElements
204
       RefMesh % Elements(i) % Splitted = 0
205
    END DO
206

207
!   Compute the local error indicators:
208
!   -----------------------------------
209
    t = CPUTime()
210
    CALL AllocateVector( ErrorIndicator, RefMesh % NumberOfBulkElements )
211

212
    WithRecovery = ListGetLogical(Params, 'Flux Recovery', Found) 
213
    IF (WithRecovery) THEN
214
      CALL FluxRecovery(Model, Solver, RefMesh, ErrorIndicator, MaxError)
215
      RTFlux => VariableGet(RefMesh % Variables, 'RTFlux')
216
    ELSE
217
      MaxError = ComputeError( Model, ErrorIndicator, RefMesh, &
218
          Quant, Perm, InsideResidual, EdgeResidual, BoundaryResidual )
219
    END IF
220
    !
221
    !   Global error estimate:
222
    !   ----------------------
223
    ErrorEstimate = SUM( ErrorIndicator**2) 
224
    IF (Parallel) ErrorEstimate = ParallelReduction(ErrorEstimate)
225
    
226
    IF (WithRecovery) THEN
227
      ErrorEstimate =  SQRT(ErrorEstimate)
228
    ELSE
229
      n = SIZE(ErrorIndicator)
230
      IF (Parallel) n = ParallelReduction(n)
231
      ErrorEstimate =  SQRT( ErrorEstimate / n )
232
    END IF
233
    
234
    WRITE( Message, * ) 'Error computation time (cpu-secs):               ',CPUTime()-t
235
    CALL Info( Caller, Message, Level = 6 )
236
    
237
    IF(ListGetLogical(Params,'Adaptive Error Histogram',Found ) ) THEN
238
      CALL ShowVectorHistogram(ErrorIndicator,SIZE(ErrorIndicator))
239
    END IF
240
    
241
    WRITE( Message, * ) 'Max error      =                                 ',MaxError
242
    CALL Info( Caller, Message, Level = 6 )
243
    WRITE( Message, * ) 'Error estimate =                                 ',ErrorEstimate
244
    CALL Info( Caller, Message, Level = 6 )
245
    !WRITE(12,*) RefMesh % NumberOfBulkElements,ErrorEstimate,MaxError
246

247
!
248
!   Add nodal average of the h-value to the mesh variable list:
249
!   -----------------------------------------------------------
250

251
    nn = RefMesh % NumberOfNodes
252
    Var => VariableGet( RefMesh % Variables, 'Hvalue', ThisOnly=.TRUE. )
253

254
    IF ( ASSOCIATED( Var ) ) THEN
255
      Hvalue => Var % Values
256
      Var % PrimaryMesh => RefMesh
257
      IF( AdaptInit ) Hvalue = 0.0_dp
258
    ELSE
259
      CALL AllocateVector( Hvalue, nn )
260

261
      CALL VariableAdd( RefMesh % Variables, RefMesh, Solver, &
262
          'Hvalue', 1, Hvalue, Output = AdaptiveOutput )       
263

264
      Var => VariableGet( RefMesh % Variables, 'Hvalue', ThisOnly=.TRUE. )      
265
      IF(.NOT. ASSOCIATED(Var) ) THEN
266
        CALL Fatal(Caller,'Could not add variable Hvalue?')
267
      END IF
268
      Hvalue = 0.0d0
269
    END IF
270

271
    CALL AllocateVector( PrevHvalue, nn )
272
    IF( AdaptInit ) THEN
273
      PrevHValue(1:nn) = 0.0_dp
274
    ELSE
275
      PrevHvalue(1:nn) = Hvalue(1:nn)
276
    END IF
277

278
    CALL AllocateVector( Referenced, nn )
279

280
    Hvalue = 0.0d0
281
    Referenced = 0
282
    n = RefMesh % MaxElementNodes
283
    ALLOCATE( Nodes % x(n), Nodes % y(n), Nodes % z(n) )
284

285
    ! Average nodal Hvalue from the bulk elements
286
    DO i=1,RefMesh % NumberOfBulkElements
287
       RefElement => RefMesh % Elements(i)
288
       n = RefElement % TYPE % NumberOfNodes
289

290
       Nodes % x(1:n) = RefMesh % Nodes % x(RefElement % NodeIndexes)
291
       Nodes % y(1:n) = RefMesh % Nodes % y(RefElement % NodeIndexes)
292
       Nodes % z(1:n) = RefMesh % Nodes % z(RefElement % NodeIndexes)
293
       s = ElementDiameter( RefElement, Nodes )
294
       DO j=1,n
295
          k = RefMesh % Elements(i) % NodeIndexes(j)
296
          Hvalue(k) = Hvalue(k) + s
297
          Referenced(k) = Referenced(k) + 1
298
       END DO
299
    END DO
300

301
    DEALLOCATE( Nodes % x, Nodes % y, Nodes % z )
302

303
    WHERE( Referenced(1:nn) > 0 )
304
      Hvalue(1:nn) = Hvalue(1:nn) / Referenced(1:nn)
305
    END WHERE
306
    CALL ParallelAverageHvalue(  RefMesh, Hvalue )
307
    
308
!   Add estimate of the convergence with respect to h:
309
!  ----------------------------------------------------
310
    Var => VariableGet( RefMesh % Variables, 'hConvergence', ThisOnly=.TRUE. )
311

312
    IF ( ASSOCIATED( Var ) ) THEN
313
      hConvergence => Var % Values
314
      Var % PrimaryMesh => RefMesh
315
      IF( AdaptInit ) hConvergence = 1.0_dp
316
    ELSE
317
      CALL AllocateVector( hConvergence, nn )
318
      hConvergence = 1.0d0
319
      CALL VariableAdd( RefMesh % Variables, RefMesh, Solver, &
320
          'hConvergence', 1, hConvergence, Output=AdaptiveOutput )
321
      Var => VariableGet( RefMesh % Variables, 'hConvergence', ThisOnly=.TRUE. )
322
    END IF
323

324
!   Add nodal average of the computed estimate for the
325
!   solution error to the mesh variable list:
326
!   --------------------------------------------------
327
    VarName = GetVarName(Solver % Variable)
328
    nlen = LEN_TRIM(VarName)
329
    Var => VariableGet( RefMesh % Variables, &
330
         VarName(1:nlen)  // '.error', ThisOnly=.TRUE. )
331

332
    IF ( ASSOCIATED( Var ) ) THEN
333
       NodalError  => Var % Values
334
       Var % PrimaryMesh => RefMesh
335
    ELSE
336
       CALL AllocateVector( NodalError, nn )
337
       CALL VariableAdd( RefMesh % Variables, RefMesh, Solver, &
338
          VarName(1:nlen) // '.error', 1, NodalError )
339
    END IF
340

341
    Var => VariableGet( RefMesh % Variables, &
342
         VarName(1:nlen) // '.perror', ThisOnly=.TRUE. )
343

344
    IF ( ASSOCIATED( Var ) ) THEN
345
      PrevNodalError  => Var % Values
346
      Var % PrimaryMesh => RefMesh
347
      IF(AdaptInit) PrevNodalError = 0.0_dp
348
    ELSE
349
      CALL AllocateVector( PrevNodalError, RefMesh % NumberOfNodes )
350
      PrevNodalError = 0.0d0
351
      CALL VariableAdd( RefMesh % Variables, RefMesh, Solver, &
352
          VarName(1:nlen) // '.perror', 1, PrevNodalError, Output=AdaptiveOutput)
353
    END IF
354

355
    NodalError = 0.0d0
356
    Referenced = 0
357
    DO i = 1, RefMesh % NumberOfBulkElements
358
       DO j=1,RefMesh % Elements(i) % TYPE % NumberOfNodes
359
          k = RefMesh % Elements(i) % NodeIndexes(j)
360
          Referenced(k) = Referenced(k) + 1
361
          NodalError(k) = NodalError(k) + ErrorIndicator(i)
362
       END DO
363
    END DO
364

365
    WHERE( Referenced(1:nn) > 0 )
366
       NodalError(1:nn) = NodalError(1:nn) / Referenced(1:nn)
367
    END WHERE
368

369
    IF (ParEnv % PEs>1) THEN
370
      CALL ParallelAverageHvalue(  RefMesh, NodalError )
371

372
      NodalMaxError = ParallelReduction( MAXVAL(NodalError),2 )
373
      WRITE( Message, * ) 'Nodal Max error      =                           ',NodalMaxError
374
      CALL Info( Caller, Message, Level = 6 )
375
    END IF
376

377
!
378
!   Smooth error, if requested:
379
!   ---------------------------
380
    k = ListGetInteger( Params, 'Adaptive Pre Smoothing', Found )
381
    IF ( Found .AND. k > 0 ) THEN 
382
       CALL AllocateVector( eRef, nn )
383
       DO j=1,k
384
          eRef(1:nn) = NodalError(1:nn)
385
          Referenced = 0
386
          NodalError = 0
387
          DO i=1,RefMesh % NumberOfBulkElements
388
             n = RefMesh % Elements(i) % TYPE % NumberOfNodes
389
             NodalError(RefMesh % Elements(i) % NodeIndexes) = &
390
                NodalError(RefMesh % Elements(i) % NodeIndexes) + &
391
                   SUM( eRef(RefMesh % Elements(i) % NodeIndexes) ) / n
392
             Referenced( RefMesh % Elements(i) % NodeIndexes ) = &
393
                Referenced( RefMesh % Elements(i) % NodeIndexes ) + 1
394
          END DO
395
          WHERE( Referenced(1:nn) > 1 )
396
             NodalError(1:nn) = NodalError(1:nn) / Referenced(1:nn)
397
          END WHERE
398
       END DO
399
       DEALLOCATE( eRef )
400
    END IF
401

402
    DEALLOCATE( Referenced )
403
!
404
!   Add reference error to variable list:
405
!   -------------------------------------
406
    Var => VariableGet( RefMesh % Variables, &
407
         VarName(1:nlen) // '.eRef', ThisOnly=.TRUE. )
408

409
    IF ( ASSOCIATED( Var ) ) THEN
410
      eRef => Var % Values
411
      Var % PrimaryMesh => RefMesh
412
      IF( AdaptInit ) eRef(1:nn) = NodalError(1:nn)
413
    ELSE
414
      CALL AllocateVector( eRef, nn )
415
      eRef(1:nn) = NodalError(1:nn)      
416
      CALL VariableAdd( RefMesh % Variables, RefMesh, Solver, &
417
          VarName(1:nlen) // '.eRef',1,eRef, Output=AdaptiveOutput )
418
    END IF
419
!
420
!   Mesh projection may alter the values somewhat!
421
!   ----------------------------------------------
422
    eRef = MAX( eRef, 1.0d-12 )
423

424
!
425
!   Add reference h to variable list:
426
!   ---------------------------------
427
    Var => VariableGet( RefMesh % Variables, 'hRef', ThisOnly=.TRUE. )
428

429
    IF ( ASSOCIATED( Var ) ) THEN
430
      hRef => Var % Values
431
      Var % PrimaryMesh => RefMesh
432
      IF( AdaptInit ) hRef(1:nn) = Hvalue(1:nn)
433
    ELSE
434
      CALL AllocateVector( hRef, nn )
435
      hRef(1:nn) = Hvalue(1:nn)
436
      CALL VariableAdd( RefMesh % Variables, RefMesh, Solver, &
437
          'hRef', 1, hRef, Output=AdaptiveOutput)
438
    END IF
439
!
440
!   Mesh projection may alter the values somewhat!
441
!   ----------------------------------------------
442
    hRef = MAX( hRef, 1.0d-12 )
443

444
!   Check for convergence:
445
!   ----------------------
446
    ErrorLimit = ListGetConstReal( Params,'Adaptive Error Limit', Found )
447
    IF ( .NOT.Found ) ErrorLimit = 0.5d0
448
    IF (WithRecovery) THEN
449
      ConvCond = ErrorEstimate < ErrorLimit
450
      n = SIZE(ErrorIndicator)
451
      IF (Parallel) n = ParallelReduction(n)
452
      ! The target element error if the error were evenly distributed:
453
      ErrorLimit =  ErrorLimit / SQRT(REAL(n))
454
    END IF
455

456
    MaxScale = ListGetConstReal( Params,'Adaptive Max Error Scale', Found )
457
    IF(.NOT. Found) MaxScale = 1.0_dp
458

459
    AveScale = ListGetConstReal( Params,'Adaptive Average Error Scale', Found )
460
    IF(.NOT. Found) AveScale = 1.0_dp
461

462
    MaxFrac = ListGetConstReal( Params,'Adaptive Max Outlier Fraction',Found ) 
463
    IF( Found ) THEN
464
      OutFrac = OutlierFraction(ErrorIndicator,SIZE(ErrorIndicator),ErrorLimit)     
465
      WRITE( Message, * ) 'Outlier frac.  =                                 ',OutFrac
466
      CALL Info( Caller, Message, Level = 6 )
467
    ELSE
468
      ! Just make values that are omitted
469
      MaxFrac = 1.0_dp
470
      OutFrac = 0.0_dp
471
    END IF
472
            
473
    !PRINT *,'Check for convergence:'
474
    !PRINT *,'MaxError, MaxScale*ErrorLimit:', MaxError,MaxScale*ErrorLimit 
475
    !PRINT *,'MaxError < MaxScale*ErrorLimit :', MaxError < MaxScale*ErrorLimit
476
    !PRINT *,'ErrorEstimate, AveScale*ErrorLimit', ErrorEstimate, AveScale*ErrorLimit
477
    !PRINT *,'ErrorEstimate < AveScale*ErrorLimit', ErrorEstimate < AveScale*ErrorLimit
478
    !PRINT *,'OutFrac:',OutFrac,MaxFrac, OutFrac < MaxFrac
479
    !PRINT *,'Depth:',RefMesh % AdaptiveDepth, MinDepth
480
    
481
    IF( RefMesh % AdaptiveDepth > MinDepth ) THEN
482
      mError = MaxError
483
      IF(ListGetLogical(Params, 'Adaptive Use Nodal Error As Limit', Found)) mError = NodalMaxError
484
      IF (.NOT. WithRecovery) THEN
485
        ConvCond = ErrorEstimate < AveScale * ErrorLimit
486
      END IF
487
        
488
      IF ( mError < MaxScale * ErrorLimit .AND. ConvCond .AND. OutFrac < MaxFrac ) THEN
489
        FinalRef = ListGetConstReal( Params,'Adaptive Final Refinement', DoFinalRef ) 
490
        IF(DoFinalRef ) THEN      
491
          CALL Info( Caller, 'Performing one final refinement',Level=6)
492
          ErrorLimit = FinalRef * ErrorLimit 
493
        ELSE
494
          CALL Info( Caller, 'Mesh convergence limit reached. Nothing to do!', Level=6 )
495
          RefMesh % Parent % OutputActive = .FALSE.      
496
          RefMesh % AdaptiveFinished = .TRUE.
497
          IF (WithRecovery) RTFlux % SteadyConverged = 1 
498
          GOTO 10
499
        END IF
500
      END IF
501
    END IF
502

503
!
504
!   Get additional parameters:
505
!   --------------------------
506
    minH = ListGetConstReal( Params, 'Adaptive Min H', Found )
507
    maxH = ListGetConstReal( Params, 'Adaptive Max H', Found )
508

509
    MaxChangeFactor = ListGetConstReal( Params, &
510
            'Adaptive Max Change', Found )
511
    IF ( .NOT.Found .OR. MaxChangeFactor <= AEPS ) MaxChangeFactor = 3.0d0
512

513
    Coarsening = ListGetLogical( Params, 'Adaptive Coarsening', Found )
514
    IF( .NOT.Found ) Coarsening = .TRUE.
515
!
516
!   Compute local convergence of the solution with respect to h:
517
!   ------------------------------------------------------------
518

519
    WHERE( eRef(1:nn) > 0 )
520
      PrevNodalError(1:nn) = PrevNodalError(1:nn) + &
521
         LOG( HValue(1:nn) / hRef(1:nn) ) * LOG( NodalError(1:nn) / eRef(1:nn) )
522
    END WHERE
523
    CALL ParallelAverageHvalue(  RefMesh, PrevNodalError )
524

525
    PrevHvalue(1:nn) = PrevHvalue(1:nn) + LOG( HValue(1:nn) / hRef(1:nn) )**2
526
    CALL ParallelAverageHvalue(  RefMesh, PrevHvalue )
527

528
    IF ( RefMesh % AdaptiveDepth > 0 ) THEN
529
       WHERE( PrevHValue(1:nn) > 0 )
530
          hConvergence(1:nn)  = MAX( PrevNodalError(1:nn) / PrevHValue(1:nn), 0.25d0 )
531
       ELSEWHERE
532
          hConvergence(1:nn)  = 0.25d0
533
       END WHERE
534
    END IF
535
    CALL ParallelAverageHvalue(  RefMesh, hConvergence )
536

537
!   Generate the new mesh:
538
!   ----------------------
539
    IF ( ListGetLogical( Params, 'Adaptive Remesh', Found ) ) THEN
540
      t = RealTime()
541
      IF( ListGetLogical( Params,'Adaptive Remesh Use MMG', Found ) ) THEN
542
#ifdef HAVE_MMG
543
        CALL Info(Caller,'Using MMG library for mesh refinement', Level=5)        
544
        NewMesh => MMG_ReMesh( RefMesh, ErrorLimit/3, HValue, &
545
            NodalError, hConvergence, minH, maxH, MaxChangeFactor, Coarsening )         
546
#else
547
        CALL Fatal( Caller,'Remeshing requested with MMG but not compiled with!')
548
#endif          
549
      ELSE       
550
        CALL Info(Caller,'Using file I/O for mesh refinement',Level=5)
551
        NewMesh => External_ReMesh( RefMesh, ErrorLimit/3, HValue, &
552
            NodalError, hConvergence, minH, maxH, MaxChangeFactor, Coarsening )
553
      END IF
554
      RemeshTime = RealTime() - t
555
      WRITE( Message, * ) 'Remeshing time (real-secs):                      ',RemeshTime
556
      CALL Info( Caller, Message, Level=6 )
557
    ELSE
558
      NewMesh => SplitMesh( RefMesh, ErrorIndicator, ErrorLimit, &
559
          NodalError, Hvalue, hConvergence, minH, maxH, MaxChangeFactor )
560
    END IF
561

562
    Hvalue(1:nn) = PrevHValue(1:nn)
563
!   NodalError = PrevNodalError
564

565
    IF ( .NOT.ASSOCIATED( NewMesh ) ) THEN
566
      CALL Info( Caller,'Current mesh seems fine. Nothing to do.', Level=6 )
567
      IF (ASSOCIATED(RefMesh % Parent)) RefMesh % Parent % OutputActive = .FALSE.
568
      IF (WithRecovery) RTFlux % SteadyConverged = 1
569
      GOTO 10
570
    ELSE
571
      CALL PrepareMesh(Model, NewMesh, ParEnv % PEs>1)
572
    END IF
573

574
    CALL Info( Caller,'The new mesh consists of: ', Level=5 )
575
    CALL Info( Caller,'Nodal points: '&
576
        //TRIM(I2S(NewMesh % NumberOfNodes)),Level=5)
577
    CALL Info( Caller,'Bulk elements: '&
578
        //TRIM(I2S(NewMesh % NumberOfBulkElements)),Level=5)
579
    CALL Info( Caller,'Boundary elements: '&
580
        //TRIM(I2S(NewMesh % NumberOfBoundaryElements)),Level=5)
581

582
!-------------------------------------------------------------------
583

584
!   All the mesh geometry related tables are ready now,
585
!   next we update model and solver related tables:
586
!   ----------------------------------------------------
587

588
    t = CPUTime()
589

590
!   Add the new mesh to the global list of meshes:
591
!   ----------------------------------------------
592
    NewMesh % Next   => Model % Meshes 
593
    Model % Meshes   => NewMesh
594
    RefMesh % Child  => NewMesh
595
    NewMesh % Parent => RefMesh
596
    NewMesh % Child => NULL()
597

598
    NewMesh % Name = ListGetString( Params, &
599
         'Adaptive Mesh Name', Found )
600
    IF ( .NOT. Found ) NewMesh % Name = 'RefinedMesh'
601

602
    MeshNumbering = ListGetLogical( Params, &
603
        'Adaptive Mesh Numbering', Found )
604
    IF(.NOT. Found ) MeshNumbering = .TRUE.
605
    
606
    NewMesh % AdaptiveDepth = RefMesh % AdaptiveDepth + 1
607
    IF( MeshNumbering ) THEN
608
      NewMesh % Name = TRIM( NewMesh % Name ) // I2S(NewMesh % AdaptiveDepth)
609
    END IF
610

611
    IF ( ListGetLogical( Params, 'Adaptive Save Mesh', Found ) ) THEN 
612
      Nlen = LEN_TRIM(OutputPath)
613
      IF ( Nlen > 0 ) THEN
614
        Path = OutputPath(1:Nlen) // '/' // TRIM(NewMesh % Name)
615
      ELSE
616
        Path = TRIM(NewMesh % Name)
617
      END IF
618
      CALL MakeDirectory( TRIM(path) // CHAR(0) )
619
    
620
      IF( ParEnv % PEs > 1 ) THEN
621
        CALL WriteMeshToDisk2( Model, NewMesh, Path, ParEnv % MyPe )
622
      ELSE
623
        CALL WriteMeshToDisk( NewMesh, Path )
624
      END IF
625
    END IF
626
    
627
!   Initialize local variables for the new mesh:
628
!   --------------------------------------------
629
    NULLIFY( NewMesh % Variables )
630
    
631
    CALL TransferCoordAndTime( RefMesh, NewMesh ) 
632

633
    IF (WithRecovery) THEN
634
      ! The following calls SetCurrentMesh( CurrentModel, NewMesh ),
635
      ! adds the variable of the solver so that its associated with
636
      ! the new mesh
637
      CALL UpdateSolverMesh( Solver, NewMesh, .TRUE. )
638
      ! The following is needed so that the "variable loads" is also created
639
      CALL AddEquationSolution(SolverPtr, Solver % TimeOrder > 0)
640
    ELSE
641
      CALL SetCurrentMesh( Model, NewMesh )
642
    END IF
643
    
644
    ! Initialize the field variables for the new mesh. These are
645
    ! interpolated from the old meshes variables. Vector variables
646
    ! are in the variable lists in two ways: as vectors and as
647
    ! vector components. We MUST update the vectors (i.e. DOFs>1)
648
    ! first!!!!!
649
    ! -----------------------------------------------------------
650
    CALL Info(Caller,'Interpolate vectors from old mesh to new mesh!',Level=7)    
651
    Var => RefMesh % Variables
652
    DO WHILE( ASSOCIATED( Var ) )
653
      ! This cycles global variable such as time etc. 
654
      IF( SIZE( Var % Values ) == Var % DOFs ) THEN
655
        Var => Var % Next
656
        CYCLE
657
      END IF
658

659
      IF ( Var % DOFs > 1 ) THEN
660
        NewVar => VariableGet( NewMesh % Variables,Var % Name,.FALSE. )
661
        k = SIZE( NewVar % Values )
662
        IF ( ASSOCIATED( NewVar % Perm ) ) THEN
663
          k = COUNT( NewVar % Perm > 0 )
664
        END IF
665
        IF ( GetVarName( NewVar ) == 'flow solution' ) THEN
666
          NewVar % Norm = 0.0d0
667
          DO i=1,NewMesh % NumberOfNodes
668
            DO j=1,NewVar % DOFs-1
669
              NewVar % Norm = NewVar % Norm + &
670
                  NewVar % Values( NewVar % DOFs*(i-1)+j )**2
671
            END DO
672
          END DO
673
          NewVar % Norm = SQRT( NewVar % Norm / k )
674
        ELSE
675
          NewVar % Norm = SQRT( SUM(NewVar % Values**2) / k )
676
        END IF
677
      END IF
678
      Var => Var % Next
679
    END DO
680
    CALL Info(Caller,'Interpolation to new mesh done!',Level=20)    
681

682
!   Second time around, update scalar variables and
683
!   vector components:
684
!   -----------------------------------------------
685
    CALL Info(Caller,'Interpolate scalars from old mesh to new mesh!',Level=7)    
686
    Var => RefMesh % Variables
687
    DO WHILE( ASSOCIATED( Var ) )
688

689
      ! This cycles global variable such as time etc. 
690
      IF( SIZE( Var % Values ) == Var % DOFs ) THEN
691
        Var => Var % Next
692
        CYCLE
693
      END IF
694

695
      SELECT CASE( Var % Name )
696

697
      CASE( 'coordinate 1', 'coordinate 2', 'coordinate 3' )
698
        CONTINUE
699

700
      CASE('rtflux')
701
        !
702
        ! We skip the recovery variable
703
        !
704
        CONTINUE
705
        
706
      CASE DEFAULT
707
        IF (WithRecovery .AND. Var % Name == Solver % Variable % Name .OR. &
708
            WithRecovery .AND. Var % Name == Solver % Variable % Name // ' ' // 'loads') THEN
709
          ! The variable of the solver should already exist, just check this:
710
          NewVar => VariableGet( NewMesh % Variables, Var % Name, .TRUE., UnfoundFatal = .TRUE. )
711
          IF (Var % Name == Solver % Variable % Name) THEN
712
            NewVar % PrevNorm = Var % Norm
713
          END IF
714
        ELSE IF ( Var % DOFs == 1 ) THEN
715
          
716
          ! Skip the fields related to adaptivity since they are specific to each mesh 
717
          Found = .FALSE.
718
!          Found = Found .OR. INDEX( Var % Name, 'ave test' ) > 0 
719
          Found = Found .OR. INDEX( Var % Name, '.error'  ) > 0
720
          Found = Found .OR. INDEX( Var % Name, '.eref'   ) > 0
721
          Found = Found .OR. INDEX( Var % Name, '.perror' ) > 0
722
          IF ( Found ) THEN
723
            k = Solver % Variable % NameLen
724
            IF ( Var % Name(1:k) /= Solver % Variable % Name(1:k) ) THEN
725
              Var => Var % Next
726
              CYCLE
727
            END IF
728
          END IF
729

730
          IF( .NOT. AdaptiveInterp ) THEN
731
            ! Add field without interpolation
732
            NewVar => VariableGet( NewMesh % Variables, Var % Name, .TRUE. )
733
            IF(.NOT. ASSOCIATED(NewVar) ) THEN
734
              CALL VariableAddVector( NewMesh % Variables, NewMesh, Solver,&
735
                  Var % Name,Var % DOFs)
736
              NewVar => VariableGet( NewMesh % Variables, Var % Name, .TRUE. )
737
            END IF
738
            NewVar % PrevNorm = Var % Norm
739
          ELSE
740
            ! Interpolate scalar variables using automatic internal interpolation 
741
            NewVar => VariableGet( NewMesh % Variables, Var % Name, .FALSE. )
742
            k = SIZE(NewVar % Values)
743
            IF ( ASSOCIATED(NewVar % Perm) ) THEN
744
              k = COUNT(NewVar % Perm > 0)
745
            END IF
746
            NewVar % Norm = SQRT(SUM(NewVar % Values**2)/k)
747
          END IF
748
        END IF
749
      END SELECT
750
      Var => Var % Next
751
    END DO
752

753
!-------------------------------------------------------------------    
754
    WRITE( Message, * ) 'Mesh variable update time (cpu-secs):            ',CPUTime()-t
755
    CALL Info( Caller, Message, Level = 6 )
756
!-------------------------------------------------------------------    
757

758
!
759
!   Update Solver structure to use the new mesh:
760
!   ---------------------------------------------    
761
    CALL MeshStabParams( NewMesh )
762
!
763
!   Nothing computed on this mesh yet:
764
!   ----------------------------------
765
    NewMesh % SavesDone    = 0  ! start new output file
766
    NewMesh % OutputActive = .FALSE.
767
    NewMesh % Changed   = .TRUE.
768

769
!
770
!   Create matrix structures for the new mesh:
771
!   ------------------------------------------    
772
    t = CPUTime()
773

774
!
775
!   Try to account for the reordering of DOFs
776
!   due to bandwidth optimization:
777
    !   -----------------------------------------
778
    
779
    CALL Info(Caller,'Updating solver mesh to reflect the new adaptive mesh!',Level=12)
780

781
    IF (.NOT. WithRecovery) THEN
782
      CALL UpdateSolverMesh( Solver, NewMesh, .TRUE. )
783
    END IF
784
    
785
    CALL SetActiveElementsTable( Model, Solver )
786
          
787
    CALL ParallelInitMatrix( Solver, Solver % Matrix )
788

789
    IF (WithRecovery) THEN
790
      CALL Info(Caller, 'UpdateSolverMesh for RTFlux solver', Level=12)
791
      CALL UpdateSolverMesh(RTFlux % Solver, NewMesh, .TRUE.)
792
      CALL Info(Caller, 'UpdateSolverMesh ready', Level=12)
793
      CALL SetActiveElementsTable( Model, RTFlux % Solver )
794
    END IF
795
      
796
    WRITE( Message, * ) 'Matrix structures update time (cpu-secs):        ',CPUTime()-t
797
    CALL Info( Caller, Message, Level=6 )
798

799
!
800
!   Release previous meshes. Keep only the original mesh, and
801
!   the last two meshes:
802
!   ---------------------------------------------------------
803
    n = 0
804
    Mesh => RefMesh % Parent
805
    DO WHILE( ASSOCIATED(Mesh) )
806
      sMesh => Mesh % Parent
807
      n = n+1
808
      IF ( Mesh % AdaptiveDepth /= 0 ) THEN
809
        IF ( ASSOCIATED( Mesh % Parent ) ) THEN
810
          Mesh % Parent % Child => Mesh % Child                        
811
        END IF
812

813
        IF ( ASSOCIATED(Mesh % Child) ) THEN
814
          Mesh % Child % Parent => Mesh % Parent
815
          ! Eliminate the mesh to be released also from here!
816
          Mesh % Child % Next => Mesh % Next 
817
        END IF
818

819
        CALL Info(Caller,'Releasing mesh: '//TRIM(Mesh % Name),Level=8)
820
        CALL ReleaseMesh( Mesh )
821

822
        !Some solvers should be able to go from Mesh to its child!
823
        !Therefore do not nullify the pointers that enable to do that!
824
        !Mesh % Child  => NULL()
825
        !Mesh % Parent => NULL()
826
      END IF
827

828
      Mesh => sMesh
829
    END DO
830

831
!------------------------------------------------------------------------------
832

833
10  CONTINUE
834

835
    IF (.NOT. WithRecovery) THEN
836
      !   Comment the next calls, if you want to keep the edge tables:
837
      !   ------------------------------------------------------------
838
      IF(.NOT.isPelement(RefMesh % Elements(1))) THEN
839
        CALL ReleaseMeshEdgeTables( RefMesh )
840
        CALL ReleaseMeshFaceTables( RefMesh )
841
      END IF
842
    END IF
843
    
844
    IF (.NOT. ASSOCIATED(Model % Mesh, RefMesh)) CALL SetCurrentMesh( Model, RefMesh )
845
    DEALLOCATE( ErrorIndicator, PrevHvalue )
846
    
847
    IF ( RemeshTime > 0 ) THEN
848
      WRITE( Message, * ) 'Mesh refine took in total (cpu-secs):  ', CPUTIme() - TotalTime 
849
      CALL Info( Caller, Message, Level=6 )
850
      WRITE( Message, * ) 'Remeshing took in total (real-secs):   ',RemeshTime
851
      CALL Info( Caller, Message, Level=6 )
852
    END IF
853
    CALL Info( Caller,'----------- E N D   M E S H   R E F I N E M E N T --------------', Level=5 )
854
    
855

856
CONTAINS
857

858
  SUBROUTINE ShowVectorHistogram(x,n)
859
    REAL(KIND=dp) :: x(:)
860
    INTEGER :: n
861

862
    REAL(KIND=dp) :: xmin,xmax,dx
863
    INTEGER :: ncoh
864
    INTEGER, ALLOCATABLE :: cohcnt(:)
865

866
    ! Just for now do it only for 1st partition
867
    ! Later make things parallel.
868
    IF( ParEnv % MyPe /= 0) RETURN
869
        
870
    ncoh = 20
871
    ALLOCATE(cohcnt(ncoh))
872
    cohcnt = 0
873

874
    ! Find extremum errors
875
    xmin = MINVAL(x(1:n))
876
    xmax = MAXVAL(x(1:n))
877
    dx = ( xmax-xmin) / ncoh
878

879
    ! Count each error cohort
880
    DO i=1,n
881
      j = CEILING( (x(i)-xmin)/dx )
882
      j = MAX( 1, MIN( j, ncoh ) )
883
      cohcnt(j) = cohcnt(j) + 1
884
    END DO
885

886
    PRINT *,'Histogram for vector:'
887
    DO i=1,ncoh
888
      IF( cohcnt(i) > 0 ) THEN
889
        PRINT *,i,': ',cohcnt(i),' (',xmin+(i-1)*dx,' to ',xmin+i*dx,')'
890
      END IF
891
    END DO
892
       
893
  END SUBROUTINE ShowVectorHistogram
894

895

896
  ! What fraction of value in x are above xlim? 
897
  FUNCTION OutlierFraction(x,n,xlim) RESULT (f)
898
    REAL(KIND=dp) :: x(:)
899
    REAL(KIND=dp) :: xlim, f
900
    INTEGER :: n
901
    INTEGER :: np, nlim
902

903
    np = n
904
    nlim = COUNT(x(1:n) > xlim)
905

906
    IF( Parallel ) THEN
907
      np = ParallelReduction(np)
908
      nlim = ParallelReduction(nlim)
909
    END IF
910

911
    f = 1.0_dp*nlim/n
912

913
  END FUNCTION OutlierFraction
914

915
  
916
  SUBROUTINE ComputeDesiredHvalue( RefMesh, ErrorLimit, HValue, NodalError, &
917
      hConvergence, minH, maxH, MaxChange, Coarsening ) 
918
!------------------------------------------------------------------------------
919
    REAL(KIND=dp) :: NodalError(:), hConvergence(:), &
920
        ErrorLimit, minH, maxH, MaxChange, HValue(:)
921
    LOGICAL :: Coarsening
922
    TYPE(Mesh_t), POINTER :: RefMesh
923
!------------------------------------------------------------------------------
924
    INTEGER :: i,j,k,n
925
    REAL(KIND=dp) :: Lambda, hLimitScale
926
    INTEGER, ALLOCATABLE :: Hcount(:)
927
    TYPE(Matrix_t), POINTER :: A
928
!------------------------------------------------------------------------------
929

930
    hLimitScale = ListGetConstReal( Params,'Adaptive H Limit Scale', Found )
931
    IF ( .NOT.Found ) hLimitScale = 1.0d0
932

933
    
934
    DO i=1,RefMesh % NumberOfNodes      
935
      IF ( NodalError(i) < 100*AEPS ) CYCLE 
936

937
      Lambda = ( ErrorLimit / NodalError(i) ) ** ( 1.0d0 / hConvergence(i) )
938

939
      IF ( RefMesh % AdaptiveDepth < 1 ) THEN
940
        Lambda = HValue(i) * MAX( MIN( Lambda, 1.33d0), 0.75d0)
941
      ELSE
942
        Lambda = HValue(i) * MAX(MIN(Lambda, MaxChange), 1.0d0/MaxChange)
943
      END IF
944

945
      IF( .NOT.Coarsening ) Lambda = MIN( Lambda, Hvalue(i) )
946

947
      IF ( maxH > 0 ) Lambda = MIN( Lambda, maxH )
948
      IF ( minH > 0 ) Lambda = MAX( Lambda, minH )
949

950
      HValue(i) = Lambda * hLimitScale
951
    END DO
952

953
    IF(.NOT. ListCheckPresent(CurrentModel % Solver % Values,'Adaptive Element Count') ) RETURN
954

955
    
956
    BLOCK
957
      TYPE(Element_t), POINTER :: Element
958
      TYPE(Nodes_t) :: Nodes
959
      REAL(KIND=dp), ALLOCATABLE, SAVE :: Basis(:)
960
      REAL(KIND=dp) :: h, Weight, detJ, CountInteg, CountInteg0=-1.0, &
961
          CountDesired, Cfix, HScale, Vol
962
      TYPE(GaussIntegrationPoints_t) :: IP
963
      INTEGER :: dim,i,t
964
      LOGICAL :: stat
965
      INTEGER :: TimesVisited = 0
966

967
      SAVE CountInteg0, Cfix, Nodes, TimesVisited
968
      
969
      n = RefMesh % MaxElementNodes      
970
      IF ( .NOT. ASSOCIATED( Nodes % x ) ) THEN
971
        ALLOCATE( Nodes % x(n), Nodes % y(n),Nodes % z(n), Basis(n) )
972
      END IF
973
      
974
      dim = RefMesh % MeshDim
975
      TimesVisited = TimesVisited + 1
976

977
      IF(TimesVisited == 1) THEN
978
        ! These fits are experimental ones in simple geometry using MMG2D and MMG3D
979
        ! rounded to the closest integer. 
980
        IF( dim == 2 ) THEN
981
          Cfix = 2.0_dp
982
        ELSE
983
          Cfix = 12.0_dp
984
        END IF
985
      ELSE
986
        Cfix = Cfix * RefMesh % NumberOfBulkElements / CountInteg0
987
      END IF
988
              
989
      CountInteg = 0.0_dp
990
      Vol = 0.0_dp
991
      
992
      DO i=1,RefMesh % NumberOfBulkElements
993
        Element => RefMesh % Elements(i)
994
        
995
        n = Element % TYPE % NumberOfNodes
996
        Nodes % x(1:n) = RefMesh % Nodes % x(Element % NodeIndexes(1:n))
997
        Nodes % y(1:n) = RefMesh % Nodes % y(Element % NodeIndexes(1:n))
998
        Nodes % z(1:n) = RefMesh % Nodes % z(Element % NodeIndexes(1:n))
999
        
1000
        IP = GaussPoints( Element )
1001
        DO t=1,IP % n
1002
          stat = ElementInfo( Element, Nodes, IP % U(t), IP % V(t), &
1003
              IP % W(t), detJ, Basis )          
1004
          Weight = IP % s(t) * DetJ
1005
          
1006
          h = SUM(Basis(1:n)*HValue(Element % NodeIndexes(1:n)))          
1007

1008
          CountInteg = CountInteg + Weight * h**(-dim)
1009
          Vol = Vol + Weight
1010
        END DO
1011
      END DO
1012

1013
      CountInteg = Cfix * CountInteg
1014
      
1015
      CountDesired = ListGetFun(CurrentModel % Solver % Values,'Adaptive Element Count',&
1016
          1.0_dp*TimesVisited,Found)      
1017
      IF( Found .AND. CountDesired > 0.0_dp ) THEN
1018
        Hscale = (CountInteg/CountDesired)**(1.0_dp/dim)
1019
        HValue = Hscale * Hvalue
1020
        CountInteg0 = CountDesired
1021
      ELSE
1022
        Hscale = 1.0_dp
1023
        CountInteg0 = CountInteg
1024
      END IF
1025
      
1026
      IF( InfoActive(10)) THEN
1027
        PRINT *,'AdaptiveCount: ',CountInteg0, CountInteg, RefMesh % NumberOfBulkElements, &
1028
            Cfix, Hscale, Vol
1029
      END IF
1030
        
1031
    END BLOCK
1032

1033

1034
    
1035
  END SUBROUTINE ComputeDesiredHvalue
1036

1037

1038
! Compute the desired Hvalue at the interface where the adaptive error computation currently
1039
! fails. This way parallel adaptivity can be done in some way at least...
1040
!-------------------------------------------------------------------------------------------
1041
  SUBROUTINE ParallelAverageHvalue( RefMesh, HValue ) 
1042
!------------------------------------------------------------------------------
1043
    REAL(KIND=dp) :: HValue(:)
1044
    TYPE(Mesh_t), POINTER :: RefMesh
1045

1046
    INTEGER :: i,j,k,l,n,minnei,maxnei
1047
    INTEGER, POINTER :: p(:)
1048
    INTEGER, ALLOCATABLE :: Hcount(:), ip(:)
1049
    TYPE(Matrix_t), POINTER :: A
1050
!------------------------------------------------------------------------------
1051
  
1052
    IF( ParEnv % PEs == 1 ) RETURN
1053

1054
    ALLOCATE(Hcount(SIZE(Hvalue)))
1055
    Hcount = 0
1056
    
1057
    A => Solver % Matrix
1058
    p => Solver % Variable % Perm
1059

1060
    ALLOCATE(ip(A% NumberOfRows)); ip=0
1061

1062
    DO i=1,RefMesh % NumberOfNodes
1063
      j = p(i)
1064
      IF(j<=0) CYCLE
1065
      ip(j) = i
1066
    END DO
1067

1068
    DO i=1,RefMesh % NumberOfNodes
1069
      ! Deal only with interface nodes here
1070
      j = p(i)
1071
      IF(j<=0) CYCLE
1072
      IF(A % ParallelInfo % GInterface(j)) THEN
1073
        Hvalue(i) = 0.0_dp
1074
        ! Go through all connected nodes
1075
        DO l=A % Rows(j),A % Rows(j+1)-1
1076
          k = A % Cols(l)
1077
          
1078
          ! Skip oneself and other interface nodes
1079
          IF(A % ParallelInfo % GInterface(k)) CYCLE
1080

1081
          k = ip(k)
1082
          IF(k<=0 .OR. i==k) CYCLE          
1083
          
1084
          ! Add the observation
1085
          Hvalue(i) = Hvalue(i) + Hvalue(k)
1086
          Hcount(i) = Hcount(i) + 1
1087
        END DO
1088
      END IF
1089
    END DO
1090
    
1091
    ! Perform parallel summation, only interface gets summed. 
1092
    BLOCK
1093
      INTEGER, ALLOCATABLE :: Xcount(:)
1094
      REAL(KIND=dp), ALLOCATABLE :: Xvalue(:)
1095
      ALLOCATE( Xcount(A % NumberOfRows), XValue(A % NumberOfRows) )
1096
      Xcount = 0; Xvalue = 0
1097

1098
      DO i=1,RefMesh % NumberOfNodes
1099
        j = p(i)
1100
        IF (j<=0) CYCLE
1101
        Xvalue(j) = HValue(i)
1102
        Xcount(j) = HCount(i)
1103
      END DO
1104

1105
      CALL ParallelSumVector( A, Xvalue )
1106
      CALL ParallelSumVectorInt( A, Xcount ) 
1107

1108
      DO i=1,RefMesh % NumberOfNodes
1109
        j = p(i)
1110
        IF (j<=0) CYCLE
1111
        Hvalue(i) = Xvalue(j)
1112
        Hcount(i) = Xcount(j)
1113
      END DO
1114
    END BLOCK
1115

1116

1117
    maxnei = MAXVAL( Hcount )
1118
    minnei = MINVAL( Hcount, Hcount > 0 ) 
1119

1120
    maxnei = ParallelReduction(maxnei,2)
1121
    minnei = ParallelReduction(minnei,1)     
1122

1123
    CALL Info('ParallelAverageHvalue','Averaging count range is ['//TRIM(I2S(minnei)) &
1124
        //','//TRIM(I2S(maxnei))//']',Level=7)
1125
    
1126
    
1127
    ! Compute the average
1128
    n = 0
1129
    DO i=1,RefMesh % NumberOfNodes
1130
      j = p(i)
1131
      IF (j<=0) CYCLE
1132
      IF(A % ParallelInfo % GInterface(j)) THEN
1133
        IF( Hcount(i) == 0 ) THEN
1134
          n = n+1
1135
        ELSE          
1136
          Hvalue(i) = Hvalue(i) / Hcount(i) 
1137
        END IF
1138
      END IF
1139
    END DO
1140
    
1141
    n = ParallelReduction(n)
1142
    CALL Info('ParallelAverageHvalue','Nodes '//TRIM(I2S(n))//' surrounded by orphans only!')
1143

1144
    ! Check dofs that have not been defined by averaging
1145
    IF( n > 0 ) THEN
1146
      Hcount = -Hcount
1147

1148
      DO i=1,RefMesh % NumberOfNodes
1149
        j = p(i)
1150
        IF(j<=0) CYCLE
1151
        IF(A % ParallelInfo % GInterface(j)) THEN
1152
          IF(Hcount(i) == 0) THEN
1153
            DO l=A % Rows(j),A % Rows(j+1)-1
1154
              k = ip( A % Cols(l) )
1155
              IF( k<=0 ) CYCLE
1156

1157
              IF(i==k) CYCLE
1158
              ! Use only nodes that were defined by averaging for the interface.
1159
              IF(Hcount(k) < 0) THEN
1160
                Hvalue(i) = Hvalue(i) + Hvalue(k)
1161
                Hcount(i) = Hcount(i) + 1
1162
              END IF
1163
            END DO
1164
          END IF
1165
        END IF
1166
      END DO
1167

1168
      ! This is a trick to get the already computed nodes to be properly re-everaged
1169
      WHERE(Hcount<0) Hcount = 1 
1170
      
1171
      ! Perform parallel summation, only interface gets summed. 
1172
      BLOCK
1173
        INTEGER, ALLOCATABLE :: Xcount(:)
1174
        REAL(KIND=dp), ALLOCATABLE :: Xvalue(:)
1175
        ALLOCATE( Xcount(A % NumberOfRows), XValue(A % NumberOfRows) )
1176
        Xcount = 0; Xvalue = 0
1177

1178
        DO i=1,RefMesh % NumberOfNodes
1179
          j = p(i)
1180
          IF (j<=0) CYCLE
1181
          Xvalue(j) = HValue(i)
1182
          Xcount(j) = HCount(i)
1183
        END DO
1184

1185
        CALL ParallelSumVector( A, Xvalue )
1186
        CALL ParallelSumVectorInt( A, Xcount ) 
1187

1188
        DO i=1,RefMesh % NumberOfNodes
1189
          j = p(i)
1190
          IF (j<=0) CYCLE
1191
          Hvalue(i) = Xvalue(j)
1192
          Hcount(i) = Xcount(j)
1193
        END DO
1194
      END BLOCK
1195

1196
!     CALL ParallelSumVector( A, Hvalue )
1197
!     CALL ParallelSumVectorInt( A, Hcount ) 
1198
      
1199
      n = 0
1200
      DO i=1,RefMesh % NumberOfNodes
1201
        j = p(i)
1202
        IF (j<=0) CYCLE 
1203
        IF(A % ParallelInfo % GInterface(j)) THEN
1204
          IF( Hcount(i) == 0 ) THEN
1205
            n = n+1
1206
          ELSE 
1207
            Hvalue(i) = Hvalue(i) / Hcount(i) 
1208
          END IF
1209
        END IF
1210
      END DO
1211
      
1212
      CALL Info('ParallelAverageHvalue','Nodes '//TRIM(I2S(n))//' surrounded by orphans only again!')        
1213
    END IF
1214
    
1215
!    IF( InfoActive(20) ) THEN
1216
!      CALL VectorValuesRange(Hvalue,SIZE(Hvalue),'Hvalue')             
1217
!    END IF
1218
    
1219
  END SUBROUTINE ParallelAverageHvalue
1220

1221

1222
  
1223
#ifdef HAVE_MMG
1224

1225
!------------------------------------------------------------------------------
1226
  FUNCTION MMG_ReMesh( RefMesh, ErrorLimit, HValue, NodalError, &
1227
       hConvergence, minH, maxH, MaxChange, Coarsening ) RESULT( NewMesh )
1228
!------------------------------------------------------------------------------
1229
    REAL(KIND=dp) :: NodalError(:), hConvergence(:), &
1230
           ErrorLimit, minH, maxH, MaxChange, HValue(:)
1231
    LOGICAL :: Coarsening
1232
    TYPE(Mesh_t), POINTER :: NewMesh, RefMesh
1233
!------------------------------------------------------------------------------
1234
    TYPE(Mesh_t), POINTER :: Mesh, TmpMesh, GatheredMesh
1235
    INTEGER :: i,j,k,n,ierr
1236
    REAL(KIND=dp) :: Lambda
1237
    CHARACTER(LEN=MAX_NAME_LEN) :: MeshCommand, Name
1238
    CHARACTER(LEN=MAX_PATH_LEN) :: MeshInputFile
1239
    LOGICAL :: Success, Rebalance, EnforceSerial
1240
    REAL(KIND=dp) :: xmax, xmin, ymax, ymin, zmax, zmin, cscale
1241
    LOGICAL :: ScaleCoord
1242
    INTEGER :: DoerPart
1243
    REAL(KIND=dp), POINTER :: NodalVals(:,:)
1244
    CHARACTER(*), PARAMETER :: Caller = 'MMG_ReMesh'
1245
        
1246
!------------------------------------------------------------------------------
1247

1248
#if 0
1249
    ! This is just some debugging code to check that the averaging works as intended.
1250
    TYPE(Variable_t), POINTER :: Var
1251
    REAL(KIND=dp), POINTER :: AveTest(:)
1252
    TYPE(Matrix_t), POINTER :: A
1253

1254
    n = RefMesh % NumberOfNodes 
1255
    ALLOCATE(AveTest(n))
1256
    AveTest(1:n) = RefMesh % Nodes % x(1:n) + & 
1257
        RefMesh % Nodes % y(1:n) + RefMesh % Nodes % z(1:n) 
1258

1259
    A => CurrentModel % Solver % Matrix       
1260
    WHERE( A % ParallelInfo % GInterface(1:n) ) 
1261
      AveTest(1:n) = -1000.0
1262
    END WHERE
1263
      
1264
    CALL VariableAdd( RefMesh % Variables, RefMesh, Solver, &
1265
        'Ave Test', 1, AveTest, Output = .TRUE.) 
1266
    CALL ParallelAverageHvalue( RefMesh, AveTest )
1267

1268
    AveTest(1:n) = AveTest(1:n) - ( RefMesh % Nodes % x(1:n) + & 
1269
        RefMesh % Nodes % y(1:n) + RefMesh % Nodes % z(1:n) )
1270
   
1271
    IF( InfoActive(20) ) THEN
1272
      CALL VectorValuesRange(Hvalue,SIZE(Hvalue),'Ave Test')             
1273
    END IF
1274
#endif
1275

1276
    CALL ComputeDesiredHvalue( RefMesh, ErrorLimit, Hvalue, NodalError, &
1277
        hConvergence, minH, maxH, MaxChange, Coarsening ) 
1278
    CALL ParallelAverageHvalue( RefMesh, Hvalue ) 
1279
    
1280
    Var => VariableGet( RefMesh % Variables, 'Hvalue', ThisOnly=.TRUE. )      
1281

1282
    IF( RefMesh % MeshDim == 2 ) THEN
1283
      IF( ParEnv % PEs > 1 ) THEN
1284
        CALL Fatal(Caller,'2D remeshing not available in parallel!')
1285
      ELSE
1286
        CALL Info(Caller,'Calling serial remeshing routines in 2D',Level=10)
1287
        NewMesh => MMG2D_ReMesh( RefMesh, Var )
1288
      END IF
1289
    ELSE
1290
      EnforceSerial = ListGetLogical( Params,'Adaptive Remesh Serial',Found )  
1291
      IF( ParEnv % PEs > 1 .AND. .NOT. EnforceSerial ) THEN
1292
        CALL Info(Caller,'Calling parallel remeshing routines in 3D',Level=10)
1293

1294
        ScaleCoord = ListGetLogical(Params, 'Adaptive Scale Coordinates for MMG', Found)
1295

1296
        IF(ScaleCoord) THEN
1297
          xmin = ParallelReduction( MINVAL(refmesh % nodes % x), 1)
1298
          xmax = ParallelReduction( MAXVAL(refmesh % nodes % x), 2)
1299

1300
          ymin = ParallelReduction( MINVAL(refmesh % nodes % y), 1)
1301
          ymax = ParallelReduction( MAXVAL(refmesh % nodes % y), 2)
1302

1303
          zmin = ParallelReduction( MINVAL(refmesh % nodes % z), 1)
1304
          zmax = ParallelReduction( MAXVAL(refmesh % nodes % z), 2)
1305

1306
          cscale = 1 * MAX( xmax - xmin, MAX( ymax - ymin, zmax - zmin)  )
1307

1308
          refmesh % nodes % x = cscale * refmesh % nodes % x
1309
          refmesh % nodes % y = cscale * refmesh % nodes % y
1310
          refmesh % nodes % z = cscale * refmesh % nodes % z
1311
          var % values = cscale * var % values
1312
        END IF
1313

1314
        CALL DistributedRemeshParMMG(Model, RefMesh, TmpMesh,&
1315
            Params = Solver % Values, HVar = Var )
1316

1317
        IF( ScaleCoord ) THEN
1318
          cscale = 1._dp / cscale
1319

1320
          tmpmesh % nodes % x = cscale * tmpmesh % nodes % x
1321
          tmpmesh % nodes % y = cscale * tmpmesh % nodes % y
1322
          tmpmesh % nodes % z = cscale * tmpmesh % nodes % z
1323

1324
          refmesh % nodes % x = cscale * refmesh % nodes % x
1325
          refmesh % nodes % y = cscale * refmesh % nodes % y
1326
          refmesh % nodes % z = cscale * refmesh % nodes % z
1327
          var % values = cscale * var % values
1328
        END IF
1329

1330
        CALL RenumberGElems(TmpMesh)
1331

1332
        Rebalance = ListGetLogical(Model % Solver % Values, "Adaptive Rebalance", Found, DefValue = .TRUE.)
1333
        IF(Rebalance) THEN
1334
          CALL Zoltan_Interface( Model, TmpMesh, StartImbalanceTol=1.1_dp, TolChange=0.02_dp, MinElems=10 )          
1335
          NewMesh => RedistributeMesh(Model, TmpMesh, .TRUE., .FALSE.)
1336
          CALL ReleaseMesh(TmpMesh)
1337
        ELSE
1338
          NewMesh => TmpMesh          
1339
        END IF
1340
      ELSE IF( ParEnv % PEs > 1 .AND. EnforceSerial ) THEN
1341
        CALL Info(Caller,'Calling serial remeshing routines for parallel 3D run',Level=10)
1342

1343
        n = RefMesh % NumberOfBulkElements + RefMesh % NumberOfBoundaryElements
1344
        IF(.NOT. ASSOCIATED(RefMesh % Repartition)) THEN
1345
          ALLOCATE(RefMesh % Repartition(n), STAT=ierr)
1346
          IF(ierr /= 0) CALL Fatal(Caller,'Could not ALLOCATE RefMesh % Repartition')
1347
        END IF               
1348
        ! For now, set to target with all to 1st partition
1349
        DoerPart = ListGetInteger( Params,'Adaptive Remesh Owner',Found ) 
1350

1351
        RefMesh % Repartition = DoerPart + 1
1352

1353
        PRINT *,'RefMesh: ',ParEnv % MyPe, &
1354
            RefMesh % NumberOfBulkElements, RefMesh % NumberOfBoundaryElements
1355

1356
        IF( ASSOCIATED( Var % Perm ) ) THEN
1357
          DO i=1,SIZE(Var % Perm)
1358
            IF(i /= Var % Perm(i)) THEN
1359
              CALL Fatal(Caller,'H field should have unity perm at this stage!')
1360
            END IF
1361
          END DO
1362
        END IF
1363
          
1364
        ALLOCATE(NodalVals(Refmesh % NumberOfNodes,1))
1365
        NodalVals(:,1) = Var % Values
1366
        GatheredMesh => RedistributeMesh(Model, RefMesh, .TRUE., .FALSE., NodalVals)
1367
!        GatheredMesh => RedistributeMesh(Model, RefMesh, .TRUE., .FALSE.)
1368
        
1369
        PRINT *,'GatheredMesh: ',ParEnv % MyPe, &
1370
            gatheredMesh % NumberOfBulkElements, gatheredMesh % NumberOfBoundaryElements
1371

1372
        IF( ParEnv % MyPe == DoerPart ) THEN
1373
#if 0 
1374
          ! Save the gathered serial mesh for debugging porposes
1375
          CALL WriteMeshToDisk2(Model, GatheredMesh,'gathered')
1376
#endif
1377

1378
          DEALLOCATE(Var % Values)
1379
          ALLOCATE(Var % Values(GatheredMesh % NumberOfNodes))
1380
          Var % Values = NodalVals(:,1)
1381
          IF(ASSOCIATED(Var % Perm)) DEALLOCATE(Var % Perm)
1382
          ALLOCATE(Var % Perm(GatheredMesh % NumberOfNodes))
1383
          DO i=1,GatheredMesh % NumberOfNodes
1384
            Var % Perm(i) = i
1385
          END DO
1386

1387
          ! Here do the adaptive remeshing with serial MMG3D since the parallel one is not very robust. 
1388
          CALL RemeshMMG3D(Model, GatheredMesh, TmpMesh,Params = Solver % Values, &
1389
              HVar = Var, Success = Success )
1390

1391
          ! Thereafter partition the mesh in a serial manner. 
1392
          IF( ListGetString( Solver % Values,'Partitioning method',Found) == 'zoltan') THEN
1393
            CALL Zoltan_Interface( Model, TmpMesh, SerialMode = .TRUE., NoPartitions = ParEnv % PEs, &
1394
                StartImbalanceTol=1.1_dp, TolChange=0.02_dp, MinElems=10 )
1395
          ELSE
1396
            CALL PartitionMeshSerial( Model, TmpMesh, Solver % Values )
1397
          END IF
1398

1399
#if 0 
1400
          ! Save the remeshed serial mesh for debugging porposes
1401
          CALL WriteMeshToDisk2(Model, TmpMesh,'remeshed')
1402
#endif
1403
        ELSE
1404
          TmpMesh => AllocateMesh()
1405
          TmpMesh % MeshDim = RefMesh % MeshDim
1406
        END IF
1407
        !CALL RenumberGElems(TmpMesh)
1408

1409
        NewMesh => RedistributeMesh(Model, TmpMesh, .TRUE., .FALSE.)
1410
        CALL ReleaseMesh(TmpMesh)
1411
      ELSE              
1412
        CALL Info(Caller,'Calling serial remeshing routines in 3D',Level=10)
1413
        CALL RemeshMMG3D(Model, RefMesh, NewMesh,Params = Solver % Values, &
1414
            HVar = Var, Success = Success )
1415
      END IF
1416
      CALL Info(Caller,'Finished MMG remeshing',Level=20)      
1417
    END IF
1418

1419
!------------------------------------------------------------------------------
1420
  END FUNCTION MMG_Remesh
1421
!------------------------------------------------------------------------------
1422
#endif
1423

1424
  
1425
!------------------------------------------------------------------------------
1426
  FUNCTION External_ReMesh( RefMesh, ErrorLimit, HValue, NodalError, &
1427
       hConvergence, minH, maxH, MaxChange, Coarsening ) RESULT( NewMesh )
1428
!------------------------------------------------------------------------------
1429
    REAL(KIND=dp) :: NodalError(:), hConvergence(:), &
1430
           ErrorLimit, minH, maxH, MaxChange, HValue(:)
1431
    LOGICAL :: Coarsening
1432
    TYPE(Mesh_t), POINTER :: NewMesh, RefMesh
1433
!------------------------------------------------------------------------------
1434
    TYPE(Mesh_t), POINTER :: Mesh
1435
    INTEGER :: i,j,k,n,dim
1436
    REAL(KIND=dp) :: Lambda
1437
    REAL(KIND=dp), POINTER :: HValueF(:)
1438
    CHARACTER(:), ALLOCATABLE :: MeshCommand, Name, MeshInputFile
1439
    LOGICAL :: GmshFormat
1440
    LOGICAL :: GmshPosFormat
1441
!------------------------------------------------------------------------------
1442

1443
    dim = CoordinateSystemDimension()
1444

1445
    ! Create a temporal field that includes the desired mesh density where Hvalue has been computed. 
1446
    ! This is in terms of desired mesh density, currently -1 value is given if the nodal error is zero
1447
    ! implying that nothing is computed here. 
1448
    ALLOCATE(HvalueF(SIZE(HValue)))
1449
    HValueF = -1.0_dp
1450
    DO i=1,RefMesh % NumberOfNodes
1451
      IF ( NodalError(i) > 100*AEPS ) THEN
1452
        Lambda = ( ErrorLimit / NodalError(i) ) ** ( 1.0d0 / hConvergence(i) )        
1453
        IF ( RefMesh % AdaptiveDepth < 1 ) THEN
1454
          Lambda = HValue(i) * MAX( MIN( Lambda, 1.33d0), 0.75d0)
1455
        ELSE
1456
          Lambda = HValue(i) * MAX(MIN(Lambda, MaxChange), 1.0d0/MaxChange)
1457
        END IF        
1458
        IF( .NOT.Coarsening ) Lambda = MIN( Lambda, Hvalue(i) )        
1459

1460
        IF ( maxH > 0 ) Lambda = MIN( Lambda, maxH )
1461
        IF ( minH > 0 ) Lambda = MAX( Lambda, minH )        
1462
        HValueF(i) = Lambda
1463
      END IF
1464
    END DO
1465

1466
    ! Save the current mesh in Elmer mesh format 
1467
    Path = ListGetString( Params, 'Adaptive Mesh Name', Found )
1468
    nlen = LEN_TRIM(Path)
1469
   
1470
    IF ( .NOT. Found ) THEN
1471
      i = RefMesh % AdaptiveDepth + 1
1472
      Path = 'RefinedMesh'//I2S(i)
1473
    END IF
1474

1475
    nlen = LEN_TRIM(OutputPath)
1476
    IF ( nlen > 0 ) THEN
1477
      Path = OutputPath(1:nlen) // '/' // TRIM(Path)
1478
    ELSE
1479
      Path = TRIM(Path)
1480
    END IF
1481
    
1482
    GmshFormat = ListGetLogical( Params,'Adaptive Remesh Use Gmsh', Found )
1483

1484
    IF( GmshFormat ) THEN
1485

1486
      GmshPosFormat = ListGetLogical( Params,'Adaptive Remesh Gmsh Use Pos Format', Found )
1487
      ! write the bacground mesh in .pos format if user requested it.
1488
      IF( GmshPosFormat) THEN
1489
      
1490
        ! Get the coordinate scaling. This is used to scale the background mesh coordinates according to the original mesh.
1491
        MeshDim = RefMesh % MaxDim
1492

1493
        Wrk => ListGetConstRealArray( Model % Simulation,'Coordinate Scaling',Found )
1494
        CoordScale = 1.0_dp    
1495
        IF( Found ) THEN            
1496
        DO i=1, MeshDim
1497
          j = MIN( i, SIZE(Wrk,1) )
1498
          CoordScale(i) = Wrk(j,1)
1499
        END DO
1500
        WRITE(Message,'(A,3ES10.3)') 'Scaling the background mesh coordinates:',CoordScale(1:3)
1501
        CALL Info(Caller ,Message, Level=10)
1502
        END IF 
1503
        
1504
        ! write the background mesh in .pos format
1505
        CALL Info( Caller,'Saving background mesh density in gmsh .pos format' )
1506
        OPEN( 11, STATUS='UNKNOWN',FILE='gmsh_bgmesh.pos' )
1507
        WRITE( 11,* ) 'View "mesh size field" {'           
1508
        DO i=1,RefMesh % NumberOfNodes
1509
          IF(.NOT. (HValueF(i) > 0.0_dp )) CYCLE
1510
          IF (dim == 2 ) THEN
1511
            WRITE( 11,* ) 'SP(', (RefMesh % Nodes % x(i)) / CoordScale(1), &
1512
                        ', ', (RefMesh % Nodes % y(i)) / CoordScale(2), ') {', &
1513
                        HValueF(i) / MIN(CoordScale(1), CoordScale(2)), '};'
1514
          ELSE
1515
            WRITE( 11,* ) 'SP(', (RefMesh % Nodes % x(i)) / CoordScale(1), &
1516
                        ', ', (RefMesh % Nodes % y(i)) / CoordScale(2), &
1517
                        ', ', (RefMesh % Nodes % z(i)) / CoordScale(3), ') {', &
1518
                        HValueF(i) / MIN(CoordScale(1), MIN(CoordScale(2), CoordScale(3))), '};'
1519
          END IF
1520
        END DO
1521
        WRITE( 11,* ) '};'
1522
        CLOSE(11)
1523
      ELSE
1524

1525
        CALL Info( Caller,'Saving background mesh density in gmsh 2.0 (.msh) format' )
1526

1527
        ! A cludge to change the pointer and save results in Gmsh format.
1528
        BLOCK
1529
          REAL(KIND=dp), POINTER :: PtoHvalue(:)
1530
          TYPE(Variable_t), POINTER :: HVar
1531
          HVar => VariableGet( RefMesh % Variables,'Hvalue')
1532
          pToHvalue => HVar % Values
1533
          HVar % Values => HvalueF        
1534
          CALL ListAddString(Solver % Values,'Scalar Field 1','Hvalue')
1535
          CALL ListAddLogical(Solver % Values,'File Append',.FALSE.)
1536
          CALL ListAddLogical(Solver % Values,'Alter Topology',.TRUE.) 
1537
          CALL ListAddNewString(Solver % Values, 'Output File Name', 'gmsh_bgmesh.msh')        
1538
          CALL SaveGmshOutput( Model,Solver,0.0_dp,.FALSE.)
1539
          HVar % Values => PtoHvalue
1540
        END BLOCK
1541
      END IF
1542
    ELSE      
1543
      CALL Info( Caller,'Saving background mesh density in point cloud format' )
1544

1545
      OPEN( 11, STATUS='UNKNOWN', FILE='bgmesh.nodes' )
1546
      WRITE( 11,* ) COUNT( HValueF > 0.0_dp )           
1547
      DO i=1,RefMesh % NumberOfNodes
1548
        IF(.NOT. (HValueF(i) > 0.0_dp )) CYCLE
1549
        IF (dim == 2 ) THEN
1550
          WRITE(11,'(3e23.15)') RefMesh % Nodes % x(i), &
1551
              RefMesh % Nodes % y(i), HValueF(i)
1552
        ELSE
1553
          WRITE(11,'(4e23.15)') RefMesh % Nodes % x(i), &
1554
              RefMesh % Nodes % y(i), &
1555
              RefMesh % Nodes % z(i), HValueF(i)
1556
        END IF
1557
      END DO
1558
      WRITE(11,*) 0
1559
      CLOSE(11)
1560

1561
      CALL MakeDirectory( TRIM(Path) // CHAR(0) )
1562
      CALL WriteMeshToDisk( RefMesh, Path )
1563
    END IF
1564

1565
    Mesh => RefMesh
1566
    DO WHILE( ASSOCIATED( Mesh ) )
1567
      IF ( Mesh % AdaptiveDepth == 0 ) EXIT
1568
      Mesh => Mesh % Parent
1569
    END DO
1570

1571
    MeshCommand = ListGetString( Solver % Values,'Mesh Command',Found) 
1572
    IF(.NOT. Found ) THEN
1573
      IF( GmshFormat ) THEN
1574
        CALL Fatal('ReMesh','For now, provide "Mesh Command" for Gmsh meshing!')
1575
      END IF
1576

1577
      MeshInputFile = ListGetString( Params, 'Mesh Input File', Found )
1578
      IF ( .NOT. Found ) THEN
1579
        MeshInputFile = ListGetString( Model % Simulation, 'Mesh Input File' )
1580
      END IF
1581

1582
      MeshCommand = TRIM(OutputPath) // '/' // TRIM(Mesh % Name) // '/' // &
1583
          TRIM( MeshInputFile )
1584

1585
      SELECT CASE( dim )
1586
      CASE(2)
1587
        ! Legacy mesh generator from the Elmer suite.
1588
        MeshCommand = 'Mesh2D '//TRIM(MeshCommand)//' '//TRIM(Path)// ' --bgmesh=bgmesh.nodes'
1589
      CASE(3)
1590
        MeshCommand = 'Mesh3D '//TRIM(MeshCommand)//' '//TRIM(Path)//' bgmesh.nodes'
1591
      END SELECT
1592
    END IF
1593

1594
    ! Remeshing command. 
1595
    CALL Info('ReMesh','Meshing command: '//TRIM(MeshCommand),Level=10)
1596
    CALL SystemCommand( MeshCommand )
1597

1598
    ! Check if also conversion command is given. 
1599
    MeshCommand = ListGetString( Solver % Values,'Mesh Conversion Command',Found)
1600
    IF( Found ) THEN
1601
      ! add the output path to the command. 
1602
      MeshCommand = MeshCommand // ' -out ' // TRIM(Path)
1603
      CALL Info('ReMesh','Conversion command: '//TRIM(MeshCommand),Level=10)
1604
      CALL SystemCommand( MeshCommand )
1605
    END IF    
1606

1607
    ! Read the new mesh. 
1608
    NewMesh => LoadMesh2( Model, OutPutPath, Path, .FALSE., 1, 0 )
1609

1610
    ! Loading Gebhart factors is more or less obsolete. 
1611
    IF ( Solver % Variable % Name == 'temperature' ) THEN
1612
       Name = ListGetString( Model % Simulation, 'Gebhart Factors', Found )
1613
       IF ( Found ) THEN
1614
          MeshCommand = 'View ' // TRIM(OutputPath) // &
1615
                '/' // TRIM(Mesh % Name) // ' ' // TRIM(Path)
1616
          CALL SystemCommand( MeshCommand )
1617

1618
          Name = TRIM(OutputPath) // '/' // &
1619
              TRIM(Mesh % Name) // '/' // TRIM(Name)
1620
          CALL LoadGebhartFactors( NewMesh, TRIM(Name) )
1621
       END IF
1622
    END IF
1623

1624
    DEALLOCATE(HvalueF)
1625
    
1626
!------------------------------------------------------------------------------
1627
  END FUNCTION External_ReMesh
1628
!------------------------------------------------------------------------------
1629

1630

1631
!------------------------------------------------------------------------------
1632
  FUNCTION SplitMesh( RefMesh,ErrorIndicator,ErrorLimit, NodalError, &
1633
       Hvalue, hConvergence, minH, maxH, MaxChange ) RESULT(NewMesh)
1634
!------------------------------------------------------------------------------
1635
    REAL(KIND=dp) :: NodalError(:), hConvergence(:), Hvalue(:), MaxChange
1636
    TYPE(Mesh_t), POINTER :: NewMesh, RefMesh
1637
    REAL(KIND=dp) :: ErrorIndicator(:),ErrorLimit,minH,maxH
1638
!------------------------------------------------------------------------------
1639
    TYPE(Mesh_t), POINTER :: NewMesh1
1640
    REAL(KIND=dp) :: Lambda, EhConvergence
1641
    INTEGER :: i,j,k,n,MarkedElements
1642
    TYPE(Element_t), POINTER :: RefElement
1643
!------------------------------------------------------------------------------
1644

1645
    NULLIFY( NewMesh )
1646

1647
!   Determine the marked elements:
1648
!   ------------------------------
1649
    MarkedElements = 0
1650

1651
    DO i = 1,RefMesh % NumberOfBulkElements
1652
       RefElement => RefMesh % Elements(i)
1653

1654
       IF ( RefElement % TYPE % ElementCode /= 303 ) THEN
1655
          CALL Fatal( 'SplitMesh', 'Internal splitting implemented only for linear triangles.' )
1656
       END IF
1657

1658
       n = RefElement % TYPE % NumberOfNodes
1659

1660
       IF( RefMesh % AdaptiveDepth < 1 ) THEN
1661
          EhConvergence = 1.0d0 ! First round: Assume full convergence speed
1662
       ELSE
1663
          EhConvergence = SUM( hConvergence( RefElement % Nodeindexes(1:n) ) ) / n
1664
       END IF
1665

1666
       RefElement % Splitted = 0
1667
       IF( ErrorIndicator(i) > 100*AEPS ) THEN
1668
          Lambda = ( ErrorLimit / ErrorIndicator(i) ) ** ( 1.0d0 / EhConvergence )
1669
          RefElement % Splitted = MIN( MaxChange, 1.0d0/Lambda )
1670
       END IF
1671

1672
       IF ( RefElement % Splitted > 0 ) MarkedElements = MarkedElements  + 1
1673
    END DO
1674

1675
    IF ( MarkedElements == 0 ) THEN
1676
       RefMesh % Changed = .FALSE.
1677
       RETURN
1678
    END IF
1679

1680
!   Refine until all elements splitted specified times:
1681
!   ---------------------------------------------------
1682
    NewMesh => SplitOneLevel( RefMesh )
1683
    DO WHILE( .TRUE. )
1684
       MarkedElements = 0
1685
       DO i=1,NewMesh % NumberOfBulkElements
1686
          IF ( NewMesh % Elements(i) % Splitted > 0 ) THEN
1687
             MarkedElements = MarkedElements + 1
1688
          END IF
1689
       END DO
1690

1691
       IF ( MarkedElements == 0 ) EXIT
1692

1693
       NewMesh1 => SplitOneLevel( NewMesh )
1694
       CALL ReleaseMesh( NewMesh )
1695
       DEALLOCATE( NewMesh )
1696

1697
       NewMesh => NewMesh1
1698
    END DO
1699

1700
!------------------------------------------------------------------------------
1701
  END FUNCTION SplitMesh
1702
!------------------------------------------------------------------------------
1703

1704

1705
!------------------------------------------------------------------------------
1706
  FUNCTION SplitOneLevel( RefMesh ) RESULT( NewMesh )
1707
!------------------------------------------------------------------------------
1708
    IMPLICIT NONE
1709

1710
    TYPE( Mesh_t ), POINTER :: RefMesh, NewMesh
1711
!------------------------------------------------------------------------------
1712
    REAL(KIND=dp) :: t
1713
    INTEGER :: EdgeNumber,LongestEdge,Node1,Node2
1714
    INTEGER :: i,j,k,l,n,NewElCnt,NewNodeCnt,MarkedEdges
1715

1716
    TYPE(Element_t), POINTER :: RefElement,Parent,Child,Edge
1717

1718
    LOGICAL, POINTER :: EdgeSplitted(:)
1719
    INTEGER, POINTER :: MarkedOrder(:), Children(:,:)
1720

1721
    TYPE(PElementDefs_t), POINTER :: PD
1722
    REAL(KIND=dp) :: x1, x2, y1, y2, EdgeLength, MaxLength
1723
!------------------------------------------------------------------------------
1724

1725
    t = CPUTime()
1726
    CALL FindMeshEdges( RefMesh )
1727
    WRITE( Message, * ) 'Find mesh edges time (cpu-secs):                 ',CPUTime()-t
1728
    CALL Info( 'SplitOneLevel', Message, Level=6 )
1729

1730
!   RGB Refinement:
1731
!   ---------------
1732
    t = CPUTime()
1733
    CALL AllocateVector( EdgeSplitted, RefMesh % NumberOfEdges )
1734
    MarkedEdges = RGBRefinement( EdgeSplitted,RefMesh )
1735
    WRITE( Message, * ) 'RGB Refinement time (cpu-secs):                  ',CPUTime()-t
1736
    CALL Info( 'SplitOneLevel', Message, Level=6 )
1737

1738
!   Initialize the new mesh:
1739
!   ------------------------
1740
    NewMesh => AllocateMesh()
1741
    NewMesh % MaxElementNodes = 3
1742
    NewMesh % MaxElementDOFs  = 3
1743
    NewMesh % MeshDim = RefMesh % MeshDim
1744

1745
!   Create node tables for the new mesh:
1746
!   ------------------------------------    
1747
    t = CPUTime()
1748
    NewMesh % NumberOfNodes = RefMesh % NumberOfNodes + MarkedEdges
1749
    CALL AllocateVector( NewMesh % Nodes % x, NewMesh % NumberOfNodes )
1750
    CALL AllocateVector( NewMesh % Nodes % y, NewMesh % NumberOfNodes )
1751
    CALL AllocateVector( NewMesh % Nodes % z, NewMesh % NumberOfNodes )
1752

1753
!   Add old nodes to the new mesh:
1754
!   ------------------------------    
1755
    NewMesh % Nodes % x(1:RefMesh % NumberOfNodes) = &
1756
               RefMesh % Nodes % x(1:RefMesh % NumberOfNodes)
1757
    NewMesh % Nodes % y(1:RefMesh % NumberOfNodes) = &
1758
               RefMesh % Nodes % y(1:RefMesh % NumberOfNodes)
1759
    NewMesh % Nodes % z(1:RefMesh % NumberOfNodes) = &
1760
               RefMesh % Nodes % z(1:RefMesh % NumberOfNodes)
1761

1762
!   Add new nodes to the new mesh:
1763
!   ------------------------------    
1764
    NewNodeCnt = RefMesh % NumberOfNodes
1765
    DO i = 1,RefMesh % NumberOfEdges
1766
       IF ( EdgeSplitted(i) ) THEN
1767
          Node1 = RefMesh % Edges(i) % NodeIndexes(1)
1768
          Node2 = RefMesh % Edges(i) % NodeIndexes(2)
1769
          x1 = RefMesh % Nodes % x(Node1)
1770
          x2 = RefMesh % Nodes % x(Node2)
1771
          y1 = RefMesh % Nodes % y(Node1)
1772
          y2 = RefMesh % Nodes % y(Node2)
1773
          NewNodeCnt = NewNodeCnt + 1
1774
          NewMesh % Nodes % x(NewNodeCnt) = (x1+x2) / 2.0d0
1775
          NewMesh % Nodes % y(NewNodeCnt) = (y1+y2) / 2.0d0
1776
          NewMesh % Nodes % z(NewNodeCnt) = 0.0d0
1777
       END IF
1778
    END DO
1779
    WRITE( Message, * ) 'Node tables generation time (cpu-secs):          ',CPUTime()-t
1780
    CALL Info( 'SplitOneLevel', Message, Level=6 )
1781

1782
!   Count the new number of bulk elements:
1783
!   --------------------------------------
1784
    CALL AllocateVector( MarkedOrder, RefMesh % NumberOfEdges )
1785
    MarkedOrder = 0
1786

1787
    k = 0
1788
    NewElCnt = 0
1789
    DO i = 1,RefMesh % NumberOfBulkElements
1790
       MarkedEdges = 0
1791
       DO j = 1,3
1792
          EdgeNumber = RefMesh % Elements(i) % EdgeIndexes(j)
1793
          IF( EdgeSplitted(EdgeNumber) ) THEN
1794
             MarkedEdges = MarkedEdges + 1
1795
             IF ( MarkedOrder(EdgeNumber) == 0 ) THEN
1796
                k = k + 1
1797
                MarkedOrder(EdgeNumber) = k + RefMesh % NumberOfNodes
1798
             END IF
1799
          END IF
1800
       END DO
1801
       NewElCnt = NewElCnt + MarkedEdges + 1
1802
    END DO
1803
    NewMesh % NumberOfBulkElements = NewElCnt
1804
!
1805
!   Count the new number of boundary elements:
1806
!   ------------------------------------------
1807
    NewElCnt = 0
1808
    DO i = RefMesh % NumberOfBulkElements+1,RefMesh % NumberOfBulkElements+&
1809
         RefMesh % NumberOfBoundaryElements
1810

1811
       RefElement => RefMesh % Elements(i) % BoundaryInfo % Left
1812
       IF ( .NOT.ASSOCIATED( RefElement) ) &
1813
            RefElement => RefMesh % Elements(i) % BoundaryInfo % Right
1814

1815
       IF ( ASSOCIATED( RefElement ) ) THEN
1816
          NULLIFY( Edge )
1817

1818
          DO j=1,3
1819
             Edge => RefMesh % Edges(RefElement % EdgeIndexes(j))
1820

1821
             IF ( Edge % NodeIndexes(1) == RefMesh % Elements(i) % NodeIndexes(1) .AND. &
1822
                  Edge % NodeIndexes(2) == RefMesh % Elements(i) % NodeIndexes(2) .OR.  &
1823
                  Edge % NodeIndexes(2) == RefMesh % Elements(i) % NodeIndexes(1) .AND. &
1824
                  Edge % NodeIndexes(1) == RefMesh % Elements(i) % NodeIndexes(2) ) EXIT
1825
          END DO
1826
   
1827
          IF ( EdgeSplitted( RefElement % EdgeIndexes(j) ) ) THEN
1828
             NewElCnt = NewElCnt + 2
1829
          ELSE
1830
             NewElCnt = NewElCnt + 1
1831
          END IF
1832
       ELSE
1833
          NewElCnt = NewElCnt + 1
1834
       END IF
1835
    END DO
1836

1837
    NewMesh % NumberOfBoundaryElements = NewElCnt
1838

1839
!   Allocate element tables:
1840
!   ------------------------
1841
    t = CPUTime()
1842
    CALL AllocateVector( NewMesh % Elements, NewMesh % NumberOfBulkElements + &
1843
         NewMesh % NumberOfBoundaryElements )
1844

1845
    CALL AllocateArray( Children, RefMesh % NumberOfBulkElements + &
1846
             RefMesh % NumberOfBoundaryElements, 4 )
1847
    Children = 0
1848

1849
!   Find the new bulk elements:
1850
!   ---------------------------
1851
    NewElCnt    = 0
1852
    DO i = 1,RefMesh % NumberOfBulkElements
1853
       RefElement => RefMesh % Elements(i)
1854
       n = RefElement % TYPE % NumberOfNodes
1855

1856
       MarkedEdges = 0
1857
       DO j = 1,3
1858
          EdgeNumber = RefElement % EdgeIndexes(j)
1859
          IF ( EdgeSplitted(EdgeNumber) ) THEN
1860
             MarkedEdges = MarkedEdges + 1
1861
          END IF
1862
       END DO
1863

1864
!      Make elements for the new mesh:
1865
!      --------------------------------
1866
       SELECT CASE(MarkedEdges)
1867
       CASE(0)
1868
!         Just copy of the old one:
1869
!         -------------------------
1870
          NewElCnt = NewElCnt + 1
1871
          NewMesh % Elements(NewElCnt) = RefElement
1872
          NewMesh % Elements(NewElCnt) % ElementIndex = NewElCnt
1873
          CALL AllocateVector( NewMesh % Elements(NewElCnt) % NodeIndexes,n )
1874
          NewMesh % Elements(NewElCnt) % NodeIndexes(1:n) = &
1875
               RefElement % NodeIndexes(1:n)
1876

1877
          Children(i,1) = NewElCnt
1878
          
1879
!-------------------------------------------------------------------------
1880
       CASE(1)
1881
!         Bisect the longest edge to give two triangles:
1882
!         ----------------------------------------------
1883
          DO j = 1,3
1884
             EdgeNumber = RefElement % EdgeIndexes(j)
1885
             IF ( EdgeSplitted( EdgeNumber ) ) EXIT
1886
          END DO
1887
            
1888
!         Find node (k) opposite to the splitted edge:
1889
!         --------------------------------------------
1890
          DO k = 1,3
1891
             IF ( RefElement % NodeIndexes(k) /= &
1892
                  RefMesh % Edges(EdgeNumber) % NodeIndexes(1) .AND. &
1893
                  RefElement % NodeIndexes(k) /= &
1894
                  RefMesh % Edges(EdgeNumber) % NodeIndexes(2) ) EXIT
1895
          END DO
1896

1897
!         New element 1
1898
!         -------------
1899
          NewElCnt = NewElCnt + 1
1900
          NewMesh % Elements(NewElCnt) = RefElement
1901
          NewMesh % Elements(NewElCnt) % ElementIndex = NewElCnt
1902
          CALL AllocateVector( NewMesh % Elements(NewElCnt) % NodeIndexes,n )
1903

1904
          NewMesh % Elements(NewElCnt) % NodeIndexes(1) = &
1905
               RefElement % NodeIndexes(k)
1906
          NewMesh % Elements(NewElCnt) % NodeIndexes(2) = &
1907
               RefMesh % Edges(EdgeNumber) % NodeIndexes(1)
1908

1909
          NewMesh % Elements(NewElCnt) % NodeIndexes(3) = &
1910
               MarkedOrder(RefElement % EdgeIndexes(j))
1911

1912
          Children(i,1) = NewElCnt
1913

1914
!         New element 2
1915
!----------------------------------------------------
1916
          NewElCnt = NewElCnt + 1
1917
          NewMesh % Elements(NewElCnt) = RefElement
1918
          NewMesh % Elements(NewElCnt) % ElementIndex = NewElCnt
1919
          CALL AllocateVector( NewMesh % Elements(NewElCnt) % NodeIndexes,n )
1920

1921
          NewMesh % Elements(NewElCnt) % NodeIndexes(1) = &
1922
               RefElement % NodeIndexes(k)
1923

1924
          NewMesh % Elements(NewElCnt) % NodeIndexes(2) = &
1925
               MarkedOrder(RefElement % EdgeIndexes(j))
1926

1927
          NewMesh % Elements(NewElCnt) % NodeIndexes(3) = &
1928
               RefMesh % Edges(EdgeNumber) % NodeIndexes(2)
1929

1930
          Children(i,2) = NewElCnt
1931

1932
!-------------------------------------------------------------------------
1933
       CASE(2)
1934
!         Bisect two of the edges to give three new elements:
1935
!         ---------------------------------------------------
1936

1937
!         Find the edge NOT splitted:
1938
!         ---------------------------
1939
          DO j = 1,3
1940
             EdgeNumber = RefElement % EdgeIndexes(j)
1941
             IF ( .NOT.EdgeSplitted( EdgeNumber ) ) EXIT             
1942
          END DO
1943

1944
!         Find node (k) opposite to the edge NOT splitted:
1945
!         ------------------------------------------------
1946
          DO k = 1,3
1947
             IF (RefElement % NodeIndexes(k) /= &
1948
                  RefMesh % Edges(EdgeNumber) % NodeIndexes(1) .AND. &
1949
                  RefElement % NodeIndexes(k) /= &
1950
                  RefMesh % Edges(EdgeNumber) % NodeIndexes(2) ) EXIT
1951
          END DO
1952

1953
!         New element 1
1954
!----------------------------------------------------
1955
          NewElCnt = NewElCnt + 1
1956
          NewMesh % Elements(NewElCnt) = RefElement
1957
          NewMesh % Elements(NewElCnt) % ElementIndex = NewElCnt
1958
          CALL AllocateVector( NewMesh % Elements(NewElCnt) % NodeIndexes,n )
1959

1960
          NewMesh % Elements(NewElCnt) % NodeIndexes(1) = &
1961
               RefElement % NodeIndexes(k)
1962

1963
          l = 1
1964
          DO k = 1,3
1965
             IF ( k /= j ) THEN
1966
                l = l + 1
1967
                NewMesh % Elements(NewElCnt) % NodeIndexes(l) = &
1968
                     MarkedOrder(RefElement % EdgeIndexes(k))
1969
             END IF
1970
          END DO
1971

1972
          Children(i,1) = NewElCnt
1973

1974
!         New element 2
1975
!----------------------------------------------------
1976
          NewElCnt = NewElCnt + 1
1977
          NewMesh % Elements(NewElCnt) = RefElement
1978
          NewMesh % Elements(NewElCnt) % ElementIndex = NewElCnt
1979
          CALL AllocateVector( NewMesh % Elements(NewElCnt) % NodeIndexes,n )
1980

1981
          l = 0
1982
          DO k = 1,3
1983
             IF ( k /= j ) THEN
1984
                l = l + 1
1985
                NewMesh % Elements(NewElCnt) % NodeIndexes(l) = &
1986
                     MarkedOrder(RefElement % EdgeIndexes(k))
1987
             END IF
1988
          END DO
1989

1990
          MaxLength = 0.0d0
1991
          DO k = 1,3
1992
             IF ( k /= j ) THEN
1993
                EdgeNumber = RefElement % EdgeIndexes(k)
1994
                Node1 = RefMesh % Edges( EdgeNumber ) % NodeIndexes(1)
1995
                Node2 = RefMesh % Edges( EdgeNumber ) % NodeIndexes(2)
1996
                x1 = RefMesh % Nodes % x( Node1 )
1997
                x2 = RefMesh % Nodes % x( Node2 )
1998
                y1 = RefMesh % Nodes % y( Node1 )
1999
                y2 = RefMesh % Nodes % y( Node2 )
2000
                EdgeLength = SQRT((x2-x1)**2+(y2-y1)**2)
2001
                IF (EdgeLength >= MaxLength) THEN
2002
                   MaxLength = EdgeLength
2003
                   LongestEdge = k
2004
                END IF
2005
             END IF
2006
          END DO
2007
          k = LongestEdge
2008
          IF ( k <= 0 .OR. k > 3 ) PRINT *,'longest edge:',k
2009

2010
          IF ( RefMesh % Edges(RefElement % EdgeIndexes(j)) % NodeIndexes(1) ==  &
2011
               RefMesh % Edges(RefElement % EdgeIndexes(k)) % NodeIndexes(1) .OR.&
2012
               RefMesh % Edges(RefElement % EdgeIndexes(j)) % NodeIndexes(1) ==  &
2013
               RefMesh % Edges(RefElement % EdgeIndexes(k)) % NodeIndexes(2) ) THEN
2014
             NewMesh % Elements(NewElCnt) % NodeIndexes(3) = &
2015
                  RefMesh % Edges(RefElement % EdgeIndexes(j)) % NodeIndexes(2)
2016
          ELSE
2017
             NewMesh % Elements(NewElCnt) % NodeIndexes(3) = &
2018
                  RefMesh % Edges(RefElement % EdgeIndexes(j)) % NodeIndexes(1)
2019
          END IF
2020

2021
          Children(i,2) = NewElCnt
2022

2023
!         New element 3
2024
!----------------------------------------------------
2025
          NewElCnt = NewElCnt + 1
2026
          NewMesh % Elements(NewElCnt) = RefElement
2027
          NewMesh % Elements(NewElCnt) % ElementIndex = NewElCnt
2028
          CALL AllocateVector( NewMesh % Elements(NewElCnt) % NodeIndexes,n )
2029

2030
          DO j = 1,3
2031
             EdgeNumber = RefElement % EdgeIndexes(j)
2032
             IF ( .NOT.EdgeSplitted( EdgeNumber ) ) EXIT             
2033
          END DO
2034

2035
          DO k = 1,2
2036
             NewMesh % Elements(NewElCnt) % NodeIndexes(k) = &
2037
                  RefMesh % Edges(EdgeNumber) % NodeIndexes(k)
2038
          END DO
2039

2040
          NewMesh % Elements(NewElCnt) % NodeIndexes(3) = &
2041
               MarkedOrder(RefElement % EdgeIndexes(LongestEdge))
2042

2043
          Children(i,3) = NewElCnt
2044

2045
!-------------------------------------------------------------------------
2046
       CASE(3)
2047
!         Bisect all the edges to give four new elements:
2048
!         -----------------------------------------------
2049

2050
!         New element 1
2051
!----------------------------------------------------
2052
          NewElCnt = NewElCnt + 1
2053
          NewMesh % Elements(NewElCnt) = RefElement
2054
          NewMesh % Elements(NewElCnt) % ElementIndex = NewElCnt
2055
          CALL AllocateVector( NewMesh % Elements(NewElCnt) % NodeIndexes,n )
2056

2057
          NewMesh % Elements(NewElCnt) % NodeIndexes(1) = &
2058
               RefElement % NodeIndexes(1)
2059

2060
          j = RefElement % EdgeIndexes(1)
2061
          NewMesh % Elements(NewElCnt) % NodeIndexes(2) = MarkedOrder(j)
2062

2063
          j = RefElement % EdgeIndexes(3)
2064
          NewMesh % Elements(NewElCnt) % NodeIndexes(3) = MarkedOrder(j)
2065

2066
          Children(i,1) = NewElCnt
2067

2068
!         New element 2
2069
!----------------------------------------------------
2070
          NewElCnt = NewElCnt + 1
2071
          NewMesh % Elements(NewElCnt) = RefElement
2072
          NewMesh % Elements(NewElCnt) % ElementIndex = NewElCnt
2073
          CALL AllocateVector( NewMesh % Elements(NewElCnt) % NodeIndexes,n )
2074

2075
          NewMesh % Elements(NewElCnt) % NodeIndexes(1) = &
2076
               RefElement % NodeIndexes(2)
2077

2078
          j = RefElement % EdgeIndexes(2)
2079
          NewMesh % Elements(NewElCnt) % NodeIndexes(2) = MarkedOrder(j)
2080

2081
          j = RefElement % EdgeIndexes(1)
2082
          NewMesh % Elements(NewElCnt) % NodeIndexes(3) = MarkedOrder(j)
2083

2084
          Children(i,2) = NewElCnt
2085

2086
!         New element 3
2087
!----------------------------------------------------
2088
          NewElCnt = NewElCnt + 1
2089
          NewMesh % Elements(NewElCnt) = RefElement
2090
          NewMesh % Elements(NewElCnt) % ElementIndex = NewElCnt
2091
          CALL AllocateVector( NewMesh % Elements(NewElCnt) % NodeIndexes,n )
2092

2093
          NewMesh % Elements(NewElCnt) % NodeIndexes(1) = &
2094
               RefElement % NodeIndexes(3)
2095

2096
          j = RefElement % EdgeIndexes(3)
2097
          NewMesh % Elements(NewElCnt) % NodeIndexes(2) = MarkedOrder(j)
2098

2099
          j = RefElement % EdgeIndexes(2)
2100
          NewMesh % Elements(NewElCnt) % NodeIndexes(3) = MarkedOrder(j)
2101

2102
          Children(i,3) = NewElCnt
2103

2104
!         New element 4
2105
!----------------------------------------------------
2106
          NewElCnt = NewElCnt + 1
2107
          NewMesh % Elements(NewElCnt) = RefElement
2108
          NewMesh % Elements(NewElCnt) % ElementIndex = NewElCnt
2109
          CALL AllocateVector( NewMesh % Elements(NewElCnt) % NodeIndexes,n )
2110

2111
          DO j=1,n
2112
             NewMesh % Elements(NewElCnt) % NodeIndexes(j) = &
2113
                  MarkedOrder( RefElement % EdgeIndexes(j) )
2114
          END DO
2115

2116
          Children(i,4) = NewElCnt
2117
!----------------------------------------------------
2118
       END SELECT
2119

2120
!----------------------------------------------------
2121
       DO j=1,4
2122
          k = Children(i,j)
2123
          IF ( k > 0 ) THEN
2124
             NewMesh % Elements(k) % Splitted = RefElement % Splitted-1
2125
          END IF
2126
       END DO
2127
    END DO
2128

2129

2130
    WRITE( Message, * ) 'Bulk element tables generation time (cpu-secs):  ',CPUTime()-t
2131
    CALL Info( 'SplitOneLevel', Message, Level=6 )
2132
    
2133
!
2134
!   Update boundary elements:
2135
!   -------------------------
2136
    t = CPUTime()
2137
    NewElCnt = NewMesh % NumberOfBulkElements
2138
    DO j = RefMesh % NumberOfBulkElements + 1, &
2139
       RefMesh % NumberOfBulkElements + &
2140
          RefMesh % NumberOfBoundaryElements
2141

2142
       RefElement => RefMesh % Elements(j) % BoundaryInfo % Left
2143
       IF ( .NOT.ASSOCIATED( RefElement) ) &
2144
            RefElement => RefMesh % Elements(j) % BoundaryInfo % Right
2145

2146
       IF ( ASSOCIATED( RefElement ) ) THEN
2147
          NULLIFY( Edge )
2148
          DO i=1,3
2149
             Edge => RefMesh % Edges(RefElement % EdgeIndexes(i))
2150
             IF ( Edge % NodeIndexes(1) == RefMesh % Elements(j) % NodeIndexes(1) .AND. &
2151
                  Edge % NodeIndexes(2) == RefMesh % Elements(j) % NodeIndexes(2) .OR.  &
2152
                  Edge % NodeIndexes(2) == RefMesh % Elements(j) % NodeIndexes(1) .AND. &
2153
                  Edge % NodeIndexes(1) == RefMesh % Elements(j) % NodeIndexes(2) ) EXIT
2154
          END DO
2155
          EdgeNumber = RefElement % EdgeIndexes(i)
2156

2157
          RefElement => RefMesh % Elements(j)
2158
          n = RefElement % TYPE % NumberOfNodes
2159
            
2160
          IF ( EdgeSplitted(EdgeNumber) ) THEN
2161
!
2162
!            New element 1:
2163
!            --------------
2164
             NewElCnt = NewElCnt + 1
2165
             NewMesh % Elements(NewElCnt) = RefElement
2166
             NewMesh % Elements(NewElCnt) % ElementIndex = NewElCnt
2167
             CALL AllocateVector( NewMesh % Elements(NewElCnt) % NodeIndexes,n )
2168
             NewMesh % Elements(NewElCnt) % NodeIndexes(1) = &
2169
                  RefElement % NodeIndexes(1)
2170
             NewMesh % Elements(NewElCnt) % NodeIndexes(2) = &
2171
                  MarkedOrder(EdgeNumber)
2172

2173
             ALLOCATE( NewMesh % Elements(NewElCnt) % BoundaryInfo )
2174
             NewMesh % Elements(NewElCnt) % BoundaryInfo = &
2175
                  RefElement % BoundaryInfo
2176

2177
             NULLIFY( NewMesh % Elements(NewElCnt) % &
2178
                Boundaryinfo % RadiationFactors )
2179
               
2180
             CALL SetParents( NewMesh % Elements(NewElCnt), &
2181
                  NewMesh, Children, Edge )
2182

2183
             Children(j,1) = NewElCnt
2184
               
2185
!
2186
!            New element 2:
2187
!            --------------
2188
             NewElCnt = NewElCnt + 1
2189
             NewMesh % Elements(NewElCnt) = RefElement
2190
             NewMesh % Elements(NewElCnt) % ElementIndex = NewElCnt
2191
             CALL AllocateVector( NewMesh % Elements(NewElCnt) % NodeIndexes,n )
2192
             NewMesh % Elements(NewElCnt) % NodeIndexes(1) = &
2193
                  MarkedOrder(EdgeNumber)
2194
             NewMesh % Elements(NewElCnt) % NodeIndexes(2) = &
2195
                  RefElement % NodeIndexes(2)
2196

2197
             ALLOCATE( NewMesh % Elements(NewElCnt) % BoundaryInfo )
2198
             NewMesh % Elements(NewElCnt) % BoundaryInfo = &
2199
                  RefElement % BoundaryInfo
2200

2201
             NULLIFY( NewMesh % Elements(NewElCnt) % &
2202
                Boundaryinfo % RadiationFactors )
2203

2204
             CALL SetParents( NewMesh % Elements(NewElCnt), &
2205
                  NewMesh, Children, Edge )
2206

2207
             Children(j,2) = NewElCnt
2208
          ELSE
2209
!
2210
!            New element 1:
2211
!            --------------
2212
             NewElCnt = NewElCnt + 1
2213
             NewMesh % Elements(NewElCnt) = RefElement
2214
             NewMesh % Elements(NewElCnt) % ElementIndex = NewElCnt
2215
             CALL AllocateVector( NewMesh % Elements(NewElCnt) % NodeIndexes,n )
2216
             NewMesh % Elements(NewElCnt) % NodeIndexes = &
2217
                  RefElement % NodeIndexes
2218

2219
             ALLOCATE( NewMesh % Elements(NewElCnt) % BoundaryInfo )
2220

2221
             NewMesh % Elements(NewElCnt) % BoundaryInfo = &
2222
                  RefElement % BoundaryInfo
2223

2224
             NULLIFY( NewMesh % Elements(NewElCnt) % &
2225
                Boundaryinfo % RadiationFactors )
2226
            
2227
             CALL SetParents( NewMesh % Elements(NewElCnt), &
2228
                  NewMesh, Children, Edge )
2229

2230
             Children(j,1) = NewElCnt
2231
          END IF
2232
       ELSE
2233
!
2234
!         New element 1, this is point element:
2235
!         -------------------------------------
2236
          NewElCnt = NewElCnt + 1
2237
          RefElement => RefMesh % Elements(j)
2238
          n = RefElement % TYPE % NumberOfNodes
2239

2240
          NewMesh % Elements(NewElCnt) = RefElement
2241
          NewMesh % Elements(NewElCnt) % ElementIndex = NewElCnt
2242
          CALL AllocateVector( NewMesh % Elements(NewElCnt) % NodeIndexes,n )
2243
          NewMesh % Elements(NewElCnt) % NodeIndexes = &
2244
               RefElement % NodeIndexes
2245
               
2246
          ALLOCATE( NewMesh % Elements(NewElCnt) % BoundaryInfo )
2247

2248
          NewMesh % Elements(NewElCnt) % BoundaryInfo = &
2249
               RefElement % BoundaryInfo
2250
 
2251
          NULLIFY( NewMesh % Elements(NewElCnt) % &
2252
             Boundaryinfo % RadiationFactors )
2253

2254
          NULLIFY( NewMesh % Elements(NewElCnt) % BoundaryInfo % Left )
2255
          NULLIFY( NewMesh % Elements(NewElCnt) % BoundaryInfo % Right )
2256

2257
          Children(j,1) = NewElCnt
2258
       END IF
2259
    END DO
2260

2261
    NewMesh % MaxBDOFs = RefMesh % MaxBDOFs
2262
    DO i = 1,NewMesh % NumberOfBulkElements+NewMesh % NumberOfBoundaryElements
2263
      RefElement => NewMesh % Elements(i)
2264
      NULLIFY( RefElement % PDefs )
2265
      NULLIFY( RefElement % DGIndexes )
2266
      NULLIFY( RefElement % EdgeIndexes )
2267
      NULLIFY( RefElement % FaceIndexes )
2268
      NULLIFY( RefElement % BubbleIndexes )
2269
      IF ( RefElement % BDOFs > 0 ) THEN
2270
        ALLOCATE( RefElement % BubbleIndexes(RefElement % BDOFs) )
2271
        DO j=1,RefElement % BDOFs
2272
          RefElement % BubbleIndexes(j) = NewMesh % MaxBDOFs*(i-1)+j
2273
        END DO
2274
      END IF
2275

2276
      IF ( ASSOCIATED(RefElement % PDefs) ) THEN
2277
        PD => RefElement % PDefs
2278
        CALL AllocatePDefinitions(RefElement)
2279
        RefElement % PDefs = PD
2280
      END IF
2281
    END DO
2282

2283
!
2284
!   Update Gebhart factors, if present and the current solver
2285
!   is a heat equation solver:
2286
!   ------------------------------------------------------------
2287
    IF ( ListGetLogical( Solver % Values, 'Radiation Solver', Found ) ) THEN
2288
      CALL UpdateGebhartFactors( RefMesh, NewMesh, Children )
2289
    END IF
2290
      
2291
    WRITE( Message, * ) 'Bndry element tables generation time (cpu-secs): ',CPUTime()-t
2292
    CALL Info( 'SplitOneLevel', Message, Level=6 )
2293

2294
    DEALLOCATE( EdgeSplitted, MarkedOrder, Children )
2295
!------------------------------------------------------------------------------
2296
  END FUNCTION SplitOneLevel
2297
!------------------------------------------------------------------------------
2298

2299

2300
!------------------------------------------------------------------------------
2301
  FUNCTION RGBRefinement(  EdgeSplitted,RefMesh ) RESULT(MarkedEdges)
2302
!------------------------------------------------------------------------------
2303
    IMPLICIT NONE
2304

2305
    LOGICAL :: EdgeSplitted(:)
2306
    INTEGER :: MarkedEdges
2307
    TYPE(Mesh_t), POINTER :: RefMesh
2308
!------------------------------------------------------------------------------
2309
    LOGICAL :: MarkedEdgesFound
2310
    INTEGER :: i,j,EdgeNumber,HangingNodes,RGBIterations,Node1,Node2,&
2311
         LongestEdge
2312
    REAL(KIND=dp) :: x1,y1,x2,y2,EdgeLength,MaxLength
2313
!------------------------------------------------------------------------------
2314
    EdgeSplitted = .FALSE.
2315

2316
!   Mark all three edges of the marked elements (RED refinement):
2317
!   -------------------------------------------------------------
2318
!     DO i = 1,RefMesh % NumberOfBulkElements
2319
!        IF ( RefMesh % Elements(i) % Splitted > 0 ) THEN
2320
!           DO j = 1,3
2321
!              EdgeNumber = RefMesh % Elements(i) % EdgeIndexes(j)
2322
!              EdgeSplitted( EdgeNumber ) = .TRUE.
2323
!           END DO
2324
!        END IF
2325
!     END DO
2326

2327
!   Mark the longest edges of the marked elements (GREEN refinement):
2328
!   -----------------------------------------------------------------
2329
    DO i = 1,RefMesh % NumberOfBulkElements
2330
       IF ( RefMesh % Elements(i) % Splitted > 0 ) THEN
2331
          MaxLength   = 0.0D0
2332
          LongestEdge = 0
2333
          DO j = 1,3
2334
             EdgeNumber = RefMesh % Elements(i) % EdgeIndexes(j)
2335
             Node1 = RefMesh % Edges( EdgeNumber ) % NodeIndexes(1)
2336
             Node2 = RefMesh % Edges( EdgeNumber ) % NodeIndexes(2)
2337
             x1 = RefMesh % Nodes % x( Node1 )
2338
             x2 = RefMesh % Nodes % x( Node2 )
2339
             y1 = RefMesh % Nodes % y( Node1 )
2340
             y2 = RefMesh % Nodes % y( Node2 )
2341
             EdgeLength = SQRT((x2-x1)**2+(y2-y1)**2)
2342
             IF (EdgeLength >= MaxLength) THEN
2343
                MaxLength = EdgeLength
2344
                LongestEdge = EdgeNumber
2345
             END IF
2346
          END DO
2347
          EdgeSplitted( LongestEdge ) = .TRUE.
2348
       END IF
2349
    END DO
2350

2351
    MarkedEdges = 0
2352
    DO i = 1,RefMesh % NumberOfEdges
2353
       IF ( EdgeSplitted(i) ) THEN
2354
          MarkedEdges = MarkedEdges + 1
2355
       END IF
2356
    END DO
2357

2358
!   Mark longest edges until we have a RGB-refinement:
2359
!   --------------------------------------------------
2360
    RGBiterations = 0
2361
    DO WHILE( .TRUE. )
2362
       HangingNodes = 0
2363
       RGBiterations = RGBiterations+1
2364
       DO i = 1,RefMesh % NumberOfBulkElements
2365
            
2366
!         Check for marked edges and find the longest edge:
2367
!         -------------------------------------------------
2368
          MarkedEdgesFound = .FALSE.
2369
          LongestEdge      = 0
2370
          MaxLength        = 0.0d0
2371
          DO j = 1,3
2372
             EdgeNumber = RefMesh % Elements(i) % EdgeIndexes(j)
2373
             MarkedEdgesFound = MarkedEdgesFound.OR.EdgeSplitted(EdgeNumber)
2374
             Node1 = RefMesh % Edges(EdgeNumber) % NodeIndexes(1)
2375
             Node2 = RefMesh % Edges(EdgeNumber) % NodeIndexes(2)
2376
             x1 = RefMesh % Nodes % x( Node1 )
2377
             x2 = RefMesh % Nodes % x( Node2 )
2378
             y1 = RefMesh % Nodes % y( Node1 )
2379
             y2 = RefMesh % Nodes % y( Node2 )
2380
             EdgeLength = SQRT((x2-x1)**2+(y2-y1)**2)
2381
             IF (EdgeLength >= MaxLength) THEN
2382
                MaxLength = EdgeLength
2383
                LongestEdge = EdgeNumber
2384
             END IF
2385
          END DO
2386
          
2387
!         If there are marked edges, the longest edge must be one of them:
2388
!         ----------------------------------------------------------------
2389
          IF ( MarkedEdgesFound.AND.(.NOT.EdgeSplitted(LongestEdge)) ) THEN
2390
             HangingNodes = HangingNodes + 1
2391
             EdgeSplitted( LongestEdge ) = .TRUE.
2392
          END IF
2393
       END DO
2394

2395
       IF( HangingNodes > 0) THEN
2396
          WRITE( Message, * ) 'RGB ',RGBiterations,' : ',HangingNodes,' new nodes'
2397
          CALL Info( 'RGBRefinement', Message, Level=6 )
2398
          MarkedEdges = MarkedEdges + HangingNodes
2399
       ELSE
2400
          EXIT
2401
       END IF
2402
    END DO
2403
!------------------------------------------------------------------------------
2404
  END FUNCTION RGBRefinement
2405
!------------------------------------------------------------------------------
2406

2407

2408
!------------------------------------------------------------------------------
2409
! Find the parent elements to the splitted boundary element
2410
! among the children of the original parent element:
2411
! ---------------------------------------------------------
2412
!------------------------------------------------------------------------------
2413
  SUBROUTINE SetParents( Element, Mesh, Children, Edge )
2414
!------------------------------------------------------------------------------
2415
    TYPE(Element_t) :: Element
2416
    TYPE(Element_t), POINTER :: Edge
2417

2418
    INTEGER :: Children(:,:)
2419
    TYPE(Mesh_t), POINTER :: Mesh
2420

2421
    INTEGER j,k,l,n,i0,j0,k0
2422

2423
    TYPE(Element_t), POINTER :: Child
2424

2425
    n = Element % TYPE % NumberOfNodes
2426

2427
    k = Edge % BoundaryInfo % Left % ElementIndex
2428
    NULLIFY( Child )
2429
    DO l=1,4
2430
       IF ( Children(k,l)>0 ) THEN
2431
          Child => Mesh % Elements( Children(k,l) )
2432
          i0 = 0
2433
          DO j0=1,n
2434
             DO k0=1,Child % TYPE % NumberOfNodes
2435
                IF ( Child % NodeIndexes(k0) == Element % NodeIndexes(j0) ) THEN
2436
                   i0 = i0 + 1 
2437
                   EXIT
2438
                END IF
2439
             END DO
2440
          END DO
2441
          IF ( i0 == n ) EXIT
2442
       END IF
2443
    END DO
2444

2445
    IF ( l > 4 ) STOP 'Adaptive: parent 1 not found'
2446
        
2447
    Element % BoundaryInfo % Left  => Child
2448
    NULLIFY( Element % BoundaryInfo % Right )
2449
        
2450
    NULLIFY( Child )
2451
    IF ( ASSOCIATED(Edge % BoundaryInfo % Right) ) THEN
2452
       k = Edge % BoundaryInfo % Right % ElementIndex
2453
       DO l=1,4
2454
          IF ( Children(k,l)>0 ) THEN
2455
             Child => Mesh % Elements( Children(k,l) )
2456
             i0 = 0
2457
             DO j0=1,n
2458
                DO k0=1,Child % TYPE % NumberOfNodes
2459
                   IF ( Child % NodeIndexes(k0) == Element % NodeIndexes(j0) ) THEN
2460
                      i0 = i0 + 1 
2461
                      EXIT
2462
                   END IF
2463
                END DO
2464
             END DO
2465
             IF ( i0 == n ) EXIT
2466
          END IF
2467
       END DO
2468
           
2469
       Element % BoundaryInfo % Right => Child
2470
    END IF
2471
!------------------------------------------------------------------------------
2472
  END SUBROUTINE SetParents
2473
!------------------------------------------------------------------------------
2474

2475

2476
!------------------------------------------------------------------------------
2477
  SUBROUTINE UpdateGebhartFactors( RefMesh,NewMesh,Children ) 
2478
!------------------------------------------------------------------------------
2479
    TYPE(Mesh_t), POINTER :: RefMesh,NewMesh
2480
    INTEGER :: Children(:,:)
2481
!------------------------------------------------------------------------------
2482
    INTEGER :: i,j,k,n,NewFactors,TARGET
2483
    REAL(KIND=dp) :: AreaParent,AreaChild
2484
    TYPE(Factors_t), POINTER :: Factors,ChildFactors
2485
!------------------------------------------------------------------------------
2486
!
2487
!   Count numbers of factors for the new boundary elements:
2488
!   -------------------------------------------------------
2489
    DO i=RefMesh % NumberOfBulkElements+1,RefMesh % NumberOfBulkElements + &
2490
         RefMesh % NumberOfBoundaryElements
2491

2492
       Factors => RefMesh % Elements(i) % BoundaryInfo % RadiationFactors
2493
       IF ( .NOT. ASSOCIATED( Factors ) ) CYCLE
2494

2495
       NewFactors = 0
2496
       DO k=1,Factors % NumberOfFactors
2497
          TARGET = Factors % Elements(k)
2498
          IF ( Children(TARGET,2) > 0 ) THEN
2499
             NewFactors = NewFactors + 2
2500
          ELSE
2501
             NewFactors = NewFactors + 1
2502
          END IF
2503
       END DO
2504

2505
       IF (.NOT.ASSOCIATED(NewMesh % Elements(Children(i,1)) % &
2506
                BoundaryInfo % RadiationFactors) ) &
2507
         ALLOCATE(NewMesh % Elements(Children(i,1)) % BoundaryInfo % RadiationFactors)
2508

2509
       NewMesh % Elements(Children(i,1)) % BoundaryInfo % &
2510
            RadiationFactors % NumberOfFactors = NewFactors
2511

2512
       IF ( Children(i,2) > 0 ) THEN
2513
          IF (.NOT.ASSOCIATED(NewMesh % Elements(Children(i,2)) % &
2514
                BoundaryInfo % RadiationFactors) ) &
2515
            ALLOCATE(NewMesh % Elements(Children(i,2)) % BoundaryInfo % RadiationFactors)
2516

2517
          NewMesh % Elements(Children(i,2)) % BoundaryInfo % &
2518
               RadiationFactors % NumberOfFactors = NewFactors
2519
       END IF
2520
    END DO
2521

2522
!
2523
!   Update the factors:
2524
!   --------------------
2525
    DO i=RefMesh % NumberOfBulkElements+1,RefMesh % NumberOfBulkElements + &
2526
         RefMesh % NumberOfBoundaryElements
2527

2528
       Factors => RefMesh % Elements(i) % BoundaryInfo % RadiationFactors
2529
       IF ( .NOT. ASSOCIATED( Factors ) ) CYCLE
2530

2531
       AreaParent = ElementArea( RefMesh, RefMesh % Elements(i), &
2532
            RefMesh % Elements(i) % TYPE % NumberOfNodes )
2533

2534
       n = Children(i,1)
2535

2536
       AreaChild  = ElementArea( NewMesh, NewMesh % Elements(n), &
2537
            NewMesh % Elements(n) % TYPE % NumberOfNodes )
2538

2539
       ChildFactors => NewMesh % Elements(n) % BoundaryInfo % RadiationFactors
2540

2541
       CALL UpdateChildFactors( AreaParent, Factors, &
2542
            AreaChild, ChildFactors, Children )
2543

2544
       n = Children(i,2)
2545

2546
       IF ( n > 0 ) THEN
2547
          AreaChild = ElementArea( NewMesh, NewMesh % Elements(n), &
2548
               NewMesh % Elements(n) % TYPE % NumberOfNodes )
2549

2550
          ChildFactors => NewMesh % Elements(n) % &
2551
               BoundaryInfo % RadiationFactors
2552

2553
          CALL UpdateChildFactors( AreaParent, Factors, &
2554
               AreaChild, ChildFactors, Children )
2555
       END IF
2556
    END DO
2557
!------------------------------------------------------------------------------
2558
  END SUBROUTINE UpdateGebhartFactors
2559
!------------------------------------------------------------------------------
2560

2561

2562
!------------------------------------------------------------------------------
2563
  SUBROUTINE UpdateChildFactors( Area, Factors, AreaNew, NewFactors,Children )
2564
!------------------------------------------------------------------------------
2565
    REAL(KIND=dp) :: Area, AreaNew
2566
    INTEGER :: Children(:,:)
2567
    TYPE(Factors_t), POINTER :: Factors, NewFactors
2568
!------------------------------------------------------------------------------
2569
    INTEGER k,n,TARGET,NEW
2570
!------------------------------------------------------------------------------
2571
    ALLOCATE(NewFactors % Factors(NewFactors % NumberOfFactors))
2572
    ALLOCATE(NewFactors % Elements(NewFactors % NumberOfFactors))
2573

2574
    NEW = 0
2575
    DO k=1,Factors % NumberOfFactors
2576
       TARGET = Factors % Elements(k)
2577
       n = Children(TARGET,1)
2578

2579
       NEW = NEW + 1
2580
       NewFactors % Elements(NEW) = n
2581
       NewFactors % Factors(NEW)  = AreaNew * Factors % Factors(k) / Area
2582
            
2583
       n = Children(TARGET,2)
2584
       IF ( n > 0 ) THEN
2585
          NEW = NEW + 1
2586
          NewFactors % Elements(NEW) = n
2587
          NewFactors % Factors(NEW)  = AreaNew * Factors % Factors(k) / Area
2588
       END IF
2589
    END DO
2590
!------------------------------------------------------------------------------
2591
  END SUBROUTINE UpdateChildFactors
2592
!------------------------------------------------------------------------------
2593
 
2594
!------------------------------------------------------------------------------
2595
 END SUBROUTINE RefineMesh
2596
!------------------------------------------------------------------------------
2597

2598

2599
!------------------------------------------------------------------------------
2600
  FUNCTION ComputeError( Model, ErrorIndicator, RefMesh,  &
2601
       Quant, Perm, InsideResidual, EdgeResidual, BoundaryResidual ) RESULT(MaxError)
2602
!------------------------------------------------------------------------------
2603
    USE CRSMatrix
2604
    IMPLICIT NONE
2605

2606
    TYPE(Mesh_t), POINTER :: RefMesh
2607
    TYPE(Model_t) :: Model
2608
    INTEGER :: Perm(:)
2609
    REAL(KIND=dp) :: ErrorIndicator(:), Quant(:), MaxError
2610

2611
    INTERFACE
2612
       SUBROUTINE BoundaryResidual( Model,Edge,Mesh,Quant,Perm,Gnorm, Indicator)
2613
         USE Types
2614
         TYPE(Element_t), POINTER :: Edge
2615
         TYPE(Model_t) :: Model
2616
         TYPE(Mesh_t), POINTER :: Mesh
2617
         REAL(KIND=dp) :: Quant(:), Indicator(2), Gnorm
2618
         INTEGER :: Perm(:)
2619
       END SUBROUTINE BoundaryResidual
2620

2621
       SUBROUTINE EdgeResidual( Model,Edge,Mesh,Quant,Perm, Indicator)
2622
         USE Types
2623
         TYPE(Element_t), POINTER :: Edge
2624
         TYPE(Model_t) :: Model
2625
         TYPE(Mesh_t), POINTER :: Mesh
2626
         REAL(KIND=dp) :: Quant(:), Indicator(2)
2627
         INTEGER :: Perm(:)
2628
       END SUBROUTINE EdgeResidual
2629

2630
       SUBROUTINE InsideResidual( Model,Element,Mesh,Quant,Perm,Fnorm, Indicator)
2631
         USE Types
2632
         TYPE(Element_t), POINTER :: Element
2633
         TYPE(Model_t) :: Model
2634
         TYPE(Mesh_t), POINTER :: Mesh
2635
         REAL(KIND=dp) :: Quant(:), Indicator(2), Fnorm
2636
         INTEGER :: Perm(:)
2637
       END SUBROUTINE InsideResidual
2638
    END INTERFACE
2639
!------------------------------------------------------------------------------
2640
    TYPE(Element_t), POINTER :: Edge, Face, Boundary, Element
2641
    INTEGER :: i, j, k, Parent
2642
    REAL(KIND=dp), POINTER :: TempIndicator(:,:)
2643
    REAL(KIND=dp) :: LocalIndicator(2), Fnorm, LocalFnorm,s,s1,s2
2644
!------------------------------------------------------------------------------
2645
    CALL FindMeshEdges( RefMesh )
2646

2647
    Fnorm = 0.0_dp
2648
    ErrorIndicator = 0.0_dp
2649

2650
    CALL AllocateArray(TempIndicator,2,SIZE(ErrorIndicator))
2651
    TempIndicator = 0.0_dp
2652
!
2653
!   Bulk equation residuals:
2654
!   ------------------------
2655
    DO i=1,RefMesh % NumberOfBulkElements
2656
       Element => RefMesh % Elements(i)
2657
       CurrentModel % CurrentElement => Element
2658

2659
       CALL InsideResidual( Model, Element, &
2660
             RefMesh, Quant, Perm, LocalFnorm, LocalIndicator )
2661

2662
       Fnorm = Fnorm + LocalFnorm
2663
       TempIndicator(:,i) = TempIndicator(:,i) + LocalIndicator
2664
    END DO
2665

2666

2667
    SELECT CASE( CoordinateSystemDimension())
2668
    CASE(2)
2669
!
2670
!   Edge jumps (2D):
2671
!   ----------------
2672
    DO i = 1,RefMesh % NumberOfEdges
2673
       Edge => RefMesh % Edges(i)
2674
       CurrentModel % CurrentElement => Edge
2675

2676
       IF ( .NOT. ASSOCIATED( Edge % BoundaryInfo ) ) CYCLE
2677

2678
       IF ( ASSOCIATED( Edge % BoundaryInfo % Right ) ) THEN
2679
          CALL EdgeResidual( Model, Edge, RefMesh, Quant, Perm, LocalIndicator )
2680

2681
          Parent = Edge % BoundaryInfo % Left % ElementIndex
2682
          TempIndicator( :,Parent ) = &
2683
               TempIndicator( :,Parent ) + LocalIndicator
2684
          
2685
          Parent = Edge % BoundaryInfo % Right % ElementIndex
2686
          TempIndicator( :,Parent ) = &
2687
               TempIndicator( :,Parent ) + LocalIndicator
2688
       END IF
2689
    END DO
2690

2691
    CASE(3)
2692
!
2693
!   Face jumps (3D):
2694
!   ----------------
2695

2696
    DO i = 1,RefMesh % NumberOfFaces
2697
       Face => RefMesh % Faces(i)
2698
       CurrentModel % CurrentElement => Face
2699

2700
       IF ( .NOT. ASSOCIATED( Face % BoundaryInfo ) ) CYCLE
2701

2702
       IF ( ASSOCIATED( Face % BoundaryInfo % Right ) ) THEN
2703
          CALL EdgeResidual( Model, Face, RefMesh, Quant, Perm, LocalIndicator )
2704

2705
          Parent = Face % BoundaryInfo % Left % ElementIndex
2706
          TempIndicator( :,Parent ) = TempIndicator( :,Parent ) + LocalIndicator
2707
          
2708
          Parent = Face % BoundaryInfo % Right % ElementIndex
2709
          TempIndicator( :,Parent ) = TempIndicator( :,Parent ) + LocalIndicator
2710
       END IF
2711
    END DO
2712
    END SELECT
2713

2714
!
2715
!   Boundary condition residuals:
2716
!   -----------------------------
2717
    DO i = RefMesh % NumberOfBulkElements + 1,  &
2718
           RefMesh % NumberOfBulkElements + RefMesh % NumberOfBoundaryElements
2719

2720
       Boundary => RefMesh % Elements(i)
2721
       CurrentModel % CurrentElement => Boundary
2722

2723
       IF ( Boundary % Type % ElementCode == 101 ) CYCLE
2724

2725
       CALL BoundaryResidual( Model, Boundary, &
2726
             RefMesh, Quant, Perm, LocalFnorm, LocalIndicator )
2727

2728
       Fnorm = Fnorm + LocalFnorm
2729

2730
       IF ( ASSOCIATED( Boundary % BoundaryInfo % Left) ) THEN
2731
         Parent = Boundary % BoundaryInfo % Left % ElementIndex
2732
         IF ( Parent > 0 ) TempIndicator( :,Parent ) = &
2733
              TempIndicator( :,Parent ) + LocalIndicator
2734
       END IF
2735
          
2736
       IF ( ASSOCIATED( Boundary % BoundaryInfo % RIght) ) THEN
2737
         Parent = Boundary % BoundaryInfo % Right % ElementIndex
2738
         IF ( Parent > 0 ) TempIndicator( :,Parent ) = &
2739
              TempIndicator( :,Parent ) + LocalIndicator
2740
       END IF
2741
    END DO
2742

2743

2744
    s1 = SUM(TempIndicator(2,:))
2745
    S2 = SUM(TempIndicator(1,:))
2746
    
2747
    IF(ParEnv % PEs > 1 ) THEN
2748
      s1 = ParallelReduction(s1)
2749
      s2 = ParallelReduction(s2)
2750
      Fnorm = ParallelReduction(Fnorm) 
2751
    END IF
2752
          
2753
    s = SQRT(s1) / SQRT(s2)
2754
    ErrorIndicator = SQRT( TempIndicator(1,:)/(2*s) + s*TempIndicator(2,:)/2 )
2755

2756
    IF ( Fnorm > AEPS ) THEN
2757
       ErrorIndicator = ErrorIndicator / SQRT(Fnorm)
2758
    END IF
2759

2760
    MaxError = MAXVAL( ErrorIndicator )
2761
    IF(ParEnv % PEs>1) MaxError = ParallelReduction(MaxError,2)
2762
        
2763
    DEALLOCATE( TempIndicator )
2764
!------------------------------------------------------------------------------
2765
  END FUNCTION ComputeError
2766
!------------------------------------------------------------------------------
2767

2768
!------------------------------------------------------------------------------
2769
SUBROUTINE FluxRecovery(Model, Solver, Mesh, ErrorIndicator, MaxError)
2770
!------------------------------------------------------------------------------
2771
! Flux recovery & a posteriori error estimation
2772
! The author of the flux recovery part: [email protected]
2773
!------------------------------------------------------------------------------
2774
  IMPLICIT NONE
2775
  TYPE(Model_t) :: Model
2776
  TYPE(Solver_t) :: Solver
2777
  TYPE(Mesh_t), POINTER :: Mesh
2778
  REAL(KIND=dp) :: ErrorIndicator(:), MaxError
2779
!------------------------------------------------------------------------------
2780
  TYPE(Element_t),POINTER :: Element
2781
  TYPE(Variable_t), POINTER :: RTFlux, NodalLoads
2782
  TYPE(ValueList_t), POINTER :: RTSolverPars
2783
  TYPE(Element_t), POINTER :: Parent, Face
2784
  INTEGER, ALLOCATABLE :: VisitsCounter(:)
2785
  INTEGER, ALLOCATABLE, SAVE :: Indices(:), RT_Indices(:)
2786
  INTEGER, POINTER :: FaceMap(:,:)
2787
  LOGICAL, SAVE :: AllocationsDone = .FALSE.
2788
  LOGICAL :: Found, UseReactions, Parallel, PostSmoothing, BDM
2789
  LOGICAL :: ReverseSign(3), OrientationsMatch
2790
  INTEGER :: n, nb, nd, t, active  
2791
  INTEGER :: i, j, k, m, p, ni, nj, nd_rt, istat, ActiveFaceId
2792
  INTEGER :: i_r, j_r
2793
  REAL(KIND=dp) :: UK(3), s, R1, R2, w1, w2, detw, r_i, r_j, hK
2794
  REAL(KIND=dp) :: LinFun(8)   ! The size corresponds to RT_1(K)
2795
  REAL(KIND=dp) :: Err, SolNorm, SolNormEst, Est, APostEst, Est_K
2796
  REAL(KIND=dp), ALLOCATABLE :: RTFluxPost(:), ReactionWeights(:)
2797
  CHARACTER(LEN=:), ALLOCATABLE :: matpar_name, force_name 
2798
!------------------------------------------------------------------------------  
2799

2800
  ! We need a variable 'RTFlux' that is constructed as an approximation in RT_1/BDM
2801
  !
2802
  RTFlux => VariableGet(Mesh % Variables, 'RTFlux')
2803

2804
  IF (ASSOCIATED(RTFlux)) THEN
2805
    IF (.NOT. ASSOCIATED(RTFlux % Solver)) &
2806
        CALL Fatal('FluxRecovery', 'RTFlux variable is not associated with any solver')
2807

2808
    RTSolverPars => GetSolverParams(RTFlux % Solver)
2809
    BDM = GetLogical(RTSolverPars, 'Second Kind Basis', Found)
2810

2811
    ALLOCATE(CHARACTER(MAX_NAME_LEN)::force_name)
2812
    ALLOCATE(CHARACTER(MAX_NAME_LEN)::matpar_name)
2813

2814
    matpar_name = ListGetString(RTSolverPars, 'Material Parameter Name', Found)
2815
    force_name = ListGetString(RTSolverPars, 'Body Force Name', Found)
2816

2817
    n = SIZE(RTFlux % Values)
2818
    ALLOCATE( VisitsCounter(n), STAT=istat )
2819
    VisitsCounter = 0
2820
    RTFlux % Values(:) = 0.0d0
2821

2822
    IF (.NOT. AllocationsDone) THEN
2823
      N = Mesh % MaxElementDOFs
2824
      ALLOCATE(Indices(N), RT_Indices(N), STAT=istat)
2825
      AllocationsDone = .TRUE.
2826
    ELSE
2827
      IF (SIZE(Indices) < Mesh % MaxElementDOFs) THEN
2828
        CALL Fatal('FluxRecovery', 'Mesh changed, the maximun counts of element DOFs are different?')
2829
      END IF
2830
    END IF
2831
  ELSE
2832
    CALL Fatal('FluxRecovery', 'RTFlux variable cannot be found')
2833
  END IF
2834

2835
  
2836
  !------------------------------------------------------------------------------
2837
  ! Step I: 
2838
  ! Compute the values of linear functionals (that is, DOFs) to obtain the elementwise
2839
  ! representation of the flux (stress resultant) in RT_1(K) (or BDM(K)). We add the computed values
2840
  ! to entries of the variable which is suitable for defining a conforming approximation
2841
  ! in RT_1. A later averaging may be applied to ensure that a conforming approximation
2842
  ! in RT_1 is obtained.
2843
  !------------------------------------------------------------------------------
2844
  Active = GetNOFActive()
2845
  DO K=1,Active
2846
    Element => GetActiveElement(K)
2847
    IF (.NOT.(GetElementFamily(Element) == 3)) CYCLE
2848

2849
    n  = Element % Type % NumberOfNodes
2850
    nd = MGetElementDOFs(Indices, Element, Solver)
2851
    nb = GetElementNOFBDOFs(Element, Solver)
2852

2853
    IF (nb > 0) CALL Fatal('FluxRecovery', 'Bubbles in Global System = True assumed')
2854
    
2855
    UK(1:nd) = Solver % Variable % Values(Solver % Variable % Perm(Indices(1:nd)))
2856

2857
    nd_rt = MGetElementDOFs(RT_Indices, Element, USolver = RTFlux % Solver)
2858

2859
    ! Compute the values of DOFs for the representation in RT_1(K)
2860
    !
2861
    CALL EstimateError(UK, Element, n, nd, LinFun = LinFun, MatParName = matpar_name, &
2862
        ForceName = force_name)
2863

2864
    DO i=1,nd_rt
2865
      j = RTFlux % Solver % Variable % Perm(RT_Indices(i))
2866
      RTFlux % Values(j) = RTFlux % Values(j) + LinFun(i)
2867
      VisitsCounter(j) = VisitsCounter(j) + 1
2868
    END DO
2869
  END DO
2870

2871
  ! Averaging and applying the know constraints related to BCs. First, average:
2872
  !
2873
  DO i=1,SIZE(RTFlux % Values)
2874
    IF (VisitsCounter(i) > 1) THEN
2875
      RTFlux % Values(i) = RTFlux % Values(i) / VisitsCounter(i)
2876
    END IF
2877
  END DO
2878

2879
  NodalLoads => NULL()
2880
  NodalLoads => VariableGet(Mesh % Variables, &      
2881
      GetVarName(Solver % Variable) // ' Loads' )
2882

2883
  IF (.NOT. ASSOCIATED(NodalLoads)) THEN
2884
    UseReactions = .FALSE.
2885
  ELSE  
2886
    ! First, check whether there are reactions caused by BCs
2887
    !
2888
    IF (ALLOCATED(Solver % Matrix % ConstrainedDOF)) THEN
2889
      IF (COUNT(Solver % Matrix % ConstrainedDOF) > 0) THEN
2890
        UseReactions = .TRUE.
2891
      ELSE
2892
        UseReactions = .FALSE.
2893
      END IF
2894
    END IF
2895
  END IF
2896

2897
  Apply_reactions: IF (UseReactions) THEN
2898
    !
2899
    ! Create data in order to average reactions 
2900
    !
2901
    ! First we tag DOFs on the model boundary by computing suitable weights
2902
    ! for averaging
2903
    ! TO DO: This should also be done for the DOFs which are associated with
2904
    !        point loads. They need not necessarily be on the boundary.
2905
    !        Smaller arrays could also be allocated by revising the code.
2906
    !
2907
    ALLOCATE(ReactionWeights(SIZE(Solver % Variable % Values)), STAT=istat)
2908
    ReactionWeights = 0
2909
    Parallel = ASSOCIATED(Mesh % ParallelInfo % GInterface)
2910

2911
!    m = 0
2912
    count_shared_nodes: DO K=1, GetNOFBoundaryElements()
2913
      Element => GetBoundaryElement(K)
2914
      IF (.NOT.(GetElementFamily(Element) == 2)) CYCLE
2915
      IF (ActiveBoundaryElement()) THEN
2916
        n = Element % Type % NumberOfNodes
2917
        nd = MGetElementDOFs(Indices)
2918
        
2919
        IF (COUNT(Solver % Matrix % ConstrainedDOF(Solver % Variable % Perm(Indices(1:n)))) == n) THEN
2920
!          m = m + 1
2921
          hK = SQRT((Mesh % Nodes % x(Indices(2)) - Mesh % Nodes % x(Indices(1)))**2 + &
2922
              (Mesh % Nodes % y(Indices(2)) - Mesh % Nodes % y(Indices(1)))**2)
2923
          DO i=1,n
2924
            ReactionWeights(Indices(i)) = ReactionWeights(Indices(i)) + hK
2925
          END DO
2926
        END IF
2927
      END IF
2928
    END DO count_shared_nodes
2929

2930
    ! Now, loop again to create the values of DOFs using the reactions. We already know what elements
2931
    ! should be inspected, so this could be made more efficient ...
2932
    !
2933
    w1 = 0.5d0 * (1.0d0 + 1.0d0/sqrt(3.0d0))
2934
    w2 = 0.5d0 * (1.0d0 - 1.0d0/sqrt(3.0d0))
2935
    detw = w1**2 - w2**2
2936
    
2937
    replace_boundary_dofs: DO K=1, GetNOFBoundaryElements()
2938
      Element => GetBoundaryElement(K)
2939
      IF (.NOT.(GetElementFamily(Element) == 2)) CYCLE
2940
      IF (ActiveBoundaryElement()) THEN
2941
        n = Element % Type % NumberOfNodes
2942
        nd = MGetElementDOFs(Indices)
2943
        
2944
        IF (COUNT(Solver % Matrix % ConstrainedDOF(Solver % Variable % Perm(Indices(1:n)))) == n) THEN
2945
          !
2946
          ! We need the parent to check the sign reversion
2947
          !
2948
          Parent => Element % BoundaryInfo % Left
2949
          IF (.NOT. ASSOCIATED(Parent)) THEN
2950
            Parent => Element % BoundaryInfo % Right
2951
          END IF
2952
          IF (.NOT. ASSOCIATED(Parent)) Call Fatal('FluxRecovery', 'A parent element is not defined')  
2953
          !
2954
          ! Identify the face representing the element among the faces of 
2955
          ! the parent element:
2956
          !
2957
          CALL PickActiveFace(Mesh, Parent, Element, Face, ActiveFaceId)
2958
          IF (ActiveFaceId == 0) Call Fatal('FluxRecovery', 'Cannot determine an element face')
2959

2960
          CALL FaceElementOrientation(Parent, ReverseSign, ActiveFaceId)
2961
          
2962
          IF (IsLeftHanded(Parent)) THEN
2963
            IF (ReverseSign(ActiveFaceId)) THEN
2964
              s = 1.0d0
2965
            ELSE
2966
              s = -1.0d0
2967
            END IF
2968
          ELSE
2969
            IF (ReverseSign(ActiveFaceId)) THEN
2970
              s = -1.0d0
2971
            ELSE
2972
              s = 1.0d0
2973
            END IF
2974
          END IF
2975
          
2976
          FaceMap => GetEdgeMap(GetElementFamily(Parent))
2977
          i_r = Element % NodeIndexes(1)
2978
          j_r = Element % NodeIndexes(2)
2979
          IF (ReverseSign(ActiveFaceId)) THEN
2980
            i = Parent % NodeIndexes(FaceMap(ActiveFaceId,2))
2981
          ELSE
2982
            i = Parent % NodeIndexes(FaceMap(ActiveFaceId,1))
2983
          END IF
2984
          OrientationsMatch = i_r == i
2985

2986
          hK = SQRT((Mesh % Nodes % x(Indices(2)) - Mesh % Nodes % x(Indices(1)))**2 + &
2987
              (Mesh % Nodes % y(Indices(2)) - Mesh % Nodes % y(Indices(1)))**2)
2988
          
2989
          IF (OrientationsMatch) THEN
2990
            R1 = s * NodalLoads % Values(Solver % Variable % Perm(i_r)) * hk / ReactionWeights(i_r)
2991
            R2 = s * NodalLoads % Values(Solver % Variable % Perm(j_r)) * hk / ReactionWeights(j_r)
2992
          ELSE
2993
            R1 = s * NodalLoads % Values(Solver % Variable % Perm(j_r)) * hK / ReactionWeights(j_r)
2994
            R2 = s * NodalLoads % Values(Solver % Variable % Perm(i_r)) * hK / ReactionWeights(i_r)
2995
          END IF
2996

2997
          nd_rt = MGetElementDOFs(RT_Indices, Parent, USolver = RTFlux % Solver)
2998

2999
!          print *, '=== the first DOF now/before', R1, &
3000
!              RTFlux % Values(RTFlux % Solver % Variable % Perm(RT_Indices(2*ActiveFaceId - 1)))
3001
!          print *, '=== the second DOF now/before', R2, &
3002
!              RTFlux % Values(RTFlux % Solver % Variable % Perm(RT_Indices(2*ActiveFaceId)))
3003

3004
          ! Finally use the reactions to replace the values of DOFs on the boundary
3005
          !
3006
          IF (BDM) THEN
3007
            RTFlux % Values(RTFlux % Solver % Variable % Perm(RT_Indices(2*ActiveFaceId - 1))) = &
3008
                1.0d0/detw * (w1*R1 - w2*R2)
3009
            RTFlux % Values(RTFlux % Solver % Variable % Perm(RT_Indices(2*ActiveFaceId))) = &
3010
                1.0d0/detw * (-w2*R1 + w1*R2)
3011
          ELSE
3012
            RTFlux % Values(RTFlux % Solver % Variable % Perm(RT_Indices(2*ActiveFaceId - 1))) = R1
3013
            RTFlux % Values(RTFlux % Solver % Variable % Perm(RT_Indices(2*ActiveFaceId))) = R2
3014
          END IF
3015
        END IF
3016
      END IF
3017
    END DO replace_boundary_dofs
3018
  END IF Apply_reactions
3019

3020
  !------------------------------------------------------------------------------
3021
  ! Step II: 
3022
  ! Equilibrate fluxes (stress resultants). This brings us back to the recovery
3023
  ! which is defined only in the local RT space.
3024
  ! Here the exact solution is also used to study the accuracy.
3025
  !------------------------------------------------------------------------------
3026

3027
  PostSmoothing = .FALSE. ! Eliminated as this doesn't seem to be beneficial
3028
  IF (PostSmoothing) THEN
3029
    ALLOCATE(RTFluxPost(SIZE(RTFlux % Values)))
3030
    RTFluxPost(:) = 0.0d0
3031
    VisitsCounter(:) = 0
3032
  END IF
3033
  
3034
  Err = 0.0d0
3035
  SolNorm = 0.0d0
3036
  Est = 0.0d0
3037
  APostEst = 0.0d0
3038
  SolNormEst = 0.0d0
3039
  
3040
  Elementwise_equilibration: DO K=1,Active
3041
    Element => GetActiveElement(K)
3042
    IF (.NOT. (GetElementFamily(Element) == 3)) CYCLE
3043

3044
    n  = Element % Type % NumberOfNodes
3045
    nd = MGetElementDOFs(Indices)
3046
    nb = GetElementNOFBDOFs()
3047
    IF (nb > 0) CALL Fatal('FluxRecovery', 'Bubbles in Global System = True assumed')
3048
    
3049
    UK(1:nd) = Solver % Variable % Values( Solver % Variable % Perm(Indices(1:nd)) )
3050

3051
    nd_rt = MGetElementDOFs(RT_Indices, Element, USolver = RTFlux % Solver)
3052

3053
    ! At the calling time, we still have the RT DOFs respecting the global continuity requirements
3054
    !
3055
    DO i=1,nd_rt
3056
      j = RTFlux % Solver % Variable % Perm(RT_Indices(i))
3057
      LinFun(i) = RTFlux % Values(j)
3058
    END DO
3059
    
3060
    CALL EstimateError(UK, Element, n, nd, PostLinFun = LinFun, &
3061
        APostEst_K = Est_K, SolNormEst=SolNormEst, MatParName = matpar_name, &
3062
        ForceName = force_name)
3063

3064
    ! The following call could be used to run a verification case:
3065
!    CALL EstimateError(UK, Element, n, nd, Err, SolNorm, Est, APostEst, PostLinFun = LinFun, &
3066
!        APostEst_K = Est_K, SolNormEst=SolNormEst, MatParName = matpar_name, &
3067
!        ForceName = force_name)
3068

3069
    ErrorIndicator(K) = SQRT(Est_K)
3070

3071
    IF (PostSmoothing) THEN
3072
      DO i=1,nd_rt
3073
        j = RTFlux % Solver % Variable % Perm(RT_Indices(i))
3074
        RTFluxPost(j) = RTFluxPost(j) + LinFun(i)
3075
        VisitsCounter(j) = VisitsCounter(j) + 1
3076
      END DO
3077
    END IF
3078
    
3079
  END DO Elementwise_equilibration
3080

3081
  IF (SolNormEst > AEPS) ErrorIndicator = (1.0d0 / SQRT(SolNormEst)) * ErrorIndicator
3082
  MaxError = MAXVAL(ErrorIndicator)
3083
  IF (ParEnv % PEs>1) MaxError = ParallelReduction(MaxError,2)
3084
  
3085
!  WRITE (*, '(A,E16.8)') 'Solution Norm = ', SQRT(ParallelReduction(SolNorm))
3086
!  WRITE (*, '(A,E16.8)') 'Error Norm = ', SQRT(ParallelReduction(Err))/SQRT(ParallelReduction(SolNorm))
3087
!  WRITE (*, '(A,E16.8)') 'Recovery Error Norm = ', SQRT(ParallelReduction(Est))/SQRT(ParallelReduction(SolNorm))
3088
!  WRITE (*, '(A,E16.8)') 'A posteriori Error = ', SQRT(ParallelReduction(APostEst))/SQRT(ParallelReduction(SolNorm))
3089
!  WRITE (*, '(A,E16.8)') 'A posteriori Error Est = ', SQRT(SUM(ErrorIndicator**2)) 
3090
!  WRITE (*, '(A,E16.8)') 'Efficiency Factor = ', SQRT(ParallelReduction(APostEst))/SQRT(ParallelReduction(Err))
3091

3092
  
3093
  !
3094
  ! The following computation would average again to obtain the recovery field in the global RT_1 space
3095
  !
3096
  post_smoothing: IF (PostSmoothing) THEN
3097

3098
    Err = 0.0d0
3099
    Est = 0.0d0
3100
    APostEst = 0.0d0
3101
    SolNorm = 0.0d0
3102

3103
    ! Post averaging
3104
    !
3105
    DO i=1,SIZE(RTFlux % Values)
3106
      IF (VisitsCounter(i) > 1) THEN
3107
        RTFluxPost(i) = RTFluxPost(i) / VisitsCounter(i)
3108
      END IF
3109
    END DO
3110

3111
    Err = 0.0d0
3112
    Est = 0.0d0
3113
    APostEst = 0.0d0
3114
    SolNorm = 0.0d0
3115

3116
    DO K=1,Active
3117
      Element => GetActiveElement(K)
3118
      IF ( .NOT. (GetElementFamily(Element) == 3) ) CYCLE
3119

3120
      n  = Element % Type % NumberOfNodes
3121
      nd = MGetElementDOFs(Indices)
3122
      nb = GetElementNOFBDOFs()
3123
      IF (nb > 0) CALL Fatal('FluxRecovery', 'Bubbles in Global System = True assumed')
3124

3125
      UK(1:nd) = Solver % Variable % Values( Solver % Variable % Perm(Indices(1:nd)) )
3126

3127
      nd_rt = MGetElementDOFs(RT_Indices, Element, USolver = RTFlux % Solver)
3128

3129
      DO i=1,nd_rt
3130
        j = RTFlux % Solver % Variable % Perm(RT_Indices(i))
3131
        LinFun(i) = RTFluxPost(j)
3132
      END DO
3133

3134
      CALL EstimateError(UK, Element, n, nd, Err, SolNorm, Est, APostEst, PostLinFun = LinFun, &
3135
          UseGiven = .TRUE., MatParName = matpar_name, ForceName = force_name)
3136

3137
    END DO
3138

3139
    WRITE (*, '(A,E16.8)') 'Solution Norm = ', SQRT(ParallelReduction(SolNorm))
3140
    WRITE (*, '(A,E16.8)') 'Error Norm = ', SQRT(ParallelReduction(Err))/SQRT(ParallelReduction(SolNorm))
3141
    WRITE (*, '(A,E16.8)') 'Estimated Error Norm = ', SQRT(ParallelReduction(Est))/SQRT(ParallelReduction(SolNorm))
3142
    WRITE (*, '(A,E16.8)') 'A posteriori Error Est = ', SQRT(ParallelReduction(APostEst))/SQRT(ParallelReduction(SolNorm))
3143
    WRITE (*, '(A,E16.8)') 'Efficiency Factor = ', SQRT(ParallelReduction(APostEst))/SQRT(ParallelReduction(Err))
3144

3145
  END IF post_smoothing
3146

3147
  
3148
CONTAINS
3149

3150
!------------------------------------------------------------------------------
3151
  SUBROUTINE EstimateError(UK, Element, n, nd, Err, SolNorm, Est, APostEst, &
3152
      LinFun, PostLinFun, UseGiven, APostEst_K, SolNormEst, &
3153
      MatParName, ForceName)
3154
!------------------------------------------------------------------------------
3155
    REAL(KIND=dp), INTENT(IN) :: UK(:)
3156
    TYPE(Element_t), POINTER, INTENT(IN) :: Element    
3157
    INTEGER, INTENT(IN) :: n, nd
3158
    ! The following four arguments are related only to running a code verification case:
3159
    REAL(KIND=dp), OPTIONAL, INTENT(INOUT) :: Err, SolNorm, Est, APostEst
3160
    REAL(KIND=dp), OPTIONAL, INTENT(INOUT) :: LinFun(8)
3161
    REAL(KIND=dp), OPTIONAL, INTENT(INOUT) :: PostLinFun(8)
3162
    LOGICAL, OPTIONAL, INTENT(IN) :: UseGiven
3163
    REAL(KIND=dp), OPTIONAL, INTENT(INOUT) :: APostEst_K, SolNormEst
3164
    CHARACTER(LEN=*) :: MatParName, ForceName
3165
!--------------------------------------------------------------------------------
3166
    TYPE(ValueList_t), POINTER :: BodyForce, Material
3167
    
3168
    REAL(KIND=dp) :: FBasis(8,3), DivFBasis(8), Mass(11,11), RHS(11), c(11), d(11), s
3169
    REAL(KIND=dp) :: Basis(nd), dBasis(nd,3), DetJ, xt, uq, vq, weight, EA, f
3170
    REAL(KIND=dp) :: U, gradU(3), Uh, gradUh(3), Nh(3), intf, divN, totflux, dc
3171
    REAL(KIND=dp) :: testfun(2), diff_coeff(n), Load(n)
3172
    REAL(KIND=dp) :: faceflux(3), faceweights(3), facedelta(3)
3173
    REAL(KIND=dp) :: w1, w2
3174
    LOGICAL :: Stat, Found
3175
    LOGICAL :: ReverseSign(3), LeftHanded, Parallel, UseLM
3176
    LOGICAL :: FirstOrderEquilibration 
3177
    INTEGER, POINTER :: EdgeMap(:,:)
3178
    INTEGER :: t, i, j, p, q, ni, nj, np, FDOFs, DOFs, ElementOrder
3179
    TYPE(GaussIntegrationPoints_t) :: IP
3180

3181
    TYPE(Nodes_t), SAVE :: Nodes
3182
!------------------------------------------------------------------------------
3183
    IF (BDM) THEN
3184
      FDOFs = 6
3185
      ElementOrder = 1
3186
      FirstOrderEquilibration = .FALSE.
3187
    ELSE
3188
      FDOFs = 8
3189
      ElementOrder = 2
3190
      FirstOrderEquilibration = .TRUE.
3191
    END IF
3192

3193
    UseLM = .TRUE.
3194
    IF (UseLM) THEN
3195
      IF (BDM) THEN
3196
        DOFs = FDOFs + 1
3197
      ELSE
3198
        DOFs = FDOFs + 3
3199
      END IF
3200
    ELSE
3201
      DOFs = FDOFs
3202
    END IF
3203
       
3204
    CALL GetElementNodes( Nodes )
3205
    IP = GaussPoints(Element, EdgeBasis=.TRUE., PReferenceElement=.TRUE., EdgeBasisDegree=2)
3206

3207
    Material => GetMaterial()
3208
    diff_coeff(1:n) = GetReal(Material, MatParName)
3209

3210
    IF (PRESENT(APostEst_K)) APostEst_K = 0.0d0
3211
    
3212
    c = 0.0d0
3213
    IF (PRESENT(UseGiven) .AND. PRESENT(PostLinFun)) THEN
3214
      IF (UseGiven) THEN
3215
        c(1:FDOFs) = PostLinFun(1:FDOFs)
3216
        GOTO 303
3217
      END IF
3218
    END IF
3219

3220
    Load = 0.0d0
3221
    BodyForce => GetBodyForce()
3222
    IF (ASSOCIATED(BodyForce)) &
3223
        Load(1:n) = GetReal(BodyForce, ForceName, Found)
3224
    
3225
    IF (PRESENT(LinFun)) THEN
3226
      !
3227
      ! Compute a local representation of the flux in the elementwise RT_1(K) space.
3228
      ! The elementwise weak formulation is of the form
3229
      !           1/k (N,v) = -(u,div v) + <u,v.n> 
3230
      !
3231
      Parallel = ASSOCIATED(Mesh % ParallelInfo % GInterface)
3232
      EdgeMap => GetEdgeMap(3)
3233
      CALL FaceElementOrientation(Element, ReverseSign)
3234

3235
      Mass = 0.0_dp
3236
      RHS = 0.0_dp
3237
      IF (BDM) THEN 
3238
        w1 = 0.5d0 * (1.0d0 + 1.0d0/sqrt(3.0d0))
3239
        w2 = 0.5d0 * (1.0d0 - 1.0d0/sqrt(3.0d0))
3240
      !ELSE
3241
      !  w1 = 1.0d0
3242
      !  w2 = 0.0d0
3243
      END IF
3244
      
3245
      DO t=1,IP % n
3246

3247
        stat = FaceElementInfo(Element, Nodes, IP % U(t), IP % V(t), &
3248
            IP % W(t), detF=detJ, Basis=Basis, FBasis=FBasis, &
3249
            DivFBasis=DivFBasis, BDM=BDM, BasisDegree=ElementOrder, ApplyPiolaTransform=.TRUE., &
3250
            LeftHanded=LeftHanded)
3251

3252
        Weight = IP % s(t) * DetJ
3253

3254
        Uh = SUM(UK(1:nd) * Basis(1:nd))
3255
        EA = SUM(diff_coeff(1:n) * Basis(1:n))
3256
        
3257
        DO p=1,FDOFs
3258
          DO q=1,FDOFs
3259
            Mass(p,q) = Mass(p,q) + SUM(FBasis(q,1:2) * FBasis(p,1:2)) * Weight / EA
3260
          END DO
3261
          RHS(p) = RHS(p) - Weight * Uh * DivFBasis(p)
3262
        END DO
3263
        
3264
        IF (UseLM) THEN
3265
          !
3266
          ! Enforce the zeroth-order equilibration by using a Lagrange multiplier
3267
          !
3268
          f = SUM(Load(1:n) * Basis(1:n))
3269
          RHS(FDOFs+1) = RHS(FDOFs+1) + Weight * f
3270

3271
          IF (FirstOrderEquilibration) THEN
3272
            !
3273
            ! Enforce the first-order equilibration by using a Lagrange multiplier
3274
            !
3275
            testfun(1) = Basis(2) - Basis(1)
3276
            testfun(2) = Basis(3) - Basis(1)
3277
            DO p=1,FDOFs
3278
              Mass(10,p)= Mass(10,p) - divFBasis(p) * testfun(1) * Weight
3279
              Mass(11,p)= Mass(11,p) - divFBasis(p) * testfun(2) * Weight
3280
              RHS(10) = RHS(10) + f * testfun(1) * Weight
3281
              RHS(11) = RHS(11) + f * testfun(2) * Weight
3282
              Mass(p,10)= Mass(p,10) - divFBasis(p) * testfun(1) * Weight
3283
              Mass(p,11)= Mass(p,11) - divFBasis(p) * testfun(2) * Weight
3284
            END DO
3285
          END IF
3286
        END IF
3287
      END DO
3288

3289
      ! Assemble the boundary terms:
3290
      ! Loop over all faces (here edges)
3291
      !
3292
      DO p=1,3
3293
        ! Check whether the sign is reversed
3294
        IF (LeftHanded) THEN
3295
          IF (ReverseSign(p)) THEN
3296
            s = 1.0d0
3297
          ELSE
3298
            s = -1.0d0
3299
          END IF
3300
        ELSE
3301
          IF (ReverseSign(p)) THEN
3302
            s = -1.0d0
3303
          ELSE
3304
            s = 1.0d0
3305
          END IF
3306
        END IF
3307
        ! Check the order of basis functions
3308
        i = EdgeMap(p,1)
3309
        j = EdgeMap(p,2)
3310
        ni = Element % NodeIndexes(i)
3311
        IF (Parallel) ni = Mesh % ParallelInfo % GlobalDOFs(ni)             
3312
        nj = Element % NodeIndexes(j)
3313
        IF (Parallel) nj = Mesh % ParallelInfo % GlobalDOFs(nj)
3314

3315
        IF (BDM) THEN
3316
          IF (nj<ni) THEN
3317
            RHS(2*p-1) = RHS(2*p-1) + s * (w2 * UK(i) + w1 * UK(j))
3318
            RHS(2*p) = RHS(2*p) + s * (w1 * UK(i) + w2 * UK(j))
3319
          ELSE
3320
            RHS(2*p-1) = RHS(2*p-1) + s * (w1 * UK(i) + w2 * UK(j))
3321
            RHS(2*p) = RHS(2*p) + s * (w2 * UK(i) + w1 * UK(j))
3322
          END IF
3323
        ELSE
3324
          ! The value of RHS vector is just the nodal DOF up to the sign
3325
          IF (nj<ni) THEN
3326
            ! The first weight function is s * basis(j)
3327
            RHS(2*p-1) = RHS(2*p-1) + s * UK(j)
3328
            RHS(2*p) = RHS(2*p) + s * UK(i)
3329
          ELSE
3330
            RHS(2*p-1) = RHS(2*p-1) + s * UK(i)
3331
            RHS(2*p) = RHS(2*p) + s * UK(j)
3332
          END IF
3333
        END IF
3334
          
3335
        IF (UseLM) THEN
3336
          Mass(FDOFs+1,2*p-1) = -s
3337
          Mass(FDOFs+1,2*p) = -s
3338
          Mass(2*p-1,FDOFs+1) = -s
3339
          Mass(2*p,FDOFs+1) = -s
3340
        END IF
3341
        
3342
      END DO
3343

3344
      ! Finally solve the local representation of the flux
3345
      CALL InvertMatrix(Mass(1:DOFs,1:DOFs), DOFs)
3346
      c(1:DOFs) = MATMUL(MASS(1:DOFs,1:DOFs), RHS(1:DOFs))
3347
      LinFun(1:FDOFs) = c(1:FDOFs)
3348
      RETURN
3349
    END IF
3350

3351
    IF (PRESENT(PostLinFun)) THEN
3352
      c(1:FDOFs) = PostLinFun(1:FDOFs)
3353

3354
      EdgeMap => GetEdgeMap(3)
3355
      CALL FaceElementOrientation(Element, ReverseSign)
3356
      Parallel = ASSOCIATED(Mesh % ParallelInfo % GInterface)
3357
      
3358
      IF (.NOT. BDM) THEN
3359
        !
3360
        ! Perform the zeroth-order equilibration (note that BDM does not seem to benefit
3361
        ! from this):
3362
        !
3363
        intf = 0.0d0
3364
        divN = 0.0d0
3365
        
3366
        DO t=1,IP % n
3367

3368
          stat = FaceElementInfo(Element, Nodes, IP % U(t), IP % V(t), &
3369
              IP % W(t), detF=detJ, Basis=Basis, FBasis=FBasis, &
3370
              DivFBasis=DivFBasis, BDM=BDM, BasisDegree=ElementOrder, ApplyPiolaTransform=.TRUE., &
3371
              LeftHanded=LeftHanded)
3372

3373
          Weight = IP % s(t) * DetJ
3374
          f = SUM(Load(1:n) * Basis(1:n))
3375
          intf = intf + f * Weight
3376
        END DO
3377

3378
        ! Loop over all faces (here edges)
3379
        !
3380
        totflux = 0.0d0
3381
        faceflux = 0.0d0
3382
        DO p=1,3
3383
          IF (LeftHanded) THEN
3384
            IF (ReverseSign(p)) THEN
3385
              s = 1.0d0
3386
            ELSE
3387
              s = -1.0d0
3388
            END IF
3389
          ELSE
3390
            IF (ReverseSign(p)) THEN
3391
              s = -1.0d0
3392
            ELSE
3393
              s = 1.0d0
3394
            END IF
3395
          END IF
3396
          totflux = totflux + s*(c(2*p-1)+c(2*p))
3397
          faceflux(p) = faceflux(p) + s*(c(2*p-1)+c(2*p))
3398
!          IF (c(2*p-1)*c(2*p) < 0.0 .AND. -c(2*p-1)*c(2*p)/(MAXVAL(ABS(c(:))))**2 > 1.0d-16) THEN
3399
!            CALL Warn('FluxRecovery', 'No consensus on flux sign')
3400
!            print *, 'c1', c(2*p-1)
3401
!            print *, 'c2', c(2*p)
3402
!          END IF
3403
        END DO
3404

3405
        DO p=1,3
3406
          FaceWeights(p) = abs(faceflux(p))/sum(abs(faceflux(:)))
3407
        END DO
3408
        
3409
!        print *, 'total flux out', totflux
3410
!        print *, 'sources tot', -intf
3411
!        print *, 'change flux out by an amount', -intf - totflux
3412
!        dc = (-intf - totflux)/6.0d0
3413
        
3414
        DO p=1,3
3415
          facedelta(p) = FaceWeights(p) * (-intf - totflux)   
3416
        END DO        
3417
        !
3418
        ! The redefinition to ensure the zeroth-order equilibration:
3419
        !
3420
        DO p=1,3
3421
          IF (LeftHanded) THEN
3422
            IF (ReverseSign(p)) THEN
3423
              s = 1.0d0
3424
            ELSE
3425
              s = -1.0d0
3426
            END IF
3427
          ELSE
3428
            IF (ReverseSign(p)) THEN
3429
              s = -1.0d0
3430
            ELSE
3431
              s = 1.0d0
3432
            END IF
3433
          END IF
3434
          
3435
!          d(2*p-1) = c(2*p-1) + s*dc
3436
!          d(2*p) = c(2*p) + s*dc
3437
          d(2*p-1) = c(2*p-1) + s*facedelta(p)*0.5d0
3438
          d(2*p) = c(2*p) + s*facedelta(p)*0.5d0
3439
        END DO
3440

3441
        post_check: IF (.TRUE.) THEN
3442
          totflux = 0.0d0
3443
          DO p=1,3
3444
            IF (LeftHanded) THEN
3445
              IF (ReverseSign(p)) THEN
3446
                s = 1.0d0
3447
              ELSE
3448
                s = -1.0d0
3449
              END IF
3450
            ELSE
3451
              IF (ReverseSign(p)) THEN
3452
                s = -1.0d0
3453
              ELSE
3454
                s = 1.0d0
3455
              END IF
3456
            END IF
3457
            
3458
            totflux = totflux + s*(d(2*p-1)+d(2*p))
3459
          END DO
3460

3461
          IF (ABS(intf + totflux) > 1.0d1 * AEPS) THEN 
3462
            print *, 'Warning: change flux out by an amount', -intf - totflux
3463
          END IF
3464
        END IF post_check
3465

3466
        c(1:6) = d(1:6)
3467
      END IF
3468

3469
      ! Apply the first-order equilibration
3470
      !
3471
      IF (FirstOrderEquilibration) THEN
3472
        Mass = 0.0_dp
3473
        RHS = 0.0_dp
3474
        
3475
        DO t=1,IP % n
3476
          
3477
          stat = FaceElementInfo(Element, Nodes, IP % U(t), IP % V(t), &
3478
              IP % W(t), detF=detJ, Basis=Basis, FBasis=FBasis, &
3479
              DivFBasis=DivFBasis, BasisDegree=2, ApplyPiolaTransform=.TRUE.)
3480

3481
          Weight = IP % s(t) * DetJ
3482
          f = SUM(Load(1:n) * Basis(1:n))
3483
          
3484
          testfun(1) = Basis(2) - Basis(1)
3485
          testfun(2) = Basis(3) - Basis(1)
3486

3487
          DO p=1,2
3488
            DO q=1,2
3489
              Mass(p,q) = Mass(p,q) - divFBasis(6+q) * testfun(p) * Weight
3490
            END DO
3491
            RHS(p) = RHS(p) + f * testfun(p) * Weight
3492
            DO q=1,6
3493
              RHS(p) = RHS(p) + c(q) * divFBasis(q) * testfun(p) * Weight
3494
            END DO
3495
          END DO
3496
        END DO
3497
        
3498
        CALL InvertMatrix(Mass(1:2,1:2),2)
3499
        c(7:8) = MATMUL(MASS(1:2,1:2), RHS(1:2))
3500
      END IF
3501
      PostLinFun(1:FDOFs) = c(1:FDOFs)
3502
    END IF
3503
    
3504
    ! Compute a posteriori estimate:
3505
303 CONTINUE
3506
    DO t=1,IP % n
3507

3508
      ! Now get the face basis functions so that we can evaluate
3509
      ! the flux in terms of them
3510
      !
3511
      stat = FaceElementInfo(Element, Nodes, IP % U(t), IP % V(t), &
3512
            IP % W(t), detF=detJ, Basis=Basis, FBasis=FBasis, &
3513
            DivFBasis=DivFBasis, dBasisdx=dBasis, BDM=BDM, BasisDegree=ElementOrder, &
3514
            ApplyPiolaTransform=.TRUE.)
3515

3516
      Weight = IP % s(t) * DetJ
3517
      
3518
      !  The exact solution:
3519
      IF (PRESENT(Err) .OR. PRESENT(SolNorm)) THEN
3520
        xt = SUM(Basis(1:n) * Nodes % x(1:n))
3521
        U = 1.0d0 - xt**2/9.0d0
3522
        gradU(:) = 0.0
3523
        gradU(1) = -2.0d0 * xt/9.0d0
3524
      END IF
3525

3526
      Uh = SUM(UK(1:nd) * Basis(1:nd))
3527
      DO i=1,2
3528
        gradUh(i) = SUM(UK(1:nd) * dBasis(1:nd,i))
3529
      END DO
3530
      gradUh(3) = 0.0d0
3531

3532
      EA = SUM(diff_coeff(1:n) * Basis(1:n))
3533
      Nh = 0.0d0
3534
      DO i=1,2
3535
        Nh(i) = SUM( c(1:FDOFs) * FBasis(1:FDOFs,i) )
3536
      END DO
3537
      
3538
      IF (.TRUE.) THEN
3539
        ! The error of flux
3540
        IF (PRESENT(Err)) Err = Err + &
3541
            SUM((EA * gradU(1:2) - EA * gradUh(1:2))**2) * Weight
3542
        IF (PRESENT(SolNorm)) SolNorm = SolNorm + &
3543
            SUM((EA * gradU(1:2))**2) * Weight
3544
        IF (PRESENT(SolNormEst)) SolNormEst = SolNormEst + &
3545
            SUM((Nh(1:2))**2) * Weight
3546
        
3547
        ! A posteriori estimate based on the recovery:
3548
        !
3549
        IF (PRESENT(Est)) Est = Est + SUM((EA * gradU(1:2) - Nh(1:2))**2) * Weight
3550
        IF (PRESENT(APostEst)) APostEst = APostEst + SUM((EA * gradUh(1:2) - Nh(1:2))**2) * Weight
3551
        IF (PRESENT(APostEst_K)) APostEst_K = APostEst_K + SUM((EA * gradUh(1:2) - Nh(1:2))**2) * Weight
3552
      ELSE
3553
        ! The error of the field
3554
        IF (PRESENT(Err)) Err = Err + &
3555
            (U - Uh)**2 * Weight
3556
        IF (PRESENT(SolNorm)) SolNorm = SolNorm + &
3557
            U**2 * Weight
3558
      END IF
3559
    END DO
3560
!------------------------------------------------------------------------------
3561
  END SUBROUTINE EstimateError
3562
!------------------------------------------------------------------------------
3563
  
3564
!------------------------------------------------------------------------------
3565
  FUNCTION IsLeftHanded(Element) RESULT(LeftHanded)
3566
!------------------------------------------------------------------------------
3567
    TYPE(Element_t), POINTER, INTENT(IN) :: Element
3568
    LOGICAL :: LeftHanded
3569
!------------------------------------------------------------------------------
3570
    TYPE(Nodes_t), SAVE :: Nodes
3571
    LOGICAL :: Stat
3572
    REAL(KIND=dp) :: Basis(Element % Type % NumberOfNOdes), DetJ
3573
    REAL(KIND=dp) :: FBasis(8,3)
3574
!------------------------------------------------------------------------------
3575
    CALL GetElementNodes(Nodes, Element)    
3576
    
3577
    stat = FaceElementInfo(Element, Nodes, 0.0d0, SQRT(3.0d0)/3.0d0, &
3578
            0.0d0, detF=detJ, Basis=Basis, FBasis=FBasis, &
3579
            BasisDegree=1, ApplyPiolaTransform=.TRUE., &
3580
            LeftHanded=LeftHanded)
3581
!------------------------------------------------------------------------------    
3582
  END FUNCTION IsLeftHanded
3583
!------------------------------------------------------------------------------  
3584

3585
!------------------------------------------------------------------------------
3586
END SUBROUTINE FluxRecovery
3587
!------------------------------------------------------------------------------
3588

3589
!------------------------------------------------------------------------------
3590
END MODULE Adaptive
3591
!-----------------------------------------------------------------------------
3592

3593
SUBROUTINE RefineMeshExt(Model,Solver,Quant,Perm,InsideResidual,EdgeResidual,BoundaryResidual)
3594
  USE adaptive
3595

3596
  IMPLICIT NONE
3597

3598
  TYPE( Model_t ) :: Model
3599
  TYPE(Solver_t), TARGET :: Solver
3600
  REAL(KIND=dp) :: Quant(:)
3601
  INTEGER :: Perm(:)
3602

3603
    INTERFACE
3604
       SUBROUTINE BoundaryResidual( Model,Edge,Mesh,Quant,Perm,Gnorm,Indicator )
3605
          USE Types
3606
          TYPE(Element_t), POINTER :: Edge
3607
          TYPE(Model_t) :: Model
3608
          TYPE(Mesh_t), POINTER :: Mesh
3609
          REAL(KIND=dp) :: Quant(:), Indicator(2), Gnorm
3610
          INTEGER :: Perm(:)
3611
       END SUBROUTINE BoundaryResidual
3612

3613

3614
       SUBROUTINE EdgeResidual( Model,Edge,Mesh,Quant,Perm, Indicator)
3615
          USE Types
3616
          TYPE(Element_t), POINTER :: Edge
3617
          TYPE(Model_t) :: Model
3618
          TYPE(Mesh_t), POINTER :: Mesh
3619
          REAL(KIND=dp) :: Quant(:), Indicator(2)
3620
          INTEGER :: Perm(:)
3621
       END SUBROUTINE EdgeResidual
3622

3623

3624
       SUBROUTINE InsideResidual( Model,Element,Mesh,Quant,Perm,Fnorm, Indicator)
3625
          USE Types
3626
          TYPE(Element_t), POINTER :: Element
3627
          TYPE(Model_t) :: Model
3628
          TYPE(Mesh_t), POINTER :: Mesh
3629
          REAL(KIND=dp) :: Quant(:), Indicator(2), Fnorm
3630
          INTEGER :: Perm(:)
3631
       END SUBROUTINE InsideResidual
3632
    END INTERFACE
3633

3634
  CALL RefineMesh( Model,Solver,Quant,Perm, &
3635
            InsideResidual, EdgeResidual, BoundaryResidual )
3636
END SUBROUTINE RefineMeshExt
3637

3638
!------------------------------------------------------------------------------
3639

3640
!> \} 
3641

3642
Product

Resources

Company
