CoCalc -- AIFlowSolve

GitHub Repository: ElmerCSC/elmerfem
Path: blob/devel/elmerice/Solvers/AIFlowSolve_nlD2.F90
³²⁰⁴ views
1
!/*****************************************************************************/
2
! *
3
! *  Elmer/Ice, a glaciological add-on to Elmer
4
! *  http://elmerice.elmerfem.org
5
! *
6
! * 
7
! *  This program is free software; you can redistribute it and/or
8
! *  modify it under the terms of the GNU General Public License
9
! *  as published by the Free Software Foundation; either version 2
10
! *  of the License, or (at your option) any later version.
11
! * 
12
! *  This program 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
15
! *  GNU General Public License for more details.
16
! *
17
! *  You should have received a copy of the GNU General Public License
18
! *  along with this program (in file fem/GPL-2); if not, write to the 
19
! *  Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, 
20
! *  Boston, MA 02110-1301, USA.
21
! *
22
! *****************************************************************************/
23
! ******************************************************************************
24
! *
25
! *  Authors:  Juha Ruokolainen, Fabien Gillet-Chaulet, Olivier Gagliardini
26
! *  Email:   Juha Ruokolainen
27
! *  Web:     http://elmerice.elmerfem.org
28
! *
29
! *  Original Date: 
30
! *       Date of modification: April 08 => non linear                      
31
! *                             May 09 => N-T (see mail Juha 20 Feb 2006) OG
32
! *                             Dec 15 =>2.5D FlowWidth O. Passalacqua
33
! * 
34
! *****************************************************************************
35
!> Module containing a solver for (primarily thermal) anisotropic flow
36
   RECURSIVE SUBROUTINE AIFlowSolver_nlD2( Model,Solver,dt,TransientSimulation )
37
!------------------------------------------------------------------------------
38

39
    USE DefUtils
40

41
    IMPLICIT NONE
42

43
!------------------------------------------------------------------------------
44
!******************************************************************************
45
!
46
!  Solve stress equations for one timestep
47
!
48
!  ARGUMENTS:
49
!
50
!  TYPE(Model_t) :: Model,  
51
!     INPUT: All model information (mesh,materials,BCs,etc...)
52
!
53
!  TYPE(Solver_t) :: Solver
54
!     INPUT: Linear equation solver options
55
!
56
!  REAL(KIND=dp) :: dt,
57
!     INPUT: Timestep size for time dependent simulations (NOTE: Not used
58
!            currently)
59
!
60
!******************************************************************************
61

62
     TYPE(Model_t)  :: Model
63
     TYPE(Solver_t), TARGET :: Solver
64

65
     LOGICAL ::  TransientSimulation
66
     REAL(KIND=dp) :: dt
67
!------------------------------------------------------------------------------
68
!    Local variables
69
!------------------------------------------------------------------------------
70
     TYPE(Solver_t), POINTER :: PSolver
71

72
     TYPE(Matrix_t),POINTER :: StiffMatrix
73

74
     INTEGER :: i, j, k, l, n, t, iter, NDeg, STDOFs, LocalNodes, istat
75
     INTEGER :: dim, comp 
76

77
     TYPE(ValueList_t),POINTER :: Material, BC, BodyForce
78
     TYPE(Nodes_t) :: ElementNodes
79
     TYPE(Element_t),POINTER :: CurrentElement
80

81
     REAL(KIND=dp) :: RelativeChange, UNorm, PrevUNorm, Gravity(3), &
82
         Normal(3), NewtonTol, NonlinearTol, s, Wn(7), MinSRInvariant
83
         
84

85
     REAL(KIND=dp)  :: NodalStresses(3,3), &
86
       NodalStrainRate(3,3),  NodalSpin(3,3)   
87

88
     REAL(KIND=dp), ALLOCATABLE :: Basis(:),ddBasisddx(:,:,:)
89
     REAL(KIND=dp), ALLOCATABLE :: dBasisdx(:,:), SlipCoeff(:,:)
90
     REAL(KIND=dp) :: u,v,w,detJ
91
     
92
     LOGICAL :: stat, CSymmetry = .FALSE., VariableFlowWidth = .FALSE., &
93
            VariableLocalFlowWidth
94
       
95
     INTEGER, PARAMETER :: INDi(1:6) = (/ 1, 2, 3, 1, 2, 3 /) ,&
96
           INDj(1:6)=(/ 1, 2, 3, 2, 3, 1 /)
97

98
     INTEGER :: NewtonIter, NonlinearIter
99

100
     TYPE(Variable_t), POINTER :: AIFlowSol, TempSol, Var, FabricVariable
101
     TYPE(Variable_t), POINTER :: SpinVar 
102
     REAL(KIND=dp), POINTER :: SpinValues(:)
103
     INTEGER, POINTER :: SpinPerm(:)
104

105
     TYPE(Variable_t), POINTER :: DevStressVar 
106
     REAL(KIND=dp), POINTER :: DSValues(:)
107
     INTEGER, POINTER :: DSPerm(:)
108

109
     TYPE(Variable_t), POINTER :: StrainRateVar 
110
     REAL(KIND=dp), POINTER :: SRValues(:)
111
     INTEGER, POINTER :: SRPerm(:)
112

113
     REAL(KIND=dp), POINTER :: Temperature(:),AIFlow(:),Work(:,:), &
114
           ForceVector(:), VonMises(:), NodalAIFlow(:), AIFlowComp(:), &
115
           FabricValues(:)
116

117
     CHARACTER(LEN=MAX_NAME_LEN) :: EquationName
118

119
     INTEGER, POINTER :: TempPerm(:), AIFlowPerm(:), NodeIndexes(:), &
120
                         FabricPerm(:)
121

122
     INTEGER :: AIFlowType
123
     LOGICAL :: GotForceBC, GotIt, NewtonLinearization = .FALSE., &
124
                NormalTangential=.FALSE., UnFoundFatal=.TRUE.
125

126
     INTEGER :: body_id,bf_id
127
     INTEGER :: old_body = -1
128
     LOGICAL :: Isotropic, AllocationsDone = .FALSE., FreeSurface, &
129
                Requal0
130

131
     REAL(KIND=dp) :: FabricGrid(4878)
132
           
133
     REAL(KIND=dp), ALLOCATABLE:: LocalMassMatrix(:,:), &
134
       LocalStiffMatrix(:,:), LoadVector(:,:), LocalForce(:), &
135
       LocalTemperature(:), Alpha(:,:), Beta(:), & 
136
       ReferenceTemperature(:), BoundaryDispl(:), K1(:), K2(:), E1(:), &
137
       E2(:), E3(:), TimeForce(:), RefS(:), RefD(:), RefSpin(:), &
138
       LocalVelo(:,:), LocalFluidity(:), LocalFlowWidth(:)
139
            
140
     INTEGER :: NumberOfBoundaryNodes
141
     INTEGER, POINTER :: BoundaryReorder(:)
142

143
     REAL(KIND=dp) :: Bu, Bv, Bw, RM(3,3)
144
     REAL(KIND=dp), POINTER :: BoundaryNormals(:,:), &
145
         BoundaryTangent1(:,:), BoundaryTangent2(:,:)
146
     CHARACTER(LEN=MAX_NAME_LEN) :: viscosityFile
147
     REAL(KIND=dp) :: Radius
148
     REAL(KIND=dp) :: at, at0
149
!------------------------------------------------------------------------------
150
     SAVE NumberOfBoundaryNodes,BoundaryReorder,BoundaryNormals, &
151
              BoundaryTangent1, BoundaryTangent2, FabricGrid, viscosityFile
152

153
     SAVE TimeForce, Basis, dBasisdx, ddBasisddx
154
     SAVE LocalMassMatrix, LocalStiffMatrix, LoadVector, &
155
       LocalForce, ElementNodes, Alpha, Beta, LocalTemperature, LocalFlowWidth, &
156
       Isotropic,AllocationsDone,ReferenceTemperature,BoundaryDispl, &
157
       NodalAIFlow, K1, K2, E1, E2, E3, Wn, MinSRInvariant, old_body, &
158
       LocalFluidity
159

160
     SAVE RefD, RefS, RefSpin, LocalVelo, SlipCoeff
161
!------------------------------------------------------------------------------
162
!  Read constants from constants section of SIF file
163
!------------------------------------------------------------------------------
164
      Wn(7) = GetConstReal( Model % Constants, 'Gas Constant', GotIt )
165
      IF (.NOT.GotIt) THEN
166
        WRITE(Message,'(A)') 'VariableGas Constant  not found. &
167
                     &Setting to 8.314'
168
        CALL INFO('AIFlowSolve', Message, level=20)
169
        Wn(7) = 8.314
170
      ELSE
171
        WRITE(Message,'(A,F10.4)') 'Gas Constant = ',   Wn(7)
172
        CALL INFO('AIFlowSolve', Message , level = 20)
173
      END IF
174
!------------------------------------------------------------------------------
175
!    Get variables needed for solution
176
!------------------------------------------------------------------------------
177
      IF ( .NOT. ASSOCIATED( Solver % Matrix ) ) RETURN
178

179
      AIFlowSol => Solver % Variable
180
      AIFlowPerm => AIFlowSol % Perm
181
      STDOFs =  AIFlowSol % DOFs
182
      AIFlow => AIFlowSol % Values
183

184
      LocalNodes = COUNT( AIFlowPerm > 0 )
185
      IF ( LocalNodes <= 0 ) RETURN
186

187
      TempSol => VariableGet( Solver % Mesh % Variables, 'Temperature' )
188
      IF ( ASSOCIATED( TempSol) ) THEN
189
        TempPerm    => TempSol % Perm
190
        Temperature => TempSol % Values
191
      END IF
192

193
      FabricVariable => VariableGet(Solver % Mesh % Variables, 'Fabric')
194
      IF ( ASSOCIATED( FabricVariable ) ) THEN
195
       FabricPerm    => FabricVariable % Perm
196
       FabricValues => FabricVariable % Values
197
      END IF
198

199
      SpinVar => VariableGet(Solver % Mesh %Variables,'Spin')
200
      IF ( ASSOCIATED( SpinVar ) ) THEN
201
      SpinPerm => SpinVar % Perm    
202
      SpinValues => SpinVar % Values  
203
      END IF
204
      
205
      StrainRateVar => VariableGet(Solver % Mesh % Variables,'StrainRate')
206
      IF ( ASSOCIATED( StrainRateVar ) ) THEN
207
      SRPerm => StrainRateVar % Perm    
208
      SRValues => StrainRateVar % Values  
209
      END IF
210

211
      DevStressVar => &
212
               VariableGet(Solver % Mesh % Variables,'DeviatoricStress')
213
      IF ( ASSOCIATED( DevStressVar ) ) THEN
214
      DSPerm => DevStressVar % Perm    
215
      DSValues => DevStressVar % Values  
216
      END IF
217
      
218
      IF ( CurrentCoordinateSystem() == AxisSymmetric ) THEN
219
      CSymmetry = .TRUE.
220
      VariableFlowWidth = .TRUE.
221
      END IF
222

223
      StiffMatrix => Solver % Matrix
224
      ForceVector => StiffMatrix % RHS
225
      UNorm = Solver % Variable % Norm
226

227
!------------------------------------------------------------------------------
228
!     Allocate some permanent storage, this is done first time only
229
!------------------------------------------------------------------------------
230
      IF ( .NOT. AllocationsDone .OR. Solver % Mesh % Changed) THEN
231
        N = Model % MaxElementNodes
232
        dim = CoordinateSystemDimension()
233

234
       IF ( AllocationsDone ) THEN
235
         DEALLOCATE( ElementNodes % x,     &
236
                     ElementNodes % y,     &
237
                     ElementNodes % z,     &
238
                     BoundaryDispl,        &
239
                     ReferenceTemperature, &
240
                     LocalTemperature,     &
241
                     LocalFlowWidth,	   &
242
                     K1,K2,E1,E2,E3,       &
243
                     LocalForce,           &
244
                     RefD, RefS, RefSpin,  &
245
                     LocalVelo,            &
246
                     Basis,ddBasisddx,dBasisdx, &
247
                     TimeForce, &
248
                     LocalMassMatrix,      &
249
                     LocalStiffMatrix,     &
250
                     LoadVector, Alpha, Beta, &
251
                     SlipCoeff, LocalFluidity )
252
       END IF
253

254
       ALLOCATE( ElementNodes % x( N ), &
255
                 ElementNodes % y( N ), &
256
                 ElementNodes % z( N ), &
257
                 BoundaryDispl( N ), &
258
                 ReferenceTemperature( N ), &
259
                 LocalTemperature( N ), &
260
                 LocalFlowWidth (N), &
261
                 K1( N ), K2( N ), E1( N ), E2( N ), E3( N ), &
262
                 LocalForce( 2*STDOFs*N ),&
263
                 RefS(2*dim*LocalNodes ),&                              
264
                 RefD(2*dim*LocalNodes ),&                              
265
                 RefSpin((2*dim-3)*LocalNodes ),&                       
266
                 LocalVelo( 3,N ),&                                     
267
                 Basis( 2*N ),ddBasisddx(1,1,1), dBasisdx( 2*N,3 ), &
268
                 TimeForce( 2*STDOFs*N ), &
269
                 LocalMassMatrix( 2*STDOFs*N,2*STDOFs*N ),  &
270
                 LocalStiffMatrix( 2*STDOFs*N,2*STDOFs*N ),  &
271
                 LoadVector( 4,N ), Alpha( 3,N ), Beta( N ), &
272
                 SlipCoeff(3,N), LocalFluidity(N), STAT=istat )
273

274
       IF ( istat /= 0 ) THEN
275
          CALL Fatal( 'AIFlowSolve', 'Memory allocation error.' )
276
       END IF
277
!------------------------------------------------------------------------------
278

279
       AllocationsDone = .TRUE.
280
      END IF
281

282
!------------------------------------------------------------------------------
283
!    Do some additional initialization, and go for it
284
!------------------------------------------------------------------------------
285

286
!------------------------------------------------------------------------------
287
      NonlinearTol = GetConstReal( Solver % Values, &
288
        'Nonlinear System Convergence Tolerance' )
289

290
      NewtonTol = GetConstReal( Solver % Values, &
291
        'Nonlinear System Newton After Tolerance' )
292

293
      NewtonIter = GetInteger( Solver % Values, &
294
        'Nonlinear System Newton After Iterations' )
295

296
      NonlinearIter = GetInteger( Solver % Values, &
297
         'Nonlinear System Max Iterations',GotIt )
298

299
      IF ( .NOT.GotIt ) NonlinearIter = 1
300
!------------------------------------------------------------------------------
301

302
!------------------------------------------------------------------------------
303

304
      EquationName = GetString( Solver % Values, 'Equation' )
305

306
      FreeSurface = .FALSE.
307

308
!------------------------------------------------------------------------------
309
      DO iter=1,NonlinearIter
310

311
       at  = CPUTime()
312
       at0 = RealTime()
313

314
       CALL Info( 'AIFlowSolve', ' ', Level=4 )
315
       CALL Info( 'AIFlowSolve', ' ', Level=4 )
316
       CALL Info( 'AIFlowSolve', &
317
                   '-------------------------------------',Level=4 )
318
       WRITE( Message, * ) 'ANISOTROPIC FLOW SOLVER ITERATION', iter
319
       CALL Info( 'AIFlowSolve', Message,Level=4 )
320
       CALL Info( 'AIFlowSolve', &
321
                   '-------------------------------------',Level=4 )
322
       CALL Info( 'AIFlowSolve', ' ', Level=4 )
323
       CALL Info( 'AIFlowSolve', 'Starting assembly...',Level=4 )
324
!------------------------------------------------------------------------------
325
       CALL DefaultInitialize()
326
!------------------------------------------------------------------------------
327
       DO t=1,Solver % NumberOFActiveElements
328

329
         IF ( RealTime() - at0 > 1.0 ) THEN
330
           WRITE(Message,'(a,i3,a)' ) '   Assembly: ',  &
331
             INT(100.0 - 100.0 * (Solver % NumberOfActiveElements-t) / &
332
             (1.0*Solver % NumberOfActiveElements)), ' % done'
333
           CALL Info( 'AIFlowSolve', Message, Level=5 )
334
           at0 = RealTime()
335
         END IF
336

337
         CurrentElement => GetActiveElement(t)
338

339
         n = GetElementNOFNodes()
340
         NodeIndexes => CurrentElement % NodeIndexes
341

342
         ElementNodes % x(1:n) = Model % Nodes % x(NodeIndexes(1:n))
343
         ElementNodes % y(1:n) = Model % Nodes % y(NodeIndexes(1:n))
344
         ElementNodes % z(1:n) = Model % Nodes % z(NodeIndexes(1:n))
345

346
         Material => GetMaterial()
347

348
!------------------------------------------------------------------------------
349
!    Read in material constants from Material section
350
!------------------------------------------------------------------------------
351
         body_id = CurrentElement % BodyId
352
         IF (body_id /= old_body) Then 
353
            old_body = body_id
354
            Call  GetMaterialDefs()
355
        END IF
356

357
        LocalFluidity(1:n) = ListGetReal( Material, &
358
                         'Fluidity Parameter', n, NodeIndexes, GotIt,UnFoundFatal=UnFoundFatal)
359
           !Previous default value: LocalFluidity(1:n) = 1.0
360

361

362
       LocalFlowWidth(1:n) = ListGetReal ( Material, &
363
                        'FlowWidth', n, NodeIndexes, GotIt)
364
       IF (.NOT. GotIt) THEN
365
         IF (CSymmetry) THEN
366
           DO i=1,n
367
             LocalFlowWidth(i) = ElementNodes % x(i)
368
           END DO          
369
         END IF
370
       ELSE
371
         VariableFlowWidth = .TRUE.
372
       END IF
373
	  
374
! Test if flow width is locally constant (infinite radius case)
375
       VariableLocalFlowWidth = .TRUE.
376
       IF (MAXVAL(LocalFlowWidth(1:n))- &
377
                MINVAL(LocalFlowWidth(1:n)) == 0.0) &
378
                VariableLocalFlowWidth = .FALSE.
379
            
380
!------------------------------------------------------------------------------
381
!        Set body forces
382
!------------------------------------------------------------------------------
383
         LoadVector = 0.0D0
384

385
         BodyForce => GetBodyForce()
386
         IF ( ASSOCIATED( BodyForce ) ) THEN
387
           LoadVector(1,1:n) = LoadVector(1,1:n) + ListGetReal( &
388
                   BodyForce, 'AIFlow Force 1', n, NodeIndexes, gotIt)
389
           LoadVector(2,1:n) = LoadVector(2,1:n) + ListGetReal( &
390
                   BodyForce, 'AIFlow Force 2', n, NodeIndexes, gotIt)
391
           LoadVector(3,1:n) = LoadVector(3,1:n) + ListGetReal( & 
392
                   BodyForce, 'AIFlow Force 3', n, NodeIndexes, gotIt)
393
         END IF
394
!------------------------------------------------------------------------------
395
!        Get element local stiffness & mass matrices
396
!------------------------------------------------------------------------------
397
         LocalTemperature = 0.0D0
398
         IF ( ASSOCIATED(TempSol) ) THEN
399
           DO i=1,n
400
             k = TempPerm(NodeIndexes(i))
401
             LocalTemperature(i) = Temperature(k)
402
           END DO
403
         ELSE
404
           LocalTemperature(1:n) = 0.0d0
405
         END IF
406

407
! fabric not needed if isotropic
408
         IF(.NOT.Isotropic) THEN
409
           K1(1:n) = FabricValues( 5 * (FabricPerm(NodeIndexes(1:n))-1) + 1 )
410
           K2(1:n) = FabricValues( 5 * (FabricPerm(NodeIndexes(1:n))-1) + 2 )
411
           E1(1:n) = FabricValues( 5 * (FabricPerm(NodeIndexes(1:n))-1) + 3 )
412
           E2(1:n) = FabricValues( 5 * (FabricPerm(NodeIndexes(1:n))-1) + 4 )
413
           E3(1:n) = FabricValues( 5 * (FabricPerm(NodeIndexes(1:n))-1) + 5 )
414
         ENDIF
415

416
         LocalVelo = 0.0d0
417
         DO i=1,STDOFs - 1
418
            LocalVelo(i,1:n) = AIFlow( STDOFs*(AIFlowPerm(NodeIndexes(1:n))-1) + i)
419
         END DO
420

421
         CALL LocalMatrix( LocalMassMatrix, LocalStiffMatrix, &
422
              LocalForce, LoadVector, K1, K2, E1, E2, E3, LocalVelo, &
423
              LocalTemperature, LocalFlowWidth, LocalFluidity, CurrentElement, n, &
424
              ElementNodes, Wn, MinSRInvariant, Isotropic, VariableFlowWidth, &
425
              VariableLocalFlowWidth)
426

427
        TimeForce = 0.0d0
428
         CALL NSCondensate(N, N,STDOFs-1,LocalStiffMatrix,LocalForce,TimeForce )
429
!------------------------------------------------------------------------------
430
!        Update global matrices from local matrices 
431
!------------------------------------------------------------------------------
432
         CALL DefaultUpdateEquations( LocalStiffMatrix, LocalForce )
433
      END DO
434
      CALL DefaultFinishBulkAssembly()
435
      CALL Info( 'AIFlowSolve', 'Assembly done', Level=4 )
436

437
!------------------------------------------------------------------------------
438
!     Neumann & Newton boundary conditions
439
!------------------------------------------------------------------------------
440
      DO t = 1, Model % NumberOFBoundaryElements
441

442
        CurrentElement => GetBoundaryElement(t)
443
              
444
        IF ( GetElementFamily() == 101 ) CYCLE
445
        IF ( .NOT. ActiveBoundaryElement() ) CYCLE
446

447
        n = GetElementNOFNodes()
448
        NodeIndexes => CurrentElement % NodeIndexes
449

450
        ElementNodes % x(1:n) = Model % Nodes % x(NodeIndexes(1:n))
451
        ElementNodes % y(1:n) = Model % Nodes % y(NodeIndexes(1:n))
452
        ElementNodes % z(1:n) = Model % Nodes % z(NodeIndexes(1:n))
453

454
        BC => GetBC()
455
        IF ( ASSOCIATED( BC ) ) THEN
456
            LoadVector = 0.0D0
457
            Alpha      = 0.0D0
458
            Beta       = 0.0D0
459
!------------------------------------------------------------------------------
460
!           Force in given direction BC: \tau\cdot n = F
461
!------------------------------------------------------------------------------
462
            GotForceBC = .FALSE.
463

464
            LoadVector(1,1:n) = &
465
                    ListGetReal( BC, 'Force 1', n, NodeIndexes,  GotIt )
466
            GotForceBC = GotForceBC .OR. gotIt
467

468
            LoadVector(2,1:n) = & 
469
                     ListGetReal( BC, 'Force 2', n, NodeIndexes, GotIt )
470
            GotForceBC = GotForceBC .OR. gotIt
471

472
            LoadVector(3,1:n) = &
473
                     ListGetReal( BC, 'Force 3', n, NodeIndexes, GotIt )
474
            GotForceBC = GotForceBC .OR. gotIt
475

476
            Beta(1:n) = &
477
                ListGetReal( BC, 'Normal Force', n, NodeIndexes, GotIt )
478
            GotForceBC = GotForceBC .OR. gotIt
479

480
!------------------------------------------------------------------------------
481
!             slip boundary condition BC: \tau\cdot n = R_k u_k
482
!------------------------------------------------------------------------------
483

484
              SlipCoeff = 0.0d0
485
              SlipCoeff(1,1:n) =  ListGetReal( BC, &
486
                    'AIFlow Slip Coeff 1', n, NodeIndexes, GotIt )
487
              GotForceBC = GotForceBC .OR. gotIt
488

489
              SlipCoeff(2,1:n) =  ListGetReal( BC, &
490
                    'AIFlow Slip Coeff 2', n, NodeIndexes, GotIt )
491
              GotForceBC = GotForceBC .OR. gotIt
492

493
              SlipCoeff(3,1:n) =  ListGetReal( BC, &
494
                    'AIFlow Slip Coeff 3', n, NodeIndexes, GotIt )
495
              GotForceBC = GotForceBC .OR. gotIt
496

497
              NormalTangential = ListGetLogical( BC, &
498
                     'Normal-Tangential AIFlow', GotIt )
499
               
500
            IF ( .NOT.GotForceBC ) CYCLE
501

502
!------------------------------------------------------------------------------
503
            CALL LocalMatrixBoundary( LocalStiffMatrix, LocalForce, &
504
                 LoadVector, Alpha, Beta, SlipCoeff, NormalTangential, &
505
                 CurrentElement, n, ElementNodes, VariableFlowWidth, &
506
                 VariableLocalFlowWidth, LocalFlowWidth )
507
!------------------------------------------------------------------------------
508
!------------------------------------------------------------------------------
509
!           Update global matrices from local matrices (will also affect
510
!           LocalStiffMatrix and LocalForce if transientsimulation is on).
511
!------------------------------------------------------------------------------
512
            CALL DefaultUpdateEquations( LocalStiffMatrix, LocalForce )
513
!------------------------------------------------------------------------------
514
         END IF
515
      END DO
516
!------------------------------------------------------------------------------
517
      CALL DefaultFinishBoundaryAssembly()
518
      CALL DefaultFinishAssembly()
519

520
!------------------------------------------------------------------------------
521
!     Dirichlet boundary conditions
522
!------------------------------------------------------------------------------
523
      CALL DefaultDirichletBCs()
524
!------------------------------------------------------------------------------
525

526
      CALL Info( 'AIFlowSolve', 'Set boundaries done', Level=4 )
527

528
!------------------------------------------------------------------------------
529
!     Solve the system and check for convergence
530
!------------------------------------------------------------------------------
531
      PrevUNorm = UNorm
532

533
      UNorm = DefaultSolve()
534

535
      IF ( PrevUNorm + UNorm /= 0.0d0 ) THEN
536
         RelativeChange = 2.0d0 * ABS( PrevUNorm - UNorm) / ( PrevUnorm + UNorm)
537
      ELSE
538
         RelativeChange = 0.0d0
539
      END IF
540

541
      WRITE( Message, * ) 'Result Norm   : ',UNorm, PrevUNorm
542
      CALL Info( 'AIFlowSolve', Message, Level=4 )
543
      WRITE( Message, * ) 'Relative Change : ',RelativeChange
544
      CALL Info( 'AIFlowSolve', Message, Level=4 )
545

546
!------------------------------------------------------------------------------
547
      IF ( RelativeChange < NewtonTol .OR. &
548
             iter > NewtonIter ) NewtonLinearization = .TRUE.
549

550
      IF ( RelativeChange < NonLinearTol ) EXIT
551

552
!------------------------------------------------------------------------------
553
    END DO ! of nonlinear iter
554
!------------------------------------------------------------------------------
555

556
!------------------------------------------------------------------------------
557
!   Compute the StrainRate, Spin  and deviatoric Stress
558
!   Nodal values      
559
!------------------------------------------------------------------------------
560

561
     IF ((ASSOCIATED( StrainRateVar)).OR.(ASSOCIATED(DevStressVar))&
562
      .OR.(ASSOCIATED(SpinVar))) THEN
563
       RefD=0.
564
       RefS=0.
565
       RefSpin=0.
566
       IF (ASSOCIATED(StrainRateVar)) SRValues = 0.
567
       IF (ASSOCIATED(devStressVar)) DSValues = 0.
568
       IF (ASSOCIATED(SPinVar)) SpinValues = 0.
569

570
      DO t=1,Solver % NumberOFActiveElements
571

572
         CurrentElement => GetActiveElement(t)
573
         n = GetElementNOFNodes()
574
         NodeIndexes => CurrentElement % NodeIndexes
575

576
         body_id = CurrentElement % BodyId
577
         dim = CoordinateSystemDimension()
578

579
!------------------------------------------------------------------------------
580
!    Read in material constants from Material section
581
!------------------------------------------------------------------------------
582
        IF (body_id /= old_body) Then 
583
              old_body = body_id
584
              Call  GetMaterialDefs()
585
        END IF
586

587
        LocalFluidity(1:n) = ListGetReal( Material, &
588
                      'Fluidity Parameter', n, NodeIndexes, GotIt,&
589
                       UnFoundFatal=UnFoundFatal)
590
         !Previous default value: LocalFluidity(1:n) = 1.0
591

592
       LocalFlowWidth(1:n) = ListGetReal ( Material, &
593
                        'FlowWidth', n, NodeIndexes, GotIt)
594
       IF (.NOT. GotIt .AND. CSymmetry) THEN
595
          DO i=1,n
596
            LocalFlowWidth(i) = ElementNodes % x(i)
597
          END DO         
598
       END IF
599
                 
600
         ElementNodes % x(1:n) = Model % Nodes % x(NodeIndexes(1:n))
601
         ElementNodes % y(1:n) = Model % Nodes % y(NodeIndexes(1:n))
602
         ElementNodes % z(1:n) = Model % Nodes % z(NodeIndexes(1:n))
603

604
! n nodale values of the temperature
605
         
606
         LocalTemperature = 0.0D0
607
         IF ( ASSOCIATED(TempSol) ) THEN
608
           DO i=1,n
609
             k = TempPerm(NodeIndexes(i))
610
             LocalTemperature(i) = Temperature(k)
611
           END DO
612
         ELSE
613
           LocalTemperature(1:n) = 0.0d0
614
         END IF
615

616
! n nodales values of the 5 fabric parameters, not needed if isotropic
617
         IF(.NOT.Isotropic) Then
618
           K1(1:n) = FabricValues( 5 * (FabricPerm(NodeIndexes(1:n))-1) + 1 ) 
619
           K2(1:n) = FabricValues( 5 * (FabricPerm(NodeIndexes(1:n))-1) + 2 )
620
           E1(1:n) = FabricValues( 5 * (FabricPerm(NodeIndexes(1:n))-1) + 3 )
621
           E2(1:n) = FabricValues( 5 * (FabricPerm(NodeIndexes(1:n))-1) + 4 )
622
           E3(1:n) = FabricValues( 5 * (FabricPerm(NodeIndexes(1:n))-1) + 5 )
623
         END IF
624

625
! 2D U,V,p    STDOFs=3
626
! 3D U,V,W,p  STDOFs=4
627
         LocalVelo = 0.0d0
628
         DO i=1,STDOFs - 1
629
            LocalVelo(i,1:n) = AIFlow( STDOFs*(AIFlowPerm(NodeIndexes(1:n))-1) + i)
630
         END DO
631
         
632
! Go for all nodes of the element        
633
         Do i=1,n
634

635
! u, v, w local coord of node i
636
           u = CurrentElement % Type % NodeU(i)
637
           v = CurrentElement % Type % NodeV(i)
638
           w = CurrentElement % Type % NodeW(i)
639
       
640
            stat = ElementInfo(CurrentElement,ELementNodes,u,v,w,detJ, &
641
               Basis,dBasisdx,ddBasisddx,.FALSE.,.FALSE.)
642

643
! Variable flow width case when R=0 strain, stress not calculated exactly in
644
! x=0 (I agree it is not very nice, better solution ???)
645
       Requal0 = .False.
646
       IF (( VariableFlowWidth ) .And. & 
647
          (SUM(LocalFlowWidth (1:n) * Basis(1:n)) == 0.0)) THEN 
648
           Requal0 = .True.
649
            u= u + 0.0001 
650
            stat = ElementInfo(CurrentElement,ELementNodes,u,v,w,detJ, &
651
                     Basis,dBasisdx,ddBasisddx,.FALSE.,.FALSE.)
652
        END IF
653

654
           CALL LocalSD(NodalStresses, NodalStrainRate, NodalSpin, & 
655
                 LocalVelo, LocalTemperature, LocalFluidity,  &
656
                LocalFlowWidth, K1, K2, E1, E2, E3, Basis, dBasisdx, &
657
                CurrentElement, n, ElementNodes, dim, Wn, &
658
                MinSRInvariant, Isotropic, VariableFlowWidth, &
659
                VariableLocalFlowWidth)
660
                
661
        IF (Requal0)   NodalSpin = 0. 
662
        
663
           IF (ASSOCIATED(StrainRateVar)) &
664
             RefD(2*dim*(SRPerm(NodeIndexes(i))-1)+1 : &
665
                                      2*dim*SRPerm(NodeIndexes(i))) &
666
             =RefD(2*dim*(SRPerm(NodeIndexes(i))-1)+1 : &
667
                                      2*dim*SRPerm(NodeIndexes(i))) + 1.
668

669
          IF (ASSOCIATED(DevStressVar)) &
670
            RefS(2*dim*(DSPerm(NodeIndexes(i))-1)+1 : &
671
                                      2*dim*DSPerm(NodeIndexes(i))) &
672
            =RefS(2*dim*(DSPerm(NodeIndexes(i))-1)+1 :  &
673
                                      2*dim*DSPerm(NodeIndexes(i))) + 1.
674

675
          IF (ASSOCIATED(SpinVar)) &
676
            RefSpin((2*dim-3)*(SpinPerm(NodeIndexes(i))-1)+1 :  &
677
                                (2*dim-3)*SpinPerm(NodeIndexes(i))) &
678
            =RefSpin((2*dim-3)*(SpinPerm(NodeIndexes(i))-1)+1 :  &
679
                                (2*dim-3)*SpinPerm(NodeIndexes(i))) + 1.
680

681

682
           IF (ASSOCIATED(StrainRateVar)) THEN
683
             comp=0
684
             DO j=1,2*dim
685
               comp=comp+1
686
               SRValues(2*dim*(SRPerm(NodeIndexes(i))-1)+comp)=&
687
               SRValues(2*dim*(SRPerm(NodeIndexes(i))-1)+comp) + &
688
                NodalStrainRate(INDi(j),INDj(j))
689
             END DO
690
           END IF
691

692
           IF (ASSOCIATED(DevStressVar)) THEN
693
             comp=0
694
             DO j=1,2*dim
695
               comp=comp+1
696
               DSValues(2*dim*(DSPerm(NodeIndexes(i))-1)+comp)=&
697
                DSValues(2*dim*(DSPerm(NodeIndexes(i))-1)+comp) + &
698
                NodalStresses(INDi(j),INDj(j))
699
             END DO
700
           END IF
701

702
           IF (ASSOCIATED(SpinVar)) THEN
703
             comp=0
704
             DO j=1,(2*dim-3)
705
             comp=comp+1
706
             SpinValues((2*dim-3)*(SpinPerm(NodeIndexes(i))-1)+comp)=&
707
             SPinValues((2*dim-3)*(SpinPerm(NodeIndexes(i))-1)+comp) + &
708
             NodalSpin(INDi(j+3),INDj(j+3))
709
             END DO
710
           END IF
711

712
          END DO
713
        END DO
714

715
        IF (ASSOCIATED(StrainRateVar)) THEN
716
              WHERE(RefD > 0.)
717
                SRVAlues = SRValues / RefD
718
              END WHERE
719
        END IF
720
        
721
        IF (ASSOCIATED(DevStressVar)) THEN
722
           WHERE(RefS > 0.)
723
                DSVAlues = DSValues / RefS
724
           END WHERE
725
        END IF
726

727
        IF (ASSOCIATED(SpinVar)) THEN
728
            WHERE(RefSpin > 0.)
729
                SpinVAlues = SpinValues / RefSpin
730
            END WHERE
731
        END IF
732
        
733
    END IF
734

735
!------------------------------------------------------------------------------
736
!  END  Compute the StrainRate and Deviatoric Stress
737
!------------------------------------------------------------------------------
738
      
739
CONTAINS
740

741
      SUBROUTINE GetMaterialDefs()
742
      ! check if we are isotropic or not
743
      Isotropic = ListGetLogical( Material , 'Isotropic',Gotit )
744
      IF (.NOT.Gotit) Then
745
          Isotropic = .False.
746
           WRITE(Message,'(A)') 'Isotropic set to False'
747
           CALL INFO('AIFlowSolve', Message, Level = 20)
748
      ELSE
749
           IF ( (ASSOCIATED( FabricVariable )).AND.Isotropic ) Then
750
              WRITE(Message,'(A)') 'Be careful Isotropic is true &
751
                           & and Fabric is defined!'
752
              CALL INFO('AIFlowSolve', Message, Level = 1)
753
           END IF
754
      END IF
755

756
      IF (.NOT.Isotropic) Then
757
        ! Get the viscosity file and store the viscosities into FabricGrid
758
         viscosityFile = ListGetString( Material ,'Viscosity File',GotIt,UnFoundFatal )
759
         OPEN( 1, File = viscosityFile)
760
         DO i=1,813
761
             READ( 1, '(6(e14.8))' ) FabricGrid( 6*(i-1)+1:6*(i-1)+6 )
762
         END DO
763
         CLOSE(1)
764
      ENDIF
765

766
      Wn(2) = ListGetConstReal( Material , 'Powerlaw Exponent', GotIt,UnFoundFatal=UnFoundFatal)
767
         !Previous default value: Wn(2) = 1.0
768
      WRITE(Message,'(A,F10.4)') 'Powerlaw Exponent = ',   Wn(2)
769
      CALL INFO('AIFlowSolve', Message, Level = 20)
770

771
      Wn(3) = ListGetConstReal( Material, 'Activation Energy 1', GotIt,UnFoundFatal=UnFoundFatal)
772
         !Previous default value: Wn(3) = 1.0
773
      WRITE(Message,'(A,F10.4)') 'Activation Energy 1 = ',   Wn(3)
774
      CALL INFO('AIFlowSolve', Message, Level = 20)
775

776
      Wn(4) = ListGetConstReal( Material, 'Activation Energy 2', GotIt,UnFoundFatal=UnFoundFatal)
777
         !Previous default value: Wn(4) = 1.0
778
      WRITE(Message,'(A,F10.4)') 'Activation Energy 2 = ',   Wn(4)
779
      CALL INFO('AIFlowSolve', Message, Level = 20)
780

781
      Wn(5) = ListGetConstReal(Material, 'Reference Temperature', GotIt,UnFoundFatal=UnFoundFatal)
782
         !Previous default value: Wn(5) = -10.0 (Celsius)
783
      WRITE(Message,'(A,F10.4)') 'Reference Temperature = ',   Wn(5)
784
      CALL INFO('AIFlowSolve', Message, Level = 20)
785
      
786
      Wn(6) = ListGetConstReal( Material, 'Limit Temperature', GotIt,UnFoundFatal=UnFoundFatal)
787
         !Previous default value: Wn(6) = -10.0 (Celsius)
788
      WRITE(Message,'(A,F10.4)') 'Limit Temperature = ',   Wn(6)
789
      CALL INFO('AIFlowSolve', Message, Level = 20)
790

791
! Get the Minimum value of the Effective Strain rate 
792
      MinSRInvariant = 100.0*AEPS
793

794
      IF ( Wn(2) > 1.0  ) THEN
795
        MinSRInvariant =  &
796
             ListGetConstReal( Material, 'Min Second Invariant', GotIt )
797
        IF (.NOT.GotIt) THEN
798
          WRITE(Message,'(A)') 'Variable Min Second Invariant not &
799
                    &found. Setting to 100.0*AEPS )'
800
          CALL INFO('AIFlowSolve', Message, Level = 20)
801
        ELSE
802
          WRITE(Message,'(A,E14.8)') 'Min Second Invariant = ', MinSRInvariant
803
          CALL INFO('AIFlowSolve', Message, Level = 20)
804
        END IF
805
      END IF
806

807
!------------------------------------------------------------------------------
808
      END SUBROUTINE GetMaterialDefs
809
!------------------------------------------------------------------------------
810

811
!------------------------------------------------------------------------------
812
      SUBROUTINE LocalMatrix( MassMatrix, StiffMatrix, ForceVector, &
813
              LoadVector, NodalK1, NodalK2, NodalEuler1, NodalEuler2, &
814
              NodalEuler3, NodalVelo, NodalTemperature, NodalFlowWidth, &
815
              NodalFluidity, Element, n, Nodes, Wn, MinSRInvariant, Isotropic, &
816
              VariableFlowWidth, VariableLocalFlowWidth )
817
                       
818
!------------------------------------------------------------------------------
819

820
     REAL(KIND=dp) :: StiffMatrix(:,:), MassMatrix(:,:)
821
     REAL(KIND=dp) :: LoadVector(:,:), NodalVelo(:,:)
822
     REAL(KIND=dp) :: Wn(7), MinSRInvariant
823
     REAL(KIND=dp), DIMENSION(:) :: ForceVector, NodalK1, NodalK2, &
824
             NodalEuler1, NodalEuler2, NodalEuler3, NodalTemperature, &
825
             NodalFluidity, NodalFlowWidth
826
     TYPE(Nodes_t) :: Nodes
827
     TYPE(Element_t) :: Element
828
     LOGICAL :: Isotropic, VariableFlowWidth, VariableLocalFlowWidth
829
     INTEGER :: n
830
!------------------------------------------------------------------------------
831
!
832
     REAL(KIND=dp) :: Basis(2*n),ddBasisddx(1,1,1)
833
     REAL(KIND=dp) :: dBasisdx(2*n,3),SqrtElementMetric
834

835
     REAL(KIND=dp) :: Force(3), K1, K2, Euler1, Euler2, Euler3
836
     Real(kind=dp) :: Bg, BGlenT
837

838
     REAL(KIND=dp), DIMENSION(4,4) :: A,M
839
     REAL(KIND=dp) :: Load(3),Temperature,  C(6,6)
840
     REAL(KIND=dp) :: nn, ss,pp, LGrad(3,3), SR(3,3), Stress(3,3), D(6), epsi
841
	 INTEGER :: INDi(6),INDj(6)
842

843

844
     INTEGER :: i, j, k, p, q, t, dim, NBasis, ind(3)
845

846
     REAL(KIND=dp) :: s,u,v,w, Radius, B(6,3), G(3,6), FW
847
  
848
     REAL(KIND=dp) :: dDispldx(3,3), ai(3), Angle(3), a2(6)
849
     TYPE(GaussIntegrationPoints_t), TARGET :: IntegStuff
850

851
     INTEGER :: N_Integ
852

853
     REAL(KIND=dp), DIMENSION(:), POINTER :: U_Integ,V_Integ,W_Integ,S_Integ
854

855
     LOGICAL :: stat
856

857
     INTERFACE
858
      Subroutine R2Ro(a2,dim,ai,angle)
859
         USE Types
860
         REAL(KIND=dp),intent(in) :: a2(6)
861
         Integer :: dim
862
         REAL(KIND=dp),intent(out) :: ai(3), Angle(3)
863
      End Subroutine R2Ro
864
                 
865
      Subroutine OPILGGE_ai_nl(ai,Angle,etaI,eta36)
866
          USE Types
867
          REAL(kind=dp), INTENT(in),  DIMENSION(3)   :: ai
868
          REAL(kind=dp), INTENT(in),  DIMENSION(3)   :: Angle
869
          REAL(kind=dp), INTENT(in),  DIMENSION(:)   :: etaI
870
          REAL(kind=dp), INTENT(out), DIMENSION(6,6) :: eta36
871
        END SUBROUTINE OPILGGE_ai_nl
872

873
      END INTERFACE
874
!------------------------------------------------------------------------------
875
      dim = CoordinateSystemDimension()
876

877
      ForceVector = 0.0D0
878
      StiffMatrix = 0.0D0
879
      MassMatrix  = 0.0D0
880

881
!    
882
!    Integration stuff
883
!    
884
      NBasis = 2*n
885
      IntegStuff = GaussPoints( Element, Element % Type % GaussPoints2 )
886

887
      U_Integ => IntegStuff % u
888
      V_Integ => IntegStuff % v
889
      W_Integ => IntegStuff % w
890
      S_Integ => IntegStuff % s
891
      N_Integ =  IntegStuff % n
892
!
893
!   Now we start integrating
894
!
895
      DO t=1,N_Integ
896

897
      u = U_Integ(t)
898
      v = V_Integ(t)
899
      w = W_Integ(t)
900

901
!------------------------------------------------------------------------------
902
!     Basis function values & derivatives at the integration point
903
!------------------------------------------------------------------------------
904
      stat = ElementInfo( Element,Nodes,u,v,w,SqrtElementMetric, &
905
            Basis,dBasisdx,ddBasisddx,.FALSE.,.TRUE. )
906

907
      s = SqrtElementMetric * S_Integ(t)
908
!------------------------------------------------------------------------------
909
!  
910
!     Force at integration point
911
!   
912
      Force = 0.0d0
913
      DO i=1,dim
914
         Force(i) = SUM( LoadVector(i,1:n)*Basis(1:n))
915
      END DO
916
!
917
!     Temperature at the integration point
918
!
919
      Temperature = SUM( NodalTemperature(1:n)*Basis(1:n) )
920
      Wn(1) = SUM( NodalFluidity(1:n)*Basis(1:n) )
921
      Bg=BGlenT(Temperature,Wn)
922
      ss=1.0_dp
923

924
! if not isotropic use GOLF
925
      C = 0.0_dp
926
      IF (.NOT.Isotropic) Then
927
         a2(1) = SUM( NodalK1(1:n) * Basis(1:n) ) 
928
         a2(2) = SUM( NodalK2(1:n) * Basis(1:n) ) 
929
         a2(3) = 1.d0 - a2(1) - a2(2)
930
         a2(4) = SUM( NodalEuler1(1:n) * Basis(1:n) )
931
         a2(5) = SUM( NodalEuler2(1:n) * Basis(1:n) )
932
         a2(6) = SUM( NodalEuler3(1:n) * Basis(1:n) )
933
      
934
         CALL R2Ro(a2,dim,ai,angle)
935
         CALL OPILGGE_ai_nl(ai,Angle,FabricGrid,C)
936
! else use isotropic law
937
      ELSE
938
          Do i=1,3
939
            C(i,i)=2.0_dp 
940
          End do
941
          Do i=4,6
942
            C(i,i)=1.0_dp
943
          End do
944
      ENDIF
945
      
946
      FW = SUM( NodalFlowWidth(1:n) * Basis(1:n) )
947
      Radius = FW / (SUM( NodalFlowWidth(1:n) * dBasisdx(1:n,1)) )
948
      IF (.NOT. VariableLocalFlowWidth) Radius = 10e7
949
      IF (VariableFlowWidth )  s = s * FW 
950

951
!
952
! Case non-linear
953
! -----------------------------
954

955
      IF ( Wn(2) > 1.0 ) THEN
956

957
         Bg=Bg**(1.0/Wn(2))
958

959
        LGrad = MATMUL( NodalVelo(:,1:n), dBasisdx(1:n,:) )
960
        SR = 0.5 * ( LGrad + TRANSPOSE(LGrad) )
961
      
962
        IF ( VariableFlowWidth ) THEN 
963
          SR(1,3) = 0.0
964
          SR(2,3) = 0.0
965
          SR(3,1) = 0.0
966
          SR(3,2) = 0.0
967
          SR(3,3) = 0.0
968
          IF ( Radius > 10*AEPS ) THEN
969
            SR(3,3) = SUM( Nodalvelo(1,1:n) * Basis(1:n) ) /Radius
970
          END IF
971
          epsi = SR(1,1)+SR(2,2)+SR(3,3)
972
          DO i=1,3   
973
            SR(i,i) = SR(i,i) - epsi/3.0
974
          END DO
975
        ELSE
976
          epsi = SR(1,1)+SR(2,2)+SR(3,3)
977
          DO i=1,dim 
978
            SR(i,i) = SR(i,i) - epsi/dim
979
          END DO
980
        END IF
981

982
    ! Compute the invariant 
983
        nn = (1.0 - Wn(2))/(2.0*Wn(2))
984
        ss = 0.0_dp
985
        DO i = 1, 3
986
          DO j = 1, 3
987
            ss = ss + SR(i,j)**2.
988
          END DO
989
        END DO
990
        IF (ss < MinSRInvariant ) ss = MinSRInvariant
991
        ss = (2.*ss)**nn
992
      END IF
993

994
! Non relative viscosity matrix
995
       C = C * ss/Bg
996

997
!
998
!    Loop over basis functions (of both unknowns and weights)
999
!
1000
      A = 0.0d0
1001
      M = 0.0d0
1002
      B = 0.0d0
1003

1004
      DO p=1,NBasis
1005

1006
       G = 0.0d0
1007

1008
       IF ( VariableFlowWidth ) THEN
1009
          G(1,1) = dBasisdx(p,1)
1010
          G(1,3) = Basis(p) / Radius
1011
          G(1,4) = dBasisdx(p,2)
1012
          G(2,2) = dBasisdx(p,2)
1013
          G(2,4) = dBasisdx(p,1)
1014
       ELSE
1015
          G(1,1) = dBasisdx(p,1)
1016
          G(2,2) = dBasisdx(p,2)
1017
          G(3,3) = dBasisdx(p,3)
1018
          G(1,4) = dBasisdx(p,2)
1019
          G(2,4) = dBasisdx(p,1)
1020
          G(2,5) = dBasisdx(p,3)
1021
          G(3,5) = dBasisdx(p,2)
1022
          G(1,6) = dBasisdx(p,3)
1023
          G(3,6) = dBasisdx(p,1)
1024
       END IF
1025

1026
       G = MATMUL( G, C )
1027
 
1028
       DO q=1,NBasis
1029

1030
         B = 0.0d0
1031

1032
         IF ( VariableFlowWidth ) THEN
1033
            B(1,1) = dBasisdx(q,1)
1034
            B(2,2) = dBasisdx(q,2)
1035
            B(3,1) = Basis(q) / Radius
1036
            B(4,1) = dBasisdx(q,2)
1037
            B(4,2) = dBasisdx(q,1)
1038
         ELSE
1039
            B(1,1) = dBasisdx(q,1)
1040
            B(2,2) = dBasisdx(q,2)
1041
            B(3,3) = dBasisdx(q,3)
1042
            B(4,1) = dBasisdx(q,2)
1043
            B(4,2) = dBasisdx(q,1)
1044
            B(5,2) = dBasisdx(q,3)
1045
            B(5,3) = dBasisdx(q,2)
1046
            B(6,1) = dBasisdx(q,3)
1047
            B(6,3) = dBasisdx(q,1)
1048
         END IF
1049

1050
         A(1:3,1:3) = MATMUL( G, B )
1051

1052
! Pressure gradient
1053
         DO i=1,dim
1054
            A(i,dim+1) = -dBasisdx(p,i) * Basis(q)
1055
         END DO
1056
         IF ( VariableFlowWidth ) A(1,dim+1) =  A(1,dim+1) - Basis(p) * Basis(q) / Radius
1057

1058
! Continuity equation:
1059
         DO i=1,dim
1060
            A(dim+1,i) = dBasisdx(q,i) * Basis(p)
1061
         END DO
1062
         IF ( VariableFlowWidth ) A(dim+1,1) =  A(dim+1,1) + Basis(p) * Basis(q) / Radius
1063
         A(dim+1, dim+1) = 0.0d0
1064

1065
! Add nodal matrix to element matrix
1066
         DO i=1,dim+1
1067
            DO j=1,dim+1
1068
               StiffMatrix( (dim+1)*(p-1)+i,(dim+1)*(q-1)+j ) =  &
1069
                    StiffMatrix( (dim+1)*(p-1)+i,(dim+1)*(q-1)+j ) + s*A(i,j)
1070
            END DO
1071
         END DO
1072

1073
       END DO
1074

1075
! The righthand side...
1076
        Load = 0.0d0
1077
  
1078
        DO i=1,dim
1079
           Load(i) = Load(i) + Force(i) * Basis(p)
1080
        END DO
1081

1082
        DO i=1,dim
1083
           ForceVector((dim+1)*(p-1)+i) = ForceVector((dim+1)*(p-1)+i) + s*Load(i)
1084
        END DO
1085
      END DO
1086

1087
      END DO 
1088
!------------------------------------------------------------------------------
1089
      END SUBROUTINE LocalMatrix
1090
!------------------------------------------------------------------------------
1091

1092
!------------------------------------------------------------------------------
1093
      SUBROUTINE LocalMatrixBoundary( BoundaryMatrix, BoundaryVector, &
1094
                 LoadVector, NodalAlpha, NodalBeta, NodalSlipCoeff, & 
1095
                 NormalTangential, Element, n, Nodes, &
1096
                 VariableFlowWidth, VariableLocalFlowWidth, NodalFlowWidth )
1097
                      
1098
!------------------------------------------------------------------------------
1099
     REAL(KIND=dp) :: BoundaryMatrix(:,:),BoundaryVector(:)
1100
     REAL(KIND=dp) :: NodalAlpha(:,:),NodalBeta(:),LoadVector(:,:)
1101
     REAL(KIND=dp) :: NodalSlipCoeff(:,:), NodalFlowWidth(:)
1102
     TYPE(Element_t),POINTER  :: Element
1103
     TYPE(Nodes_t)    :: Nodes
1104
     LOGICAL :: NormalTangential
1105
     INTEGER, POINTER :: NodeIndexes(:)
1106
     INTEGER :: n
1107
!------------------------------------------------------------------------------
1108
     REAL(KIND=dp) :: Basis(n),ddBasisddx(1,1,1)
1109
     REAL(KIND=dp) :: dBasisdx(n,3),SqrtElementMetric, FW
1110

1111
     REAL(KIND=dp) :: u,v,w,s
1112
     REAL(KIND=dp) :: Force(3),Alpha(3),Beta,Normal(3)
1113
     REAL(KIND=dp), POINTER :: U_Integ(:),V_Integ(:),W_Integ(:),S_Integ(:)
1114
     REAL(KIND=dp) :: Tangent(3),Tangent2(3),Vect(3), SlipCoeff
1115
     REAL(KIND=dp) :: Up,Vp,Wp
1116
     INTEGER :: i,t,q,p,dim,N_Integ, c
1117

1118
     LOGICAL :: stat, VariableFlowWidth, VariableLocalFlowWidth
1119
     LOGICAL,SAVE :: AllocationDone=.False.
1120
    
1121
     TYPE(GaussIntegrationPoints_t), TARGET :: IntegStuff
1122
!------------------------------------------------------------------------------
1123

1124
      dim = CoordinateSystemDimension()
1125
      c=dim+1
1126

1127
      BoundaryVector = 0.0D0
1128
      BoundaryMatrix = 0.0D0
1129
!
1130
!  Integration stuff
1131
!
1132
      IntegStuff = GaussPoints( element )
1133
      U_Integ => IntegStuff % u
1134
      V_Integ => IntegStuff % v
1135
      W_Integ => IntegStuff % w
1136
      S_Integ => IntegStuff % s
1137
      N_Integ =  IntegStuff % n
1138

1139
    NodeIndexes => Element % NodeIndexes
1140
                   
1141
     NodalFlowWidth(1:n) = ListGetReal ( Material, &
1142
                        'FlowWidth', n, NodeIndexes, GotIt)
1143
       IF (.NOT. GotIt .AND. CSymmetry) THEN
1144
           DO i=1,n
1145
             NodalFlowWidth(i) = Nodes % x(i)
1146
           END DO          
1147
       END IF
1148
!
1149
!  Now we start integrating
1150
!
1151
      DO t=1,N_Integ
1152

1153
       u = U_Integ(t)
1154
       v = V_Integ(t)
1155
       w = W_Integ(t)
1156

1157
!------------------------------------------------------------------------------
1158
!    Basis function values & derivatives at the integration point
1159
!------------------------------------------------------------------------------
1160
       stat = ElementInfo( Element, Nodes, u, v, w, SqrtElementMetric, &
1161
                Basis, dBasisdx, ddBasisddx, .FALSE. )
1162

1163
       FW = SUM( NodalFlowWidth(1:n) * Basis(1:n) )
1164

1165
       s = SqrtElementMetric * S_Integ(t)
1166

1167
       IF (.NOT. VariableLocalFlowWidth ) Radius = 10e7
1168
       IF ( VariableFlowWidth )    s= s * FW
1169

1170
!------------------------------------------------------------------------------
1171
       Force = 0.0D0
1172
       DO i=1,dim
1173
         Force(i) = SUM( LoadVector(i,1:n)*Basis(1:n) )
1174
         Alpha(i) = SUM( NodalAlpha(i,1:n)*Basis(1:n) )
1175
       END DO
1176

1177
       Normal = NormalVector( Element,Nodes,u,v,.TRUE. )
1178
       Force = Force + SUM( NodalBeta(1:n)*Basis(1:n) ) * Normal
1179

1180
       SELECT CASE( Element % TYPE % DIMENSION )
1181
       CASE(1)
1182
        Tangent(1) =  Normal(2)
1183
        Tangent(2) = -Normal(1)
1184
        Tangent(3) =  0.0d0
1185
       CASE(2)
1186
        CALL TangentDirections( Normal, Tangent, Tangent2 ) 
1187
       END SELECT
1188
  
1189
       IF ( ANY( NodalSlipCoeff(:,:) /= 0.0d0 ) ) THEN
1190
         DO p=1,n
1191
           DO q=1,n
1192
             DO i=1,DIM
1193
              SlipCoeff = SUM( NodalSlipCoeff(i,1:n) * Basis(1:n) )
1194
  
1195
              IF (NormalTangential ) THEN
1196
                SELECT CASE(i)
1197
                   CASE(1)
1198
                     Vect = Normal
1199
                   CASE(2)
1200
                     Vect = Tangent
1201
                   CASE(3)
1202
                     Vect = Tangent2
1203
                END SELECT
1204
  
1205
                DO j=1,DIM
1206
                   DO k=1,DIM
1207
                      BoundaryMatrix( (p-1)*c+j,(q-1)*c+k ) = &
1208
                         BoundaryMatrix( (p-1)*c+j,(q-1)*c+k ) + &
1209
                          s * SlipCoeff * Basis(q) * Basis(p) * Vect(j) * Vect(k)
1210
                   END DO
1211
                END DO
1212
               ELSE
1213
                 BoundaryMatrix( (p-1)*c+i,(q-1)*c+i ) = &
1214
                     BoundaryMatrix( (p-1)*c+i,(q-1)*c+i ) + &
1215
                          s * SlipCoeff * Basis(q) * Basis(p)
1216
               END IF
1217
             END DO
1218
           END DO
1219
         END DO
1220
       END IF
1221

1222

1223

1224
      DO p=1,N
1225
       DO q=1,N
1226
         DO i=1,dim
1227
           BoundaryMatrix((p-1)*(dim+1)+i,(q-1)*(dim+1)+i) =  &
1228
             BoundaryMatrix((p-1)*(dim+1)+i,(q-1)*(dim+1)+i) + &
1229
               s * Alpha(i) * Basis(q) * Basis(p)
1230
         END DO
1231
       END DO
1232
      END DO
1233

1234
      DO q=1,N
1235
       DO i=1,dim
1236
         BoundaryVector((q-1)*(dim+1)+i) = BoundaryVector((q-1)*(dim+1)+i) + &
1237
                   s * Basis(q) * Force(i)
1238
       END DO
1239
      END DO
1240

1241
      END DO
1242
!------------------------------------------------------------------------------
1243
      END SUBROUTINE LocalMatrixBoundary
1244
!------------------------------------------------------------------------------
1245

1246

1247
!------------------------------------------------------------------------------
1248
      SUBROUTINE LocalSD( Stress, StrainRate, Spin, &
1249
        NodalVelo, NodalTemp, NodalFluidity, NodalFlowWidth, &
1250
        NodalK1, NodalK2, NodalE1, NodalE2, NodalE3, &
1251
        Basis, dBasisdx, Element, n,  Nodes, dim,  Wn, MinSRInvariant, &
1252
        Isotropic, VariableFlowWidth, VariableLocalFlowWidth )
1253
!------------------------------------------------------------------------------
1254
!    Subroutine to compute the nodal Strain-Rate, Stress, ...
1255
!------------------------------------------------------------------------------
1256
     INTEGER :: n, dim
1257
     INTEGER :: INDi(6),INDj(6)
1258
     REAL(KIND=dp) :: Stress(:,:), StrainRate(:,:), Spin(:,:)
1259
     REAL(KIND=dp) :: NodalVelo(:,:), NodalTemp(:), NodalFluidity(:), &
1260
                      NodalFlowWidth(:)
1261
     REAL(KIND=dp) :: Basis(2*n), ddBasisddx(1,1,1)
1262
     REAL(KIND=dp) :: dBasisdx(2*n,3)
1263
     REAL(KIND=dp) :: detJ
1264
     REAL(KIND=dp) :: NodalK1(:), NodalK2(:)
1265
     REAL(KIND=dp) :: NodalE1(:), NodalE2(:), NodalE3(:)
1266
     REAL(KIND=dp) :: u, v, w      
1267
     REAL(KIND=dp) :: Wn(7),  D(6), MinSRInvariant
1268
     LOGICAL :: Isotropic,VariableFlowWidth, VariableLocalFlowWidth
1269
      
1270
     TYPE(Nodes_t) :: Nodes
1271
     TYPE(Element_t) :: Element
1272
!------------------------------------------------------------------------------
1273
     LOGICAL :: stat
1274
     INTEGER :: i,j,k,p,q
1275
     REAL(KIND=dp) :: LGrad(3,3), Radius, Temp, ai(3), Angle(3),a2(6)
1276
     REAL(KIND=dp) :: C(6,6), epsi
1277
     Real(kind=dp) :: Bg, BGlenT, ss, nn
1278
!------------------------------------------------------------------------------
1279
     INTERFACE
1280
      Subroutine R2Ro(a2,dim,ai,angle)
1281
         USE Types
1282
         REAL(KIND=dp),intent(in) :: a2(6)
1283
         Integer :: dim
1284
         REAL(KIND=dp),intent(out) :: ai(3), Angle(3)
1285
      End Subroutine R2Ro
1286
                 
1287
      Subroutine OPILGGE_ai_nl(ai,Angle,etaI,eta36)
1288
          USE Types
1289
          REAL(kind=dp), INTENT(in),  DIMENSION(3)   :: ai
1290
          REAL(kind=dp), INTENT(in),  DIMENSION(3)   :: Angle
1291
          REAL(kind=dp), INTENT(in),  DIMENSION(:)   :: etaI
1292
          REAL(kind=dp), INTENT(out), DIMENSION(6,6) :: eta36
1293
        END SUBROUTINE OPILGGE_ai_nl
1294
      END INTERFACE
1295
!------------------------------------------------------------------------------
1296
     
1297
!
1298
!     Temperature at the integration point
1299
      Temp = SUM( NodalTemp(1:n)*Basis(1:n) )
1300
      Wn(1) = SUM( NodalFluidity(1:n)*Basis(1:n) )
1301

1302

1303
      Stress = 0.0
1304
      StrainRate = 0.0
1305
      Spin = 0.0
1306
!
1307
!    Compute strainRate : 
1308
!    -------------------
1309

1310
      LGrad = MATMUL( NodalVelo(:,1:n), dBasisdx(1:n,:) )
1311
      
1312
      StrainRate = 0.5 * ( LGrad + TRANSPOSE(LGrad) )
1313

1314
      IF ( VariableFlowWidth ) THEN
1315

1316
        StrainRate(1,3) = 0.0
1317
        StrainRate(2,3) = 0.0
1318
        StrainRate(3,1) = 0.0
1319
        StrainRate(3,2) = 0.0
1320
        StrainRate(3,3) = 0.0
1321
        
1322
        IF (SUM( NodalFlowWidth(1:n) * dBasisdx(1:n,1)) == 0) THEN
1323
              Radius = 10e7
1324
        ELSE
1325
                Radius = SUM( NodalFlowWidth(1:n) * Basis(1:n) ) / &
1326
                (SUM( NodalFlowWidth(1:n) * dBasisdx(1:n,1)) )
1327
        END IF
1328
        
1329
        IF ( Radius > 10*AEPS ) THEN
1330
         StrainRate(3,3) = SUM( Nodalvelo(1,1:n) * Basis(1:n) ) / Radius
1331
        END IF
1332
        epsi = StrainRate(1,1)+StrainRate(2,2)+StrainRate(3,3)
1333
        DO i=1,3   
1334
          StrainRate(i,i) = StrainRate(i,i) - epsi/3.0
1335
        END DO
1336
      ELSE
1337
        epsi = StrainRate(1,1)+StrainRate(2,2)+StrainRate(3,3)
1338
        DO i=1,dim 
1339
          StrainRate(i,i) = StrainRate(i,i) - epsi/dim
1340
        END DO
1341
      END IF
1342

1343
!
1344
!    Compute Spin : 
1345
!    --------------
1346

1347
        Spin = 0.5 * ( LGrad - TRANSPOSE(LGrad) )
1348

1349
!
1350
!    Compute deviatoric stresses: 
1351
!    ----------------------------
1352

1353
      IF (.Not.Isotropic) then
1354
	    C = 0.0_dp
1355
!    Material parameters at that point
1356
        a2(1) = SUM( NodalK1(1:n) * Basis(1:n) ) 
1357
        a2(2) = SUM( NodalK2(1:n) * Basis(1:n) ) 
1358
        a2(3) = 1.d0 - a2(1) - a2(2)
1359
        a2(4) = SUM( NodalE1(1:n) * Basis(1:n) )
1360
        a2(5) = SUM( NodalE2(1:n) * Basis(1:n) )
1361
        a2(6) = SUM( NodalE3(1:n) * Basis(1:n) )
1362
      
1363
        CALL R2Ro(a2,dim,ai,Angle)
1364
        CALL OPILGGE_ai_nl(ai,Angle,FabricGrid,C)
1365
         
1366
!
1367
!    Compute deviatoric stresses: 
1368
!    ----------------------------
1369
      D(1) = StrainRate(1,1)
1370
      D(2) = StrainRate(2,2)
1371
      D(3) = StrainRate(3,3)
1372
      D(4) = 2. * StrainRate(1,2)
1373
      D(5) = 2. * StrainRate(2,3)
1374
      D(6) = 2. * StrainRate(3,1)
1375
      
1376
      INDi(1:6) = (/ 1, 2, 3, 1, 2, 3 /)
1377
      INDj(1:6) = (/ 1, 2, 3, 2, 3, 1 /)
1378
      DO k = 1, 2*dim
1379
       DO j = 1, 2*dim
1380
        Stress( INDi(k),INDj(k) ) = &
1381
        Stress( INDi(k),INDj(k) ) + C(k,j) * D(j)
1382
       END DO
1383
       IF (k > 3)  Stress( INDj(k),INDi(k) ) = Stress( INDi(k),INDj(k) )
1384
      END DO
1385

1386
	ELSE  ! ISOTROPIC CASE
1387
	     Stress=2._dp * StrainRate
1388
	END IF
1389
	   
1390

1391
! non relative viscosities
1392
    ! Glen fluidity       
1393
	 Bg=BGlenT(Temp,Wn)
1394
	 ss=1.0_dp
1395
     ! Case Non linear
1396
	 IF (Wn(2) > 1.0) THEN 
1397
	   Bg=Bg**(1.0/Wn(2))
1398
       
1399
        ss = 0.0_dp
1400
        DO i = 1, 3
1401
          DO j = 1, 3
1402
            ss = ss + StrainRate(i,j)**2
1403
          END DO
1404
        END DO
1405
        nn = (1.0 - Wn(2))/(2.0*Wn(2))
1406
        ss = (2.*ss)**nn
1407
         
1408
        IF (ss < MinSRInvariant ) ss = MinSRInvariant
1409
        
1410
		
1411
      END IF
1412
       
1413
	   Stress=Stress*ss/Bg
1414

1415

1416

1417
!------------------------------------------------------------------------------
1418
      END SUBROUTINE LocalSD      
1419
!------------------------------------------------------------------------------
1420
!        
1421
!------------------------------------------------------------------------------
1422
      END SUBROUTINE AIFlowSolver_nlD2
1423
!------------------------------------------------------------------------------
1424

1425
Product

Resources

Company
