← Back to branch summary

~ubuntu-branches/ubuntu/maverick/elmerfem/maverick

« back to all changes in this revision

Viewing changes to fem/src/LevelSet.src

  • Committer: Bazaar Package Importer
  • Author(s): Adam C. Powell, IV
  • Date: 2009-04-08 16:58:46 UTC
  • Revision ID: james.westby@ubuntu.com-20090408165846-qydtxka1cvrjw2ix
Tags: upstream-5.5.0.svn.4096.dfsg
Import upstream version 5.5.0.svn.4096.dfsg

Show diffs side-by-side

added added

removed removed

Lines of Context:
fem/src/LevelSet.src
 
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 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
! *
 
26
! *  Subroutines for moving and renormalizing the level set function, 
 
27
! *  and computing the curvature of the level set field. 
 
28
! *
 
29
! ******************************************************************************
 
30
! *
 
31
! *  Authors: Juha Ruokolainen, Peter R�back
 
32
! *  Email:   Peter.Raback@csc.fi
 
33
! *  Web:     http://www.csc.fi/elmer
 
34
! *  Address: CSC - IT Center for Science Ltd.
 
35
! *           Keilaranta 14
 
36
! *           02101 Espoo, Finland 
 
37
! *
 
38
! *  Original Date: 16.11.2005
 
39
! *
 
40
! *****************************************************************************/
 
41
 
 
42
 
 
43
 
 
44
!------------------------------------------------------------------------------
 
45
   SUBROUTINE LevelSetSolver( Model,Solver,Timestep,TransientSimulation )
 
46
!------------------------------------------------------------------------------
 
47
!******************************************************************************
 
48
!
 
49
!  Solve the advection equation for the Levelset-function using stabilization
 
50
!  or bubbless. For the accuracy it is advisable to use 2nd order time-stepping
 
51
!  and courant number smaller than one. 
 
52
!
 
53
!  ARGUMENTS:
 
54
!
 
55
!  TYPE(Model_t) :: Model,  
 
56
!     INPUT: All model information (mesh, materials, BCs, etc...)
 
57
!
 
58
!  TYPE(Solver_t) :: Solver
 
59
!     INPUT: Linear equation solver options
 
60
!
 
61
!  REAL(KIND=dp) :: Timestep,
 
62
!     INPUT: Timestep size for time dependent simulations
 
63
!
 
64
!  LOGICAL :: TransientSimulation
 
65
!     INPUT: Steady state or transient simulation
 
66
!
 
67
!******************************************************************************
 
68
     USE Types
 
69
     USE DefUtils
 
70
     USE SolverUtils
 
71
     USE MaterialModels
 
72
     USE Integration
 
73
 
 
74
     IMPLICIT NONE
 
75
!------------------------------------------------------------------------------
 
76
 
 
77
     TYPE(Model_t), TARGET :: Model
 
78
     TYPE(Solver_t) :: Solver
 
79
     REAL(KIND=dp) :: Timestep
 
80
     LOGICAL :: TransientSimulation
 
81
 
 
82
!------------------------------------------------------------------------------
 
83
!    Local variables
 
84
!------------------------------------------------------------------------------
 
85
     INTEGER :: i,j,k,n,t,iter,istat,bf_id,CoordinateSystem
 
86
 
 
87
     TYPE(Matrix_t),POINTER  :: StiffMatrix
 
88
     TYPE(Nodes_t)   :: ElementNodes
 
89
     TYPE(Element_t),POINTER :: CurrentElement
 
90
     TYPE(ValueList_t), POINTER :: Material
 
91
 
 
92
     REAL(KIND=dp) :: Norm,RelativeChange
 
93
     LOGICAL :: Stabilize = .FALSE., GotIt, AllocationsDone = .FALSE.
 
94
     INTEGER, POINTER :: SurfPerm(:), NodeIndexes(:)
 
95
     REAL(KIND=dp), POINTER :: Surface(:), ForceVector(:), Surf(:), PrevSurface(:) 
 
96
     REAL(KIND=dp), ALLOCATABLE :: LocalMassMatrix(:,:),&
 
97
       LocalStiffMatrix(:,:),Load(:),LocalForce(:),TimeForce(:)
 
98
     INTEGER :: NSDOFs,NonlinearIter,body_id,mat_id,dim
 
99
     REAL(KIND=dp) :: NonlinearTol,dt, r, ct, DsMax
 
100
     REAL(KIND=dp), ALLOCATABLE :: ElemVelo(:,:), SurfaceFlux(:)
 
101
 
 
102
     SAVE LocalMassMatrix,LocalStiffMatrix,Load, &
 
103
         TimeForce, LocalForce, ElementNodes,AllocationsDone, &
 
104
         ElemVelo, Surf, SurfaceFlux
 
105
 
 
106
     REAL(KIND=dp) :: at,totat,st,totst,CPUTime
 
107
 
 
108
 
 
109
!------------------------------------------------------------------------------
 
110
!    Get variables needed for solution
 
111
!------------------------------------------------------------------------------
 
112
 
 
113
     CALL Info( 'LevelSetSolver','-------------------------------------', Level=4 )
 
114
     CALL Info( 'LevelSetSolver', 'Solving for the levelset function', Level=4 )
 
115
     CALL Info( 'LevelSetSolver','-------------------------------------', Level=4 )
 
116
 
 
117
 
 
118
     IF ( .NOT. ASSOCIATED( Solver % Matrix ) ) RETURN
 
119
 
 
120
     CoordinateSystem = CurrentCoordinateSystem()
 
121
     dim = CoordinateSystemDimension()
 
122
 
 
123
     Surface => Solver % Variable % Values
 
124
     IF ( SIZE( Surface ) == 0 ) RETURN
 
125
     SurfPerm => Solver % Variable % Perm
 
126
     PrevSurface => Solver % Variable % PrevValues(:,1)
 
127
 
 
128
     StiffMatrix => Solver % Matrix
 
129
     ForceVector => StiffMatrix % RHS
 
130
     Norm = Solver % Variable % Norm
 
131
 
 
132
!------------------------------------------------------------------------------
 
133
!    Allocate some permanent storage, this is done first time only
 
134
!------------------------------------------------------------------------------
 
135
     IF ( .NOT. AllocationsDone ) THEN
 
136
       N = Solver % Mesh % MaxElementNodes
 
137
 
 
138
       ALLOCATE( ElemVelo(3, N), &
 
139
           ElementNodes % x( N ),   &
 
140
           ElementNodes % y( N ),   &
 
141
           ElementNodes % z( N ),   &
 
142
           LocalForce( 2*N ),       &
 
143
           TimeForce( 2*N ),        &
 
144
           LocalMassMatrix( 2*N,2*N ),  &
 
145
           LocalStiffMatrix( 2*N,2*N ), &
 
146
           Surf( N ), &
 
147
           SurfaceFlux( N ), &
 
148
           STAT=istat )
 
149
 
 
150
       IF ( istat /= 0 ) THEN
 
151
         CALL Fatal( 'LevelSetSolver', 'Memory allocation error.' )
 
152
       END IF
 
153
 
 
154
       AllocationsDone = .TRUE.
 
155
     END IF
 
156
!------------------------------------------------------------------------------
 
157
!    Do some additional initialization, and go for it
 
158
!------------------------------------------------------------------------------
 
159
 
 
160
     Stabilize = ListGetLogical( Solver % Values,'Stabilize',GotIt )
 
161
 
 
162
     NonlinearTol = ListGetConstReal( Solver % Values, &
 
163
           'Nonlinear System Convergence Tolerance',GotIt )
 
164
 
 
165
     NonlinearIter = ListGetInteger( Solver % Values, &
 
166
        'Nonlinear System Max Iterations',GotIt )
 
167
     IF ( .NOT.GotIt ) NonlinearIter = 1
 
168
 
 
169
     dt = Timestep
 
170
 
 
171
!------------------------------------------------------------------------------
 
172
 
 
173
     totat = 0.0d0
 
174
     totst = 0.0d0
 
175
 
 
176
     CALL Info( 'LevelSetSolver','-------------------------------------', Level=4 )
 
177
     CALL Info( 'LevelSetSolver', 'Updating the levelset function', Level=4 )
 
178
     CALL Info( 'LevelSetSolver','-------------------------------------', Level=4 )
 
179
 
 
180
     DO iter = 1, NonlinearIter
 
181
       
 
182
       at = CPUTime()
 
183
 
 
184
!------------------------------------------------------------------------------
 
185
       CALL InitializeToZero( StiffMatrix, ForceVector )
 
186
!------------------------------------------------------------------------------
 
187
 
 
188
!------------------------------------------------------------------------------
 
189
!      Do the assembly for bulk elements
 
190
!------------------------------------------------------------------------------
 
191
 
 
192
       DO t=1,Solver % NumberOfActiveElements
 
193
 
 
194
         CurrentElement => Solver % Mesh % Elements(Solver % ActiveElements(t))
 
195
         body_id = CurrentElement % Bodyid    
 
196
         mat_id = ListGetInteger( Model % Bodies( body_id ) % Values, 'Material' )
 
197
         Material => Model % Materials(mat_id) % Values
 
198
 
 
199
!------------------------------------------------------------------------------
 
200
!        Set the current element pointer in the model structure to
 
201
!        reflect the element being processed
 
202
!------------------------------------------------------------------------------
 
203
         Model % CurrentElement => CurrentElement
 
204
         n = CurrentElement % TYPE % NumberOfNodes
 
205
         NodeIndexes => CurrentElement % NodeIndexes
 
206
 
 
207
!------------------------------------------------------------------------------
 
208
!        Get element nodal coordinates
 
209
!------------------------------------------------------------------------------
 
210
         ElementNodes % x(1:n) = Solver % Mesh % Nodes % x(NodeIndexes)
 
211
         ElementNodes % y(1:n) = Solver % Mesh % Nodes % y(NodeIndexes)
 
212
         ElementNodes % z(1:n) = Solver % Mesh % Nodes % z(NodeIndexes)
 
213
 
 
214
         Surf = Surface( SurfPerm(NodeIndexes) )
 
215
 
 
216
!------------------------------------------------------------------------------
 
217
!         Computed velocity field
 
218
!------------------------------------------------------------------------------
 
219
 
 
220
         ElemVelo(1,1:n) = ListGetReal( Material,'Levelset Velocity 1',n,NodeIndexes,GotIt)
 
221
         ElemVelo(2,1:n) = ListGetReal( Material,'Levelset Velocity 2',n,NodeIndexes,GotIt)
 
222
         IF(dim == 3) ElemVelo(3,1:n) = ListGetReal( Material,'Levelset Velocity 3',n,NodeIndexes,GotIt)
 
223
 
 
224
!------------------------------------------------------------------------------
 
225
!          Given surface flux
 
226
!------------------------------------------------------------------------------
 
227
         bf_id = ListGetInteger( Model % Bodies(CurrentElement % BodyId) % &
 
228
             Values, 'Body Force', GotIt )
 
229
         IF ( GotIt ) THEN
 
230
           SurfaceFlux(1:n) = ListGetReal( Model % BodyForces(bf_id) % Values,  &
 
231
               'Levelset Flux',n,NodeIndexes,gotIt )
 
232
         ELSE           
 
233
           SurfaceFlux(1:n) = 0.0d0
 
234
         END IF
 
235
 
 
236
!------------------------------------------------------------------------------
 
237
         CALL LocalMatrix( LocalMassMatrix, LocalStiffMatrix, LocalForce, Surf, &
 
238
             SurfaceFlux, ElemVelo, Stabilize, CurrentElement, n, ElementNodes )
 
239
       
 
240
!------------------------------------------------------------------------------
 
241
!        If time dependent simulation add mass matrix to stiff matrix
 
242
!------------------------------------------------------------------------------
 
243
         TimeForce = 0.0_dp
 
244
         IF ( TransientSimulation ) THEN
 
245
!------------------------------------------------------------------------------
 
246
!          NOTE: This will replace LocalStiffMatrix and LocalForce with the
 
247
!                combined information...
 
248
!------------------------------------------------------------------------------
 
249
           CALL Add1stOrderTime( LocalMassMatrix, LocalStiffMatrix, &
 
250
               LocalForce, dt, n, 1, SurfPerm(NodeIndexes), Solver )
 
251
         END IF
 
252
!------------------------------------------------------------------------------
 
253
!      Update global matrix and rhs vector from local matrix & vector
 
254
!------------------------------------------------------------------------------
 
255
         IF ( .NOT.Stabilize ) THEN
 
256
           CALL Condensate( N, LocalStiffMatrix,  LocalForce, TimeForce )
 
257
         END IF
 
258
         
 
259
         CALL UpdateGlobalEquations( StiffMatrix, LocalStiffMatrix, &
 
260
             ForceVector, LocalForce, n, 1, SurfPerm(NodeIndexes) )
 
261
         
 
262
       END DO     !  Of ActiveElements
 
263
 
 
264
 
 
265
!------------------------------------------------------------------------------
 
266
!    FinishAssemebly must be called after all other assembly steps, but before
 
267
!    Dirichlet boundary settings. Actually no need to call it except for
 
268
!    transient simulations.
 
269
!------------------------------------------------------------------------------
 
270
       CALL FinishAssembly( Solver,ForceVector )
 
271
!------------------------------------------------------------------------------
 
272
!    Dirichlet boundary conditions
 
273
!------------------------------------------------------------------------------
 
274
       CALL SetDirichletBoundaries( Model, StiffMatrix, ForceVector, &
 
275
          ComponentName(Solver % Variable),1,1,SurfPerm )
 
276
!------------------------------------------------------------------------------
 
277
     CALL Info( 'LevelSetSolver', 'Assembly done', Level=4 )
 
278
     at = CPUTime() - at
 
279
     totat = totat + at
 
280
    
 
281
!------------------------------------------------------------------------------
 
282
!     Solve the system and check for convergence
 
283
!------------------------------------------------------------------------------
 
284
     st = CPUTime()
 
285
 
 
286
     CALL SolveSystem( StiffMatrix, ParMatrix, ForceVector, &
 
287
         Surface, Norm, 1, Solver )
 
288
     RelativeChange = Solver % Variable % NonlinChange
 
289
 
 
290
     st = CPUTIme()-st
 
291
     totst = totst + st
 
292
      
 
293
     IF(NonlinearIter > 1) THEN
 
294
       WRITE( Message, * ) 'Iteration   : ',iter
 
295
       CALL Info( 'LevelSetSolver', Message, Level=4 )      
 
296
     END IF
 
297
     WRITE( Message, * ) 'Result Norm   : ',Norm
 
298
     CALL Info( 'LevelSetSolver', Message, Level=4 )
 
299
     WRITE( Message, * ) 'Relative Change : ',RelativeChange
 
300
     CALL Info( 'LevelSetSolver', Message, Level=4 )
 
301
     
 
302
     IF(RelativeChange < NonlinearTol) EXIT
 
303
 
 
304
!------------------------------------------------------------------------------
 
305
   END DO ! of the nonlinear iteration
 
306
!------------------------------------------------------------------------------
 
307
 
 
308
   WRITE(Message,'(a,F8.2)') 'Assembly done in time (s):',totat
 
309
   CALL Info( 'LevelsetSolver',Message, Level=4 )
 
310
   
 
311
   WRITE(Message,'(a,F8.2)') 'Solution done in time (s):',totst
 
312
   CALL Info( 'LevelsetSolver',Message, Level=4 )
 
313
 
 
314
   DsMax = MAXVAL( ABS(Surface - PrevSurface) )
 
315
   WRITE(Message,'(a,ES12.3)') 'Maximum Levelset Change',dsmax
 
316
   CALL Info( 'LevelSetSolver',Message, Level=4 )     
 
317
   CALL ListAddConstReal(Model % Simulation,'res: LevelSet Max Change',dsmax)
 
318
 
 
319
 
 
320
!------------------------------------------------------------------------------
 
321
 
 
322
CONTAINS
 
323
 
 
324
   SUBROUTINE LocalMatrix( MassMatrix,StiffMatrix,ForceVector,  &
 
325
      Surf, Flux, NodalVelo, Stabilize,Element,n,Nodes )
 
326
!------------------------------------------------------------------------------
 
327
 
 
328
     REAL(KIND=dp), DIMENSION(:)   :: ForceVector,Surf,Flux
 
329
     REAL(KIND=dp), DIMENSION(:,:) :: MassMatrix,StiffMatrix,NodalVelo
 
330
     LOGICAL :: Stabilize
 
331
     INTEGER :: n
 
332
     TYPE(Nodes_t) :: Nodes
 
333
     TYPE(Element_t), POINTER :: Element
 
334
 
 
335
!------------------------------------------------------------------------------
 
336
!    Local variables
 
337
!------------------------------------------------------------------------------
 
338
!
 
339
     REAL(KIND=dp) :: ddBasisddx(n,3,3)
 
340
     REAL(KIND=dp) :: Basis(2*n)
 
341
     REAL(KIND=dp) :: dBasisdx(2*n,3),SQRTElementMetric
 
342
     REAL(KIND=dp) :: Velo(3),Force,Grad(3),GradAbs
 
343
     REAL(KIND=dp) :: A,M,Load, FL
 
344
     REAL(KIND=dp) :: VNorm,hK,mK
 
345
     REAL(KIND=dp) :: Lambda=1.0,Pe,Pe1,Pe2,Tau,x,y,z
 
346
 
 
347
     INTEGER :: i,j,k,c,p,q,t,dim,N_Integ,NBasis
 
348
     REAL(KIND=dp) :: s,u,v,w
 
349
     TYPE(GaussIntegrationPoints_t), TARGET :: IntegStuff
 
350
     REAL(KIND=dp) :: SU(n),SW(n)
 
351
     REAL(KIND=dp), DIMENSION(:), POINTER :: U_Integ,V_Integ,W_Integ,S_Integ
 
352
 
 
353
     LOGICAL :: stat,Bubbles
 
354
 
 
355
!------------------------------------------------------------------------------
 
356
 
 
357
     dim = CoordinateSystemDimension()
 
358
     c = dim + 1
 
359
 
 
360
     ForceVector = 0.0D0
 
361
     StiffMatrix = 0.0D0
 
362
     MassMatrix  = 0.0D0
 
363
     Load = 0.0d0
 
364
     Grad = 0.0d0
 
365
 
 
366
     NBasis = n
 
367
     Bubbles = .FALSE.
 
368
     IF ( .NOT. Stabilize ) THEN
 
369
       NBasis = 2*n
 
370
       Bubbles = .TRUE.
 
371
     END IF
 
372
 
 
373
!------------------------------------------------------------------------------
 
374
!    Integration stuff
 
375
!------------------------------------------------------------------------------
 
376
     IF ( Bubbles ) THEN
 
377
        IntegStuff = GaussPoints( element, Element % TYPE % GaussPoints2 )
 
378
     ELSE
 
379
        IntegStuff = GaussPoints( element )
 
380
     END IF
 
381
     U_Integ => IntegStuff % u
 
382
     V_Integ => IntegStuff % v
 
383
     W_Integ => IntegStuff % w
 
384
     S_Integ => IntegStuff % s
 
385
     N_Integ =  IntegStuff % n
 
386
 
 
387
!------------------------------------------------------------------------------
 
388
!    Stabilization parameters: hK, mK (take a look at Franca et.al.)
 
389
!    If there is no convection term we don t need stabilization.
 
390
!------------------------------------------------------------------------------
 
391
     IF ( Stabilize ) THEN
 
392
       hK = element % hK
 
393
       mK = element % StabilizationMK
 
394
     END IF
 
395
 
 
396
!------------------------------------------------------------------------------
 
397
!    Now we start integrating
 
398
!------------------------------------------------------------------------------
 
399
     DO t=1,N_Integ
 
400
 
 
401
       u = U_Integ(t)
 
402
       v = V_Integ(t)
 
403
       w = W_Integ(t)
 
404
 
 
405
!------------------------------------------------------------------------------
 
406
!      Basis function values & derivatives at the integration point
 
407
!------------------------------------------------------------------------------
 
408
       stat = ElementInfo( Element,Nodes,u,v,w,SQRTElementMetric, &
 
409
           Basis,dBasisdx,ddBasisddx,Stabilize,Bubbles )
 
410
 
 
411
       s = SQRTElementMetric * S_Integ(t)
 
412
!------------------------------------------------------------------------------
 
413
!      Coordinatesystem dependent info
 
414
!------------------------------------------------------------------------------
 
415
       IF ( Coordinates == AxisSymmetric .OR. &
 
416
           Coordinates == CylindricSymmetric ) THEN
 
417
         s = s * SUM( Nodes % x(1:n)*Basis(1:n) )
 
418
       END IF
 
419
 
 
420
!------------------------------------------------------------------------------
 
421
         
 
422
       FL = SUM( Flux(1:n) * Basis(1:n) )
 
423
       
 
424
       DO i=1,dim
 
425
         Grad(i) = SUM( dBasisdx(1:n,i) * Surf(1:n) )
 
426
       END DO
 
427
       GradAbs = SQRT( SUM( Grad(1:dim) * Grad(1:dim) ) )
 
428
       IF ( GradAbs > 10*AEPS ) THEN
 
429
         Grad = Grad / GradAbs
 
430
       END IF
 
431
 
 
432
!------------------------------------------------------------------------------
 
433
!      Velocity from previous iteration at the integration point
 
434
!------------------------------------------------------------------------------
 
435
 
 
436
       Velo = 0.0D0
 
437
       Velo(1) = SUM( NodalVelo(1,1:n)*Basis(1:n) )
 
438
       Velo(2) = SUM( NodalVelo(2,1:n)*Basis(1:n) )
 
439
       IF ( dim > 2 ) Velo(3) = SUM( NodalVelo(3,1:n)*Basis(1:n) )
 
440
 
 
441
 
 
442
       IF ( Stabilize ) THEN
 
443
!------------------------------------------------------------------------------
 
444
!           Stabilization parameter Tau
 
445
!------------------------------------------------------------------------------
 
446
         VNorm = SQRT( SUM(Velo(1:dim)**2) ) + FL
 
447
         Pe  = 1.0d0
 
448
         Tau = 0.0D0
 
449
         IF ( VNorm /= 0.0 ) THEN
 
450
           Tau = hK * Pe / (2 * VNorm)
 
451
           Tau = 1.0d0 / SQRT( (2.0d0 / dt)**2 + 1.0d0/Tau**2 )
 
452
         ELSE
 
453
           Tau = dt / 2.0d0 
 
454
         END IF
 
455
 
 
456
!------------------------------------------------------------------------------
 
457
!           Compute residual & stablization vectors
 
458
!------------------------------------------------------------------------------
 
459
         DO p=1,N
 
460
           SU(p) = 0.0d0
 
461
           DO i = 1,dim
 
462
             SU(p) = SU(p) + dBasisdx(p,i) * (Velo(i) + FL * Grad(i) )
 
463
           END DO
 
464
 
 
465
           SW(p) = 0.0d0
 
466
           DO i = 1,dim
 
467
             SW(p) = SW(p) + dBasisdx(p,i) * (Velo(i) + FL * Grad(i) )
 
468
           END DO
 
469
         END DO
 
470
       END IF
 
471
 
 
472
 
 
473
!------------------------------------------------------------------------------
 
474
!      Loop over basis functions of both unknowns and weights
 
475
!------------------------------------------------------------------------------
 
476
       DO p=1,NBasis
 
477
         DO q=1,NBasis
 
478
!------------------------------------------------------------------------------
 
479
!         The diffusion term
 
480
!------------------------------------------------------------------------------
 
481
 
 
482
           M = Basis(q) * Basis(p)
 
483
           A = 0.0d0
 
484
 
 
485
           IF(GradAbs > AEPS) THEN
 
486
             DO i=1,dim
 
487
               A = A + FL * (Grad(i) / GradAbs) * dBasisdx(q,i) * Basis(p)
 
488
             END DO
 
489
           END IF
 
490
 
 
491
!------------------------------------------------------------------------------
 
492
!           The convection term
 
493
!------------------------------------------------------------------------------
 
494
 
 
495
           DO i=1,dim
 
496
             A = A + Velo(i) * dBasisdx(q,i) * Basis(p)
 
497
           END DO
 
498
!------------------------------------------------------------------------------
 
499
!           Next we add the stabilization...
 
500
!------------------------------------------------------------------------------
 
501
           IF ( Stabilize ) THEN
 
502
             A = A + Tau * SU(q) * SW(p)
 
503
             M = M + Tau * Basis(q) * SW(p)
 
504
           END IF
 
505
 
 
506
           StiffMatrix(p,q) = StiffMatrix(p,q) + s * A
 
507
           MassMatrix(p,q)  = MassMatrix(p,q)  + s * M
 
508
 
 
509
         END DO
 
510
       END DO
 
511
 
 
512
!------------------------------------------------------------------------------
 
513
!      The righthand side...
 
514
!------------------------------------------------------------------------------
 
515
       Force = 0.0d0
 
516
 
 
517
       DO p=1,NBasis
 
518
          Load = Basis(p)
 
519
          IF ( Stabilize ) Load = Load + Tau * SW(p)
 
520
          ForceVector(p) = ForceVector(p) + s * Force * Load
 
521
       END DO
 
522
     END DO
 
523
 
 
524
!------------------------------------------------------------------------------
 
525
   END SUBROUTINE LocalMatrix
 
526
!------------------------------------------------------------------------------
 
527
 
 
528
 
 
529
!------------------------------------------------------------------------------
 
530
 END SUBROUTINE LevelSetSolver
 
531
!------------------------------------------------------------------------------
 
532
 
 
533
 
 
534
 
 
535
 
 
536
!------------------------------------------------------------------------------
 
537
   SUBROUTINE LevelSetDistance( Model,Solver,Timestep,TransientSimulation )
 
538
!------------------------------------------------------------------------------
 
539
!******************************************************************************
 
540
!
 
541
!  Renormalizes the levelset function using straight-forward geometric search
 
542
!  Also includes an option to do the covection at the same time as an alternavtive
 
543
!  for using a separate solver for the advection.
 
544
!
 
545
!******************************************************************************
 
546
 
 
547
     USE SolverUtils
 
548
     USE MaterialModels
 
549
     USE Integration
 
550
 
 
551
     IMPLICIT NONE
 
552
!------------------------------------------------------------------------------
 
553
 
 
554
     TYPE(Model_t), TARGET :: Model
 
555
     TYPE(Solver_t) :: Solver 
 
556
     REAL(KIND=dp) :: Timestep
 
557
     LOGICAL :: TransientSimulation
 
558
 
 
559
!------------------------------------------------------------------------------
 
560
!    Local variables
 
561
!------------------------------------------------------------------------------
 
562
     TYPE(Nodes_t)   :: ElementNodes
 
563
     TYPE(Element_t),POINTER :: CurrentElement, Element
 
564
     TYPE(ValueList_t), POINTER :: Material
 
565
     TYPE(Variable_t), POINTER :: SurfSol, DistanceSol
 
566
 
 
567
     INTEGER :: i,j,k,l,n,t,iter,istat,body_id,mat_id,bf_id,&
 
568
          CoordinateSystem, TimesVisited = 0
 
569
     REAL(KIND=dp) :: Norm,x,y,s,a,b,c,d,x0,x1,y0,y1
 
570
     INTEGER, POINTER :: NodeIndexes(:)
 
571
     LOGICAL :: GotIt, Convect, ExtractAllocated = .FALSE., DistanceAllocated=.FALSE.
 
572
     INTEGER, POINTER :: SurfPerm(:)
 
573
     REAL(KIND=dp), POINTER :: Surface(:),Distance(:), Surf(:)
 
574
     REAL(KIND=dp), ALLOCATABLE :: ZeroNodes(:,:,:), Direction(:)
 
575
     INTEGER, POINTER :: DistancePerm(:)
 
576
     INTEGER :: ZeroLevels, ReinitializeInterval, ExtractInterval
 
577
     LOGICAL :: Reinitialize, Extrct
 
578
     REAL(KIND=dp) :: Relax, dt, r, NarrowBand, DsMax
 
579
     REAL(KIND=dp), ALLOCATABLE :: ElemVelo(:,:), SurfaceFlux(:)
 
580
     REAL(KIND=dp) :: at,totat,st,totst,CPUTime
 
581
     CHARACTER(LEN=MAX_NAME_LEN) :: LevelSetVariableName
 
582
 
 
583
     SAVE ElementNodes, ElemVelo, Direction, ZeroNodes, TimesVisited, &
 
584
         Distance, DistancePerm, ExtractAllocated, DistanceAllocated
 
585
 
 
586
!------------------------------------------------------------------------------
 
587
!    Get variables needed for solution
 
588
!------------------------------------------------------------------------------
 
589
 
 
590
     TimesVisited = TimesVisited + 1
 
591
     ReinitializeInterval = ListGetInteger(Solver % Values,&
 
592
          'Reinitialize Interval',GotIt) 
 
593
     IF(.NOT. GotIt) ReinitializeInterval = 1
 
594
 
 
595
     ExtractInterval = ListGetInteger(Solver % Values,&
 
596
         'Extract Interval',GotIt) 
 
597
     IF(.NOT. GotIt) ExtractInterval = ReinitializeInterval
 
598
     
 
599
     IF( ReinitializeInterval == 0) THEN
 
600
       Reinitialize = .FALSE.
 
601
     ELSE       
 
602
       Reinitialize = ( MOD(TimesVisited, ReinitializeInterval) == 0 )
 
603
     END IF
 
604
 
 
605
     IF( ExtractInterval == 0) THEN
 
606
       Extrct = Reinitialize
 
607
     ELSE
 
608
       Extrct = Reinitialize .OR. ( MOD(TimesVisited, ExtractInterval) == 0 )
 
609
     END IF
 
610
     
 
611
     IF(.NOT. Extrct) THEN
 
612
       CALL Info( 'LevelSetDistance','--------------------------------------', Level=4 )
 
613
       CALL Info( 'LevelSetDistance','Doing nothing this time', Level=4 )
 
614
       CALL Info( 'LevelSetDistance','--------------------------------------', Level=4 )          
 
615
       RETURN
 
616
     END IF
 
617
 
 
618
     ! The variable that should be reinitialized
 
619
     LevelSetVariableName = ListGetString(Solver % Values,'LevelSet Variable',GotIt) 
 
620
     IF(.NOT. GotIT) LevelSetVariableName = 'Surface'
 
621
     SurfSol => VariableGet( Solver % Mesh % Variables, TRIM(LevelSetVariableName) )
 
622
     IF(ASSOCIATED(SurfSol)) THEN
 
623
       Surface  => SurfSol % Values
 
624
       SurfPerm => SurfSol % Perm
 
625
     ELSE
 
626
       CALL Warn('LevelSetDistance','SurfSol does not exist: '//TRIM(LevelSetVariableName))
 
627
       RETURN
 
628
     END IF
 
629
 
 
630
     CoordinateSystem = CurrentCoordinateSystem() 
 
631
     Convect = ListGetLogical(Solver % Values,'Levelset Convect',GotIt)
 
632
     NarrowBand = ListGetConstReal(Solver % Values,'Narrow Band',GotIt)
 
633
     IF(.NOT. GotIt) NarrowBand = HUGE(NarrowBand)
 
634
     dsMax = 0.0d0
 
635
     dt = Timestep
 
636
     
 
637
 
 
638
!------------------------------------------------------------------------------
 
639
!    Allocate some permanent storage, this is done first time only
 
640
!------------------------------------------------------------------------------
 
641
     IF ( .NOT. ExtractAllocated ) THEN
 
642
       N = Solver % Mesh % MaxElementNodes
 
643
       ALLOCATE( ElementNodes % x( N ), ElementNodes % y( N ), ElementNodes % z( N ),   &
 
644
           ElemVelo( 2, N), ZeroNodes(Solver % Mesh % NumberOfBulkElements,2,2), &
 
645
           STAT=istat )
 
646
       IF ( istat /= 0 ) THEN
 
647
         CALL Fatal( 'LevelSetDistance', 'Memory allocation error 1.' )
 
648
       END IF
 
649
       IF( Convect ) THEN
 
650
         ALLOCATE( Direction( Solver % Mesh % NumberOfBulkElements), STAT=istat )
 
651
       END IF
 
652
       IF ( istat /= 0 ) THEN
 
653
         CALL Fatal( 'RerormalizeSolver', 'Memory allocation error 2.' )
 
654
       END IF
 
655
       ExtractAllocated = .TRUE.
 
656
     END IF
 
657
 
 
658
!------------------------------------------------------------------------------
 
659
!    Extract the zero levelset 
 
660
!------------------------------------------------------------------------------
 
661
 
 
662
     CALL Info( 'LevelSetDistance','--------------------------------------', Level=4 )
 
663
     CALL Info( 'LevelSetDistance','Extracting the zero levelset', Level=4 )
 
664
     CALL Info( 'LevelSetDistance','--------------------------------------', Level=4 )
 
665
 
 
666
     st = CPUTime()
 
667
 
 
668
     CALL ExtractZeroLevel()
 
669
 
 
670
     st = CPUTIme()-st
 
671
     WRITE(Message,'(a,F8.2)') 'Zero level extracted in time (s):',st
 
672
     CALL Info( 'LevelSetDistance',Message, Level=4 )
 
673
 
 
674
     IF( ZeroLevels == 0) THEN
 
675
       CALL Warn('LevelSetDistance','The does not seem to be a zero level-set present, exiting...')
 
676
       RETURN
 
677
     END IF
 
678
 
 
679
     IF(.NOT. Reinitialize) THEN
 
680
       CALL Info('LevelSetDistance','Exiting without reinitialization')
 
681
       RETURN
 
682
     END IF
 
683
 
 
684
     CALL Info( 'LevelSetDistance','--------------------------------------', Level=4 )
 
685
     CALL Info( 'LevelSetDistance','Computing the signed distance function', Level=4 )
 
686
     CALL Info( 'LevelSetDistance','--------------------------------------', Level=4 )
 
687
 
 
688
!------------------------------------------------------------------------------
 
689
!    Allocate some permanent storage for computing the signed distance
 
690
!------------------------------------------------------------------------------
 
691
     IF ( .NOT. DistanceAllocated ) THEN
 
692
       
 
693
       ! The variable for computing the distance
 
694
       DistanceSol => Solver % Variable
 
695
       IF(ASSOCIATED(DistanceSol)) THEN
 
696
         DistancePerm => DistanceSol % Perm
 
697
         Distance => DistanceSol % Values
 
698
       ELSE
 
699
         ALLOCATE(Distance (SIZE(Surface)),STAT=istat)
 
700
         IF ( istat /= 0 ) THEN
 
701
           CALL Fatal( 'LevelSetDistance', 'Memory allocation error 1.' )
 
702
         END IF
 
703
         DistancePerm => SurfPerm
 
704
       END IF
 
705
       DistanceAllocated = .TRUE.
 
706
     END IF
 
707
 
 
708
!------------------------------------------------------------------------------
 
709
!    Compute the signed distance
 
710
!------------------------------------------------------------------------------
 
711
     st = CPUTIme()
 
712
     IF(Convect) THEN
 
713
       DO i=1,Solver % Mesh % NumberOfNodes
 
714
         Distance(DistancePerm(i)) =  ComputeDistanceWithDirection( Solver % Mesh % Nodes % x(i), &
 
715
             Solver % Mesh % Nodes % y(i), 0.0d0, Surface(i) )
 
716
       END DO
 
717
     ELSE
 
718
       DO i=1,Solver % Mesh % NumberOfNodes
 
719
         Distance(DistancePerm(i)) =  ComputeDistance( Solver % Mesh % Nodes % x(i), &
 
720
             Solver % Mesh % Nodes % y(i), 0.0d0 )
 
721
       END DO
 
722
       WHERE( Surface < 0 ) Distance = -Distance
 
723
     END IF
 
724
 
 
725
     n = Solver % Mesh % NumberOfNodes
 
726
     Solver % Variable % Norm = SQRT( SUM(Distance**2)/n )
 
727
 
 
728
!------------------------------------------------------------------------------
 
729
!    Apply the reinitialization to the primary levelset field
 
730
!------------------------------------------------------------------------------
 
731
 
 
732
     IF( ListGetLogical(Solver % Values,'Reinitialize Passive',GotIt) ) THEN
 
733
       CALL Info('LevelSetDistance','Reinitialization not applied to Levelset function')
 
734
     ELSE
 
735
       ! Update also the previous timesteps so that the differentials remain
 
736
       ! unchanged. Otherwise spurious effects are introduced. 
 
737
       IF(ASSOCIATED(SurfSol % PrevValues)) THEN        
 
738
         j = MIN(2, SIZE(SurfSol % PrevValues,2) )
 
739
         IF( ReinitializeInterval > j) THEN                 
 
740
           DO i=1,j
 
741
             SurfSol % PrevValues(:,i) = SurfSol % PrevValues(:,i) + Distance - Surface
 
742
           END DO
 
743
         END IF
 
744
       END IF
 
745
       Surface = Distance
 
746
     END IF
 
747
 
 
748
     st = CPUTIme()-st
 
749
     WRITE(Message,'(a,F8.2)') 'Reinitialization done in time (s):',st
 
750
     CALL Info( 'LevelSetNormalize',Message, Level=4 )
 
751
 
 
752
     IF(Convect) THEN
 
753
       WRITE(Message,'(a,ES12.3)') 'Maximum Levelset Change',dsmax
 
754
       CALL Info( 'LevelSetDistance',Message, Level=4 )     
 
755
       CALL ListAddConstReal(Model % Simulation,'res: LevelSet Max Change',dsmax)
 
756
     END IF
 
757
 
 
758
!------------------------------------------------------------------------------
 
759
 
 
760
CONTAINS
 
761
 
 
762
 
 
763
!------------------------------------------------------------------------------
 
764
  SUBROUTINE ExtractZeroLevel()
 
765
!------------------------------------------------------------------------------
 
766
 
 
767
    INTEGER :: i,j,k,n,div,onetwo,corners, &
 
768
        TriangleIndexes(3), mat_id, body_id, LocalInd(3)
 
769
    TYPE(Variable_t), POINTER :: Var
 
770
    INTEGER, POINTER :: NodeIndexes(:)
 
771
    TYPE(Element_t), POINTER :: Element
 
772
    REAL(KIND=dp) :: x0,y0,x1,y1,t,nx(3),ny(3),nz(3),srf(3),w0(3),w1(3),Value,dt, &
 
773
        r1x, r1y, r2x, r2y, aid, Maxsrf, ds1, ds2
 
774
    LOGICAL :: FileCreated = .FALSE., FileAppend, GotIt, FileSave, Found
 
775
    TYPE(ValueList_t),POINTER :: Material
 
776
    CHARACTER(LEN=MAX_NAME_LEN) :: Filename, Filename2
 
777
    
 
778
    SAVE FileCreated
 
779
!------------------------------------------------------------------------------
 
780
 
 
781
    Filename = ListGetString(Solver % Values,'Filename',FileSave )
 
782
    
 
783
    IF(FileSave) THEN
 
784
      IF(.NOT. FileCreated) THEN        
 
785
        OPEN (10, FILE=TRIM(Filename)//TRIM(".names") )
 
786
        WRITE(10,'(A,A)') 'Variables in file: ',TRIM(Filename)
 
787
        j = 1
 
788
        WRITE(10,'(I3,": ",A)') j,'timestep'
 
789
        
 
790
        Var => Model % Variables
 
791
        DO WHILE( ASSOCIATED( Var ) )          
 
792
          IF ( .NOT. Var % Output .OR. SIZE(Var % Values) == 1 .OR. (Var % DOFs /= 1) ) THEN
 
793
            Var => Var % Next        
 
794
            CYCLE
 
795
          END IF          
 
796
          j = j + 1
 
797
          WRITE(10,'(I3,": ",A)') j,TRIM(Var % Name)
 
798
          Var => Var % Next          
 
799
        END DO
 
800
        CLOSE(10)
 
801
      END IF
 
802
      
 
803
      FileAppend = ListGetLogical(Solver % Values,'File Append',GotIt)
 
804
      IF(FileCreated .AND. FileAppend) THEN 
 
805
        OPEN (10, FILE=Filename,POSITION='APPEND')
 
806
      ELSE 
 
807
        OPEN (10,FILE=Filename)
 
808
      END IF      
 
809
      FileCreated = .TRUE.
 
810
    END IF
 
811
    
 
812
    Surface => SurfSol % Values
 
813
    SurfPerm => SurfSol % Perm
 
814
    
 
815
    ZeroLevels = 0
 
816
    DO i=1,Solver % Mesh % NumberOfBulkElements
 
817
      
 
818
      Element => Solver % Mesh % Elements(i)
 
819
      n = Element % TYPE % NumberOfNodes
 
820
      NodeIndexes => Element % NodeIndexes
 
821
      
 
822
      IF ( ALL( Surface(SurfPerm(NodeIndexes)) < 0) .OR. &
 
823
          ALL( Surface(SurfPerm(NodeIndexes)) > 0) ) CYCLE
 
824
      
 
825
      corners = Element % TYPE % ElementCode / 100 
 
826
      IF(corners < 3 .OR. corners > 4) THEN
 
827
        CALL Warn('ExtractZeroLevel','Implemented only for triangles and quads')
 
828
      END IF
 
829
 
 
830
      IF( Convect ) THEN
 
831
        Model % CurrentElement => Element
 
832
 
 
833
        body_id = Element % BodyId
 
834
        mat_id = ListGetInteger( Model % Bodies(body_id) % Values, 'Material', &
 
835
            minv=1, maxv=Model % NumberOfMaterials )
 
836
        Material => Model % Materials(mat_id) % Values
 
837
 
 
838
        ElemVelo(1,1:n) = ListGetReal( Material,'Levelset Velocity 1',n,NodeIndexes,GotIt)
 
839
        ElemVelo(2,1:n) = ListGetReal( Material,'Levelset Velocity 2',n,NodeIndexes,GotIt)
 
840
      END IF
 
841
 
 
842
 
 
843
      DO div = 1,corners-2
 
844
 
 
845
        SELECT CASE (corners) 
 
846
        CASE (3)
 
847
          LocalInd(1) = 1
 
848
          LocalInd(2) = 2
 
849
          LocalInd(3) = 3
 
850
          
 
851
        CASE(4)               
 
852
          IF(div == 1) THEN
 
853
            LocalInd(1) = 1 
 
854
            LocalInd(2) = 2
 
855
            LocalInd(3) = 4
 
856
          ELSE
 
857
            LocalInd(1) = 2
 
858
            LocalInd(2) = 3
 
859
            LocalInd(3) = 4
 
860
          END IF
 
861
 
 
862
        END SELECT
 
863
        
 
864
        TriangleIndexes = NodeIndexes(LocalInd)
 
865
        srf = Surface(SurfPerm(TriangleIndexes))                       
 
866
        IF ( ALL(srf < 0) .OR. ALL( srf > 0) ) CYCLE
 
867
 
 
868
        nx = Solver % Mesh % Nodes % x(TriangleIndexes)
 
869
        ny = Solver % Mesh % Nodes % y(TriangleIndexes)
 
870
        nz = Solver % Mesh % Nodes % z(TriangleIndexes)
 
871
 
 
872
        CALL TriangleIsoLineWeights( nx,ny,nz,srf,w0,w1,found)
 
873
        IF( .NOT. Found) CYCLE
 
874
 
 
875
        x0 = SUM(w0 * nx)
 
876
        y0 = SUM(w0 * ny)
 
877
        x1 = SUM(w1 * nx)
 
878
        y1 = SUM(w1 * ny)
 
879
 
 
880
        r1x = x1 - x0
 
881
        r1y = y1 - y0 
 
882
 
 
883
        ds1 = SQRT( r1x*r1x + r1y*r1y)
 
884
        IF(ds1 < AEPS) CYCLE
 
885
 
 
886
        ZeroLevels = ZeroLevels + 1
 
887
 
 
888
        IF(Convect) THEN
 
889
 
 
890
          ! Find the value that differs most from the zero levelset
 
891
          j = 1
 
892
          MaxSrf = srf(1)
 
893
          DO k=2,3
 
894
            IF(ABS(srf(k)) > ABS(Maxsrf)) THEN
 
895
              Maxsrf = srf(k)
 
896
              j = k
 
897
            END IF
 
898
          END DO
 
899
          
 
900
          r2x = nx(j) - x0
 
901
          r2y = ny(j) - y0
 
902
          
 
903
          aid = r1x * r2y - r2x * r1y
 
904
          ds2 = SQRT( r2x*r2x + r2y*r2y)
 
905
          
 
906
          Direction( ZeroLevels ) = aid / (ds1 * ds2)
 
907
          IF( Maxsrf < 0.0) THEN
 
908
            Direction( ZeroLevels ) = -Direction( ZeroLevels )
 
909
          END IF
 
910
          
 
911
          r1x = SUM(w0 * ElemVelo(1,LocalInd)) * dt
 
912
          r1y = SUM(w0 * ElemVelo(2,LocalInd)) * dt
 
913
          r2x = SUM(w1 * ElemVelo(1,LocalInd)) * dt
 
914
          r2y = SUM(w1 * ElemVelo(2,LocalInd)) * dt
 
915
 
 
916
          ds1 = SQRT( r1x*r1x + r1y*r1y)
 
917
          ds2 = SQRT( r2x*r2x + r2y*r2y)
 
918
          dsmax = MAX(dsmax, MAX(ds1, ds2) )
 
919
            
 
920
          ZeroNodes(ZeroLevels,1,1) = x0 + r1x
 
921
          ZeroNodes(ZeroLevels,1,2) = y0 + r1y
 
922
          ZeroNodes(ZeroLevels,2,1) = x1 + r2x
 
923
          ZeroNodes(ZeroLevels,2,2) = y1 + r2y
 
924
        ELSE
 
925
          ZeroNodes(ZeroLevels,1,1) = x0
 
926
          ZeroNodes(ZeroLevels,1,2) = y0
 
927
          ZeroNodes(ZeroLevels,2,1) = x1
 
928
          ZeroNodes(ZeroLevels,2,2) = y1
 
929
        END IF
 
930
 
 
931
 
 
932
        IF(FileSave) THEN
 
933
          
 
934
          DO onetwo = 1,2
 
935
            
 
936
            WRITE(10,'(I4)',ADVANCE='NO') Solver % DoneTime
 
937
            
 
938
            Var => Model % Variables
 
939
            DO WHILE( ASSOCIATED( Var ) )
 
940
              
 
941
              IF ( .NOT. Var % Output .OR. SIZE(Var % Values) == 1 .OR. (Var % DOFs /= 1) ) THEN
 
942
                Var => Var % Next        
 
943
                CYCLE
 
944
              END IF
 
945
              
 
946
              Value = 0.0d0
 
947
              DO k=1,3
 
948
                l = TriangleIndexes(k)
 
949
                IF ( ASSOCIATED(Var % Perm) ) l = Var % Perm(l)
 
950
                IF(l > 0) THEN
 
951
                  IF(onetwo == 1) THEN
 
952
                    Value = Value + w0(k) * (Var % Values(l))
 
953
                  ELSE
 
954
                    Value = Value + w1(k) * (Var % Values(l))
 
955
                  END IF
 
956
                END IF
 
957
              END DO
 
958
              
 
959
              WRITE(10,'(ES20.11E3)',ADVANCE='NO') Value
 
960
              Var => Var % Next          
 
961
            END DO
 
962
            WRITE(10,'(A)') ' '
 
963
            
 
964
          END DO
 
965
        END IF        
 
966
      END DO
 
967
 
 
968
    END DO ! of elements
 
969
 
 
970
    
 
971
    IF(FileSave) THEN
 
972
      CLOSE(10)
 
973
    END  IF
 
974
 
 
975
!------------------------------------------------------------------------------
 
976
   END SUBROUTINE ExtractZeroLevel
 
977
!------------------------------------------------------------------------------
 
978
 
 
979
 
 
980
!------------------------------------------------------------------------------
 
981
   SUBROUTINE TriangleIsoLineWeights( NX,NY,NZ,S,w0,w1,Found )
 
982
!------------------------------------------------------------------------------
 
983
! This subroutine extracts the zero line of one triangular element 
 
984
!------------------------------------------------------------------------------
 
985
      IMPLICIT NONE
 
986
      REAL(KIND=dp) :: NX(:),NY(:),NZ(:),S(:),w0(3),w1(3)
 
987
      LOGICAL :: Found
 
988
      REAL(KIND=dp) :: t
 
989
!------------------------------------------------------------------------------
 
990
 
 
991
      Found = .TRUE.
 
992
 
 
993
      w0 = 0.0d0
 
994
      w1 = 0.0d0
 
995
 
 
996
      IF ( ABS(S(1)) < AEPS .AND. ABS(S(2)) < AEPS ) THEN
 
997
        w0(1) = 1.0d0
 
998
        w1(2) = 1.0d0
 
999
      ELSE IF ( ABS(S(1)) < AEPS .AND. ABS(S(3)) < AEPS ) THEN
 
1000
        w0(1) = 1.0d0
 
1001
        w1(3) = 1.0d0
 
1002
      ELSE IF ( ABS(S(2)) < AEPS .AND. ABS(S(3)) < AEPS ) THEN
 
1003
        w0(2) = 1.0d0
 
1004
        w1(3) = 1.0d0
 
1005
      ELSE IF ( ALL(S <= 0) .OR. ALL( S >= 0) ) THEN
 
1006
        Found = .FALSE.
 
1007
      ELSE
 
1008
        IF ( S(1) >= 0 .AND. S(2) >= 0 .OR. &
 
1009
            S(1) <= 0 .AND. S(2) <= 0 ) THEN          
 
1010
          t = -S(1) / ( S(3) - S(1) )
 
1011
          w0(3) = t
 
1012
          w0(1) = 1-t
 
1013
          t = -S(2) / ( S(3) - S(2) )
 
1014
          w1(3) = t
 
1015
          w1(2) = 1-t
 
1016
        ELSE IF ( S(1) >= 0 .AND. S(3) >= 0 .OR. &
 
1017
            S(1) <= 0 .AND. S(3) <= 0 ) THEN
 
1018
          t = -S(1) / ( S(2) - S(1) )
 
1019
          w0(2) = t
 
1020
          w0(1) = 1-t
 
1021
          t = -S(3) / ( S(2) - S(3) )
 
1022
          w1(2) = t
 
1023
          w1(3) = 1-t
 
1024
          
 
1025
        ELSE IF ( S(2) >= 0 .AND. S(3) >= 0 .OR. &
 
1026
            S(2) <= 0 .AND. S(3) <= 0 ) THEN          
 
1027
          t = -S(2) / ( S(1) - S(2) )
 
1028
          w0(1) = t
 
1029
          w0(2) = 1-t
 
1030
          t = -S(3) / ( S(1) - S(3) )
 
1031
          w1(1) = t
 
1032
          w1(3) = 1-t
 
1033
        ELSE 
 
1034
          PRINT *,'TriangleIsoLineWeights: this should not occur'
 
1035
          PRINT *,s(1),s(2),s(3)
 
1036
          STOP
 
1037
        END IF
 
1038
      END IF
 
1039
!------------------------------------------------------------------------------
 
1040
    END SUBROUTINE TriangleIsoLineWeights
 
1041
!------------------------------------------------------------------------------
 
1042
 
 
1043
 
 
1044
!------------------------------------------------------------------------------
 
1045
   FUNCTION ComputeDistance(xp,yp,zp) RESULT(dist)
 
1046
!------------------------------------------------------------------------------
 
1047
! Computes the distance from the given zero levelset given by ZeroNodes. 
 
1048
!------------------------------------------------------------------------------
 
1049
     REAL(KIND=dp) :: xp,yp,zp,dist
 
1050
!------------------------------------------------------------------------------
 
1051
     REAL(KIND=dp) :: x0,y0,x1,y1,a,b,c,d,s
 
1052
     INTEGER :: i,j,k,n
 
1053
!------------------------------------------------------------------------------
 
1054
     dist = HUGE(dist)
 
1055
     DO i=1,ZeroLevels
 
1056
       x0 = ZeroNodes(i,1,1)
 
1057
       y0 = ZeroNodes(i,1,2)
 
1058
       
 
1059
       x1 = ZeroNodes(i,2,1)
 
1060
       y1 = ZeroNodes(i,2,2)
 
1061
       
 
1062
       a = xp - x0
 
1063
       b = x0 - x1
 
1064
       d = y0 - y1
 
1065
       c = yp - y0
 
1066
       s = b**2 + d**2
 
1067
       
 
1068
       x = x0
 
1069
       y = y0
 
1070
       IF ( s > 10*AEPS ) THEN
 
1071
         s = MIN( MAX( -(a*b + c*d) / s, 0.0d0), 1.0d0 )
 
1072
         x = (1-s) * x0 + s * x1
 
1073
         y = (1-s) * y0 + s * y1
 
1074
       END IF
 
1075
       
 
1076
       dist = MIN( dist, SQRT( (xp - x)**2 + (yp - y)**2 ) )
 
1077
     END DO
 
1078
!------------------------------------------------------------------------------
 
1079
   END FUNCTION ComputeDistance
 
1080
!------------------------------------------------------------------------------
 
1081
 
 
1082
 
 
1083
!------------------------------------------------------------------------------
 
1084
   FUNCTION ComputeDistanceWithDirection(xp,yp,zp,prevdist) RESULT(mindist)
 
1085
!------------------------------------------------------------------------------
 
1086
! Computes the distance from the given zero levelset given by ZeroNodes. 
 
1087
!------------------------------------------------------------------------------
 
1088
     REAL(KIND=dp) :: xp,yp,zp,prevdist,mindist
 
1089
!------------------------------------------------------------------------------
 
1090
     REAL(KIND=dp) :: x0,y0,x1,y1,a,b,c,d,s,dist,r1x,r1y,r2x,r2y,angle,angle0
 
1091
     INTEGER :: i,j,k,n
 
1092
!------------------------------------------------------------------------------
 
1093
 
 
1094
     IF(prevdist - dsmax > narrowband) THEN
 
1095
       mindist = prevdist - dsmax
 
1096
       RETURN
 
1097
     ELSE IF(prevdist + dsmax < -narrowband) THEN
 
1098
       mindist = prevdist + dsmax
 
1099
       RETURN
 
1100
     END IF
 
1101
 
 
1102
     mindist = HUGE(mindist)
 
1103
 
 
1104
     DO i=1,ZeroLevels
 
1105
       x0 = ZeroNodes(i,1,1)
 
1106
       y0 = ZeroNodes(i,1,2)
 
1107
       
 
1108
       x1 = ZeroNodes(i,2,1)
 
1109
       y1 = ZeroNodes(i,2,2)
 
1110
       
 
1111
       a = xp - x0
 
1112
       b = x0 - x1
 
1113
       d = y0 - y1
 
1114
       c = yp - y0
 
1115
       s = b**2 + d**2
 
1116
       
 
1117
       x = x0
 
1118
       y = y0
 
1119
       IF ( s > 10*AEPS ) THEN
 
1120
         s = MIN( MAX( -(a*b + c*d) / s, 0.0d0), 1.0d0 )
 
1121
         x = (1-s) * x0 + s * x1
 
1122
         y = (1-s) * y0 + s * y1
 
1123
       END IF
 
1124
       
 
1125
       dist = SQRT( (xp - x)**2 + (yp - y)**2 ) 
 
1126
       
 
1127
 
 
1128
       IF(dist <= (ABS(mindist) + AEPS) ) THEN
 
1129
         
 
1130
         r1x = x1 - x0
 
1131
         r1y = y1 - y0
 
1132
         r2x = xp - x0
 
1133
         r2y = yp - y0
 
1134
         
 
1135
         angle = r1x * r2y - r2x * r1y
 
1136
         
 
1137
         ! Favor parents with clear angles
 
1138
         IF( dist < (ABS(mindist) - AEPS) .OR. (ABS(angle) > ABS(angle0)) ) THEN
 
1139
           IF(Direction(i) * angle < 0.0) THEN
 
1140
             mindist = -dist
 
1141
           ELSE
 
1142
             mindist = dist
 
1143
           END IF
 
1144
           angle0 = angle
 
1145
         END IF
 
1146
           
 
1147
       END IF
 
1148
 
 
1149
      END DO
 
1150
!------------------------------------------------------------------------------
 
1151
    END FUNCTION ComputeDistanceWithDirection
 
1152
!------------------------------------------------------------------------------
 
1153
 
 
1154
 
 
1155
 
 
1156
!------------------------------------------------------------------------------
 
1157
   FUNCTION Solve3rd( Element,Nodes,n,dist,which ) RESULT(val)
 
1158
!------------------------------------------------------------------------------
 
1159
 
 
1160
     REAL(KIND=dp) :: dist(3),val
 
1161
     INTEGER :: n, which
 
1162
     TYPE(Nodes_t) :: Nodes
 
1163
     TYPE(Element_t), POINTER :: Element
 
1164
!------------------------------------------------------------------------------
 
1165
 
 
1166
     REAL(KIND=dp) :: Basis(n), dBasisdx(n,3), ddBasisddx(1,1,1), detJ, &
 
1167
         a,b,c,dx,dy,s,s1,s2,t,u,x0,x1,y0,y1,cosp,ta,tb
 
1168
 
 
1169
     LOGICAL :: stat
 
1170
     INTEGER :: i,k
 
1171
!------------------------------------------------------------------------------
 
1172
 
 
1173
     stat = ElementInfo( Element, Nodes, 1.0d0/3.0d0, 1.0d0/3.0d0, 0.0d0, &
 
1174
         detJ, Basis, dBasisdx, ddBasisddx, .FALSE., .FALSE. )
 
1175
     
 
1176
     dx = 0
 
1177
     dy = 0
 
1178
     DO i=1,n
 
1179
       IF ( i /= Which ) THEN
 
1180
         dx = dx + dBasisdx(i,1) * dist(i)
 
1181
         dy = dy + dBasisdx(i,2) * dist(i)
 
1182
       END IF
 
1183
     END DO
 
1184
     
 
1185
     a = dBasisdx(Which,1)**2 + dBasisdx(Which,2)**2
 
1186
     b = 2*( dBasisdx(Which,1) * dx + dBasisdx(Which,2) * dy )
 
1187
     c = dx**2 + dy**2 - 1
 
1188
     
 
1189
     val = HUGE(val)
 
1190
     s = b**2 - 4*a*c
 
1191
     IF ( s < -1.0d-10 ) THEN
 
1192
       RETURN
 
1193
     ELSE
 
1194
       IF ( s < 0  ) THEN
 
1195
         s1 = -b / (2*a)
 
1196
         s2 = -b / (2*a)
 
1197
       ELSE
 
1198
         IF ( b>0 ) THEN
 
1199
           s1 = 2*c/(-b-SQRT(s))
 
1200
           s2 = (-b-SQRT(s))/(2*a)
 
1201
         ELSE
 
1202
           s1 = (-b+SQRT(s))/(2*a)
 
1203
           s2 = 2*c/(-b+SQRT(s))
 
1204
         END IF
 
1205
       END IF
 
1206
     END IF
 
1207
     
 
1208
     val = MAX( s1, s2 )
 
1209
     
 
1210
     SELECT CASE(Which)
 
1211
     CASE(1)
 
1212
       x0 = Nodes % x(2) - Nodes % x(1)
 
1213
       y0 = Nodes % y(2) - Nodes % y(1)
 
1214
       
 
1215
       x1 = Nodes % x(3) - Nodes % x(1)
 
1216
       y1 = Nodes % y(3) - Nodes % y(1)
 
1217
       
 
1218
       ta = Dist(2); tb = Dist(3)
 
1219
     CASE(2)
 
1220
       x0 = Nodes % x(1) - Nodes % x(2)
 
1221
       y0 = Nodes % y(1) - Nodes % y(2)
 
1222
       
 
1223
       x1 = Nodes % x(3) - Nodes % x(2)
 
1224
       y1 = Nodes % y(3) - Nodes % y(2)
 
1225
       
 
1226
       ta = Dist(1); tb = Dist(3)
 
1227
     CASE(3)
 
1228
       x0 = Nodes % x(1) - Nodes % x(3)
 
1229
       y0 = Nodes % y(1) - Nodes % y(3)
 
1230
       
 
1231
       x1 = Nodes % x(2) - Nodes % x(3)
 
1232
       y1 = Nodes % y(2) - Nodes % y(3)
 
1233
       ta = Dist(1); tb = Dist(2)
 
1234
     END SELECT
 
1235
     
 
1236
     a = SQRT( x0**2 + y0**2 )
 
1237
     b = SQRT( x1**2 + y1**2 )
 
1238
     
 
1239
     dx = dx + dBasisdx(Which,1) * val
 
1240
     dy = dy + dBasisdx(Which,2) * val
 
1241
     s = SQRT(dx**2 + dy**2)
 
1242
     
 
1243
     s1 =  (x0*y1 - y0*x1) / (a*b)
 
1244
     s2 = -(x0*dy - y0*dx) / (s*a)
 
1245
     
 
1246
     IF ( s1 < 0 ) THEN
 
1247
       s1 = -s1;
 
1248
       s2 = -s2;
 
1249
     END IF
 
1250
     
 
1251
     IF ( s2 < 0 .OR. s2 > s1 ) THEN
 
1252
       val = MIN( a + ta, b + tb )
 
1253
     END IF
 
1254
!------------------------------------------------------------------------------
 
1255
   END FUNCTION Solve3rd
 
1256
!------------------------------------------------------------------------------
 
1257
  
 
1258
!------------------------------------------------------------------------------
 
1259
 END SUBROUTINE LevelSetDistance
 
1260
!------------------------------------------------------------------------------
 
1261
 
 
1262
 
 
1263
 
 
1264
 
 
1265
!------------------------------------------------------------------------------
 
1266
   SUBROUTINE LevelSetIntegrate( Model,Solver,Timestep,TransientSimulation )
 
1267
!------------------------------------------------------------------------------
 
1268
!******************************************************************************
 
1269
!
 
1270
!  Compute the volume and area in 3D or area and line intgral in 2D over the 
 
1271
!  levelset function. This is better done within a dedicated solver since 
 
1272
!  it is crucial for the accuracy that the Heaviside and Delta function are
 
1273
!  computed at Gaussian integration points.
 
1274
!
 
1275
!******************************************************************************
 
1276
     USE SolverUtils
 
1277
     USE MaterialModels
 
1278
     USE Integration
 
1279
 
 
1280
     IMPLICIT NONE
 
1281
!------------------------------------------------------------------------------
 
1282
 
 
1283
     TYPE(Model_t), TARGET :: Model
 
1284
     TYPE(Solver_t) :: Solver
 
1285
     REAL(KIND=dp) :: Timestep
 
1286
     LOGICAL :: TransientSimulation
 
1287
 
 
1288
!------------------------------------------------------------------------------
 
1289
!    Local variables
 
1290
!------------------------------------------------------------------------------
 
1291
     TYPE(Nodes_t)   :: ElementNodes
 
1292
     TYPE(Element_t),POINTER :: CurrentElement
 
1293
     TYPE(ValueList_t), POINTER :: Material
 
1294
     TYPE(Variable_t), POINTER :: SurfSol
 
1295
 
 
1296
     INTEGER :: i,j,k,n,t,istat,bf_id
 
1297
     INTEGER, POINTER :: NodeIndexes(:), SurfPerm(:)
 
1298
     LOGICAL :: Visited = .FALSE., GotIt
 
1299
     REAL(KIND=dp), POINTER :: Surface(:), NodalSurf(:)
 
1300
     INTEGER :: body_id, dim
 
1301
     REAL(KIND=dp) :: TotVolume, TotArea, Alpha, Relax, InitVolume, dSurface, dt, &
 
1302
         Moment(3)
 
1303
     CHARACTER(LEN=MAX_NAME_LEN) :: LevelSetVariableName
 
1304
 
 
1305
     SAVE ElementNodes, Visited, NodalSurf, InitVolume
 
1306
 
 
1307
!------------------------------------------------------------------------------
 
1308
!    Get variables needed for solution
 
1309
!------------------------------------------------------------------------------
 
1310
 
 
1311
     ! The variable that should be renormalized
 
1312
     LevelSetVariableName = ListGetString(Solver % Values,'Level Set Variable',GotIt) 
 
1313
     IF(GotIt) THEN
 
1314
       SurfSol => VariableGet( Solver % Mesh % Variables, TRIM(LevelSetVariableName) )
 
1315
     ELSE  
 
1316
       SurfSol => VariableGet( Solver % Mesh % Variables, 'Surface' )
 
1317
     END IF
 
1318
     IF(ASSOCIATED(SurfSol)) THEN
 
1319
       Surface  => SurfSol % Values
 
1320
       SurfPerm => SurfSol % Perm
 
1321
     ELSE
 
1322
       CALL Warn('LevelSetIntegrate','Surface variable does not exist')
 
1323
     END IF
 
1324
 
 
1325
     IF ( ALL( SurfPerm == 0) ) THEN
 
1326
       CALL Warn('LevelSetIntegrate','Nothing to compute')
 
1327
       RETURN
 
1328
     END IF
 
1329
 
 
1330
     dim = CoordinateSystemDimension()
 
1331
 
 
1332
!------------------------------------------------------------------------------
 
1333
!    Allocate some permanent storage, this is done first time only
 
1334
!------------------------------------------------------------------------------
 
1335
     IF ( .NOT. Visited ) THEN
 
1336
       N = Solver % Mesh % MaxElementNodes
 
1337
 
 
1338
       ALLOCATE( NodalSurf( N ), &
 
1339
           ElementNodes % x( N ),   &
 
1340
           ElementNodes % y( N ),   &
 
1341
           ElementNodes % z( N ),   &
 
1342
           STAT=istat )
 
1343
 
 
1344
       IF ( istat /= 0 ) THEN
 
1345
         CALL Fatal( 'LevelSetIntegrate', 'Memory allocation error.' )
 
1346
       END IF
 
1347
     END IF
 
1348
!------------------------------------------------------------------------------
 
1349
!    Do some additional initialization, and go for it
 
1350
!------------------------------------------------------------------------------
 
1351
 
 
1352
     TotVolume = 0.0d0
 
1353
     TotArea = 0.0d0
 
1354
     Moment = 0.0d0
 
1355
     
 
1356
     Alpha = ListGetConstReal(Model % Simulation,'Levelset Bandwidth',GotIt) 
 
1357
     IF(.NOT. GotIt) Alpha = ListGetConstReal(Solver % Values,'Levelset Bandwidth')      
 
1358
       
 
1359
     CALL Info( 'LevelSetIntegrate','-------------------------------------', Level=4 )
 
1360
     CALL Info( 'LevelSetIntegrate', 'Integrating over levelset function', Level=4 )
 
1361
     CALL Info( 'LevelSetIntegrate','-------------------------------------', Level=4 )
 
1362
 
 
1363
     DO t=1,Solver % Mesh % NumberOfBulkElements
 
1364
       
 
1365
       CurrentElement => Solver % Mesh % Elements(t)
 
1366
       n = CurrentElement % TYPE % NumberOfNodes
 
1367
       NodeIndexes => CurrentElement % NodeIndexes
 
1368
       IF( ANY(SurfPerm(NodeIndexes) == 0)) CYCLE
 
1369
       
 
1370
       Model % CurrentElement => CurrentElement
 
1371
       body_id = CurrentElement % Bodyid    
 
1372
       k = ListGetInteger( Model % Bodies( body_id ) % Values, 'Material' )
 
1373
       Material => Model % Materials(k) % Values
 
1374
       
 
1375
!-----------------------------------------------------------------------------
 
1376
!        Get element nodal coordinates
 
1377
!------------------------------------------------------------------------------
 
1378
       ElementNodes % x(1:n) = Solver % Mesh % Nodes % x(NodeIndexes)
 
1379
       ElementNodes % y(1:n) = Solver % Mesh % Nodes % y(NodeIndexes)
 
1380
       ElementNodes % z(1:n) = Solver % Mesh % Nodes % z(NodeIndexes)
 
1381
         
 
1382
       NodalSurf = Surface( SurfPerm(NodeIndexes) )
 
1383
 
 
1384
       CALL HeavisideIntegrate( NodalSurf, CurrentElement, n, ElementNodes, &
 
1385
           Alpha, TotVolume, TotArea, Moment)      
 
1386
     END DO
 
1387
     
 
1388
     Moment = Moment / TotVolume
 
1389
          
 
1390
     IF(dim == 3) THEN
 
1391
       WRITE(Message,'(a,ES12.3)') 'Center 3',Moment(3)
 
1392
       CALL Info( 'LevelSetIntegrate',Message, Level=4 )
 
1393
 
 
1394
       WRITE(Message,'(a,ES12.3)') 'Inside Volume',TotVolume
 
1395
       CALL Info( 'LevelSetIntegrate',Message, Level=4 )
 
1396
     
 
1397
       WRITE(Message,'(a,ES12.3)') 'Interface Area',TotArea
 
1398
       CALL Info( 'LevelSetIntegrate',Message, Level=4 )
 
1399
       
 
1400
       CALL ListAddConstReal(Model % Simulation,'res: LevelSet Center 3',Moment(3))
 
1401
       CALL ListAddConstReal(Model % Simulation,'res: LevelSet Volume',TotVolume)
 
1402
       CALL ListAddConstReal(Model % Simulation,'res: LevelSet Area',TotArea)
 
1403
     ELSE       
 
1404
       WRITE(Message,'(a,ES12.3)') 'Inside Area',TotVolume
 
1405
       CALL Info( 'LevelSetIntegrate',Message, Level=4 )
 
1406
     
 
1407
       WRITE(Message,'(a,ES12.3)') 'Interface length',TotArea
 
1408
       CALL Info( 'LevelSetIntegrate',Message, Level=4 )
 
1409
       
 
1410
       CALL ListAddConstReal(Model % Simulation,'res: LevelSet Area',TotVolume)
 
1411
       CALL ListAddConstReal(Model % Simulation,'res: LevelSet Length',TotArea)       
 
1412
     END IF
 
1413
 
 
1414
     WRITE(Message,'(a,ES12.3)') 'Center 2',Moment(2)
 
1415
     CALL Info( 'LevelSetIntegrate',Message, Level=4 )
 
1416
 
 
1417
     WRITE(Message,'(a,ES12.3)') 'Center 1',Moment(1)
 
1418
     CALL Info( 'LevelSetIntegrate',Message, Level=4 )
 
1419
     
 
1420
     CALL ListAddConstReal(Model % Simulation,'res: LevelSet Center 2',Moment(2))
 
1421
     CALL ListAddConstReal(Model % Simulation,'res: LevelSet Center 1',Moment(1))
 
1422
 
 
1423
     IF ( ListGetLogical(Solver % Values,'Conserve Volume',GotIt) ) THEN
 
1424
       IF(.NOT. Visited) THEN
 
1425
         InitVolume = ListGetConstReal(Solver % Values,'Initial Volume',GotIt)
 
1426
         IF(.NOT. GotIt) InitVolume = TotVolume
 
1427
       END IF
 
1428
 
 
1429
       Relax = ListGetConstReal(Solver % Values,'Conserve Volume Relaxation',GotIt)
 
1430
       IF(.NOT. GotIt) Relax = 1.0d0
 
1431
      
 
1432
       dSurface = Relax * (InitVolume - TotVolume) / TotArea
 
1433
       IF(ABS(dSurface) > AEPS) THEN
 
1434
          Surface = Surface + dSurface
 
1435
       END IF
 
1436
 
 
1437
       WRITE(Message,'(a,ES12.3)') 'Levelset correction',dSurface
 
1438
       CALL Info( 'LevelSetIntegrate',Message, Level=4 )
 
1439
       CALL ListAddConstReal(Model % Simulation,'res: Levelset Correction',dSurface)
 
1440
     END IF
 
1441
 
 
1442
     Visited = .TRUE.
 
1443
 
 
1444
 
 
1445
!------------------------------------------------------------------------------
 
1446
 
 
1447
 CONTAINS
 
1448
 
 
1449
!------------------------------------------------------------------------------
 
1450
   SUBROUTINE HeavisideIntegrate( Surf,Element,n,Nodes,Alpha,Volume,Area,Moment)
 
1451
!------------------------------------------------------------------------------
 
1452
 
 
1453
     REAL(KIND=dp), DIMENSION(:)   :: Surf
 
1454
     INTEGER :: n
 
1455
     TYPE(Nodes_t) :: Nodes
 
1456
     TYPE(Element_t), POINTER :: Element
 
1457
     REAL(KIND=dp) :: Alpha, Volume, Area, Moment(3)
 
1458
 
 
1459
!------------------------------------------------------------------------------
 
1460
!    Local variables
 
1461
!------------------------------------------------------------------------------
 
1462
 
 
1463
     REAL(KIND=dp) :: Basis(n)
 
1464
     REAL(KIND=dp) :: dBasisdx(n,3),detJ
 
1465
     REAL(KIND=dp) :: Val,Grad(3),Velo(3),NormalVelo,GradAbs,Heavi,Delta
 
1466
     INTEGER :: t,N_Integ, CoordinateSystem
 
1467
     REAL(KIND=dp) :: s,u,v,w,x,y,z
 
1468
     TYPE(GaussIntegrationPoints_t), TARGET :: IntegStuff
 
1469
     REAL(KIND=dp), DIMENSION(:), POINTER :: U_Integ,V_Integ,W_Integ,S_Integ
 
1470
     LOGICAL :: stat
 
1471
 
 
1472
!------------------------------------------------------------------------------
 
1473
 
 
1474
     dim = CoordinateSystemDimension()
 
1475
     CoordinateSystem = CurrentCoordinateSystem()
 
1476
 
 
1477
!------------------------------------------------------------------------------
 
1478
!    Integration stuff
 
1479
!------------------------------------------------------------------------------
 
1480
     IntegStuff = GaussPoints( element )
 
1481
 
 
1482
     U_Integ => IntegStuff % u
 
1483
     V_Integ => IntegStuff % v
 
1484
     W_Integ => IntegStuff % w
 
1485
     S_Integ => IntegStuff % s
 
1486
     N_Integ =  IntegStuff % n
 
1487
 
 
1488
!------------------------------------------------------------------------------
 
1489
!    Now we start integrating
 
1490
!------------------------------------------------------------------------------
 
1491
     DO t=1,N_Integ
 
1492
 
 
1493
       u = U_Integ(t)
 
1494
       v = V_Integ(t)
 
1495
       w = W_Integ(t)
 
1496
 
 
1497
!------------------------------------------------------------------------------
 
1498
!      Basis function values & derivatives at the integration point
 
1499
!------------------------------------------------------------------------------
 
1500
       stat = ElementInfo( Element,Nodes,u,v,w,detJ, Basis,dBasisdx)
 
1501
 
 
1502
       s = detJ * S_Integ(t)
 
1503
 
 
1504
       x = SUM( Nodes % x(1:n) * Basis(1:n) )
 
1505
       y = SUM( Nodes % y(1:n) * Basis(1:n) )
 
1506
       IF(dim == 3) z = SUM( Nodes % z(1:n) * Basis(1:n) )
 
1507
 
 
1508
!------------------------------------------------------------------------------
 
1509
!      Coordinatesystem dependent info
 
1510
!------------------------------------------------------------------------------
 
1511
       IF ( CoordinateSystem == AxisSymmetric .OR. CoordinateSystem == CylindricSymmetric ) THEN
 
1512
         s = s * x * 2.0d0 * PI
 
1513
       END IF
 
1514
 
 
1515
!------------------------------------------------------------------------------
 
1516
       Val = SUM( Basis(1:n) * Surf(1:n) )
 
1517
 
 
1518
       IF( Val < -Alpha) THEN
 
1519
         Heavi = 0.0d0
 
1520
       ELSE IF(Val > Alpha) THEN
 
1521
         Heavi = 1.0d0
 
1522
       ELSE
 
1523
         Heavi = (1.0d0 + SIN( (Val/Alpha) * (PI/2) ) ) / 2.0d0
 
1524
         Delta = (1.0d0 + COS( (Val/Alpha) * PI ) ) / (2.0d0 * Alpha)
 
1525
         
 
1526
         DO i=1,dim
 
1527
           Grad(i) = SUM( dBasisdx(1:n,i) * Surf(1:n) )
 
1528
         END DO
 
1529
         GradAbs = SQRT( SUM( Grad(1:dim) * Grad(1:dim) ) )          
 
1530
         
 
1531
         Area = Area + s * Delta * GradAbs
 
1532
       END IF
 
1533
 
 
1534
       Volume = Volume + s * Heavi
 
1535
       Moment(1) = Moment(1) + s * Heavi * x
 
1536
       Moment(2) = Moment(2) + s * Heavi * y
 
1537
       IF(dim == 3) Moment(3) = Moment(3) + s * Heavi * z       
 
1538
 
 
1539
     END DO
 
1540
     
 
1541
!------------------------------------------------------------------------------
 
1542
   END SUBROUTINE HeavisideIntegrate
 
1543
!------------------------------------------------------------------------------
 
1544
 
 
1545
 END SUBROUTINE LevelSetIntegrate
 
1546
!------------------------------------------------------------------------------
 
1547
 
 
1548
 
 
1549
 
 
1550
!------------------------------------------------------------------------------
 
1551
   SUBROUTINE LevelSetCurvature( Model,Solver,Timestep,TransientSimulation )
 
1552
!------------------------------------------------------------------------------
 
1553
!******************************************************************************
 
1554
!
 
1555
!  Computes the curvature from the level set function. Additional diffusion may 
 
1556
!  be added in order to limit the singular curvature peaks.
 
1557
!
 
1558
!******************************************************************************
 
1559
     USE SolverUtils
 
1560
     USE Integration
 
1561
 
 
1562
     IMPLICIT NONE
 
1563
!------------------------------------------------------------------------------
 
1564
 
 
1565
     TYPE(Model_t), TARGET :: Model
 
1566
     TYPE(Solver_t) :: Solver
 
1567
     REAL(KIND=dp) :: Timestep
 
1568
     LOGICAL :: TransientSimulation
 
1569
 
 
1570
!------------------------------------------------------------------------------
 
1571
!    Local variables
 
1572
!------------------------------------------------------------------------------
 
1573
     TYPE(Matrix_t),POINTER  :: StiffMatrix
 
1574
     TYPE(Nodes_t) :: ElementNodes, ParentNodes
 
1575
     TYPE(Element_t),POINTER :: Element, Parent
 
1576
     TYPE(Variable_t), POINTER :: SurfSol
 
1577
 
 
1578
     INTEGER :: i,j,k,n,pn,t,istat,bf_id,CoordinateSystem
 
1579
     REAL(KIND=dp) :: Norm, Coeff, Diff, Alpha, Val, Delta
 
1580
     INTEGER, POINTER :: NodeIndexes(:)
 
1581
     INTEGER, POINTER :: CurvPerm(:), SurfPerm(:)
 
1582
     REAL(KIND=dp), POINTER :: Curvature(:),ForceVector(:), Curv(:),Surface(:) 
 
1583
     REAL(KIND=dp), ALLOCATABLE :: LocalStiffMatrix(:,:),LocalForce(:),Surf(:)
 
1584
     REAL(KIND=dp) :: at,totat,st,totst,CPUTime
 
1585
     CHARACTER(LEN=MAX_NAME_LEN) :: LevelSetVariableName
 
1586
     LOGICAL :: GotIt, Stat, AllocationsDone = .FALSE.
 
1587
 
 
1588
     SAVE LocalStiffMatrix,LocalForce, ElementNodes,ParentNodes,Surf,AllocationsDone
 
1589
 
 
1590
!------------------------------------------------------------------------------
 
1591
!    Get variables needed for solution
 
1592
!------------------------------------------------------------------------------
 
1593
     IF ( .NOT. ASSOCIATED( Solver % Matrix ) ) RETURN
 
1594
 
 
1595
     CoordinateSystem = CurrentCoordinateSystem()
 
1596
 
 
1597
     Curvature => Solver % Variable % Values
 
1598
     CurvPerm => Solver % Variable % Perm
 
1599
     IF ( SIZE( Curvature ) == 0 ) RETURN
 
1600
 
 
1601
     LevelSetVariableName = ListGetString(Solver % Values,'LevelSet Variable',GotIt) 
 
1602
     IF(GotIt) THEN
 
1603
       SurfSol => VariableGet( Solver % Mesh % Variables, TRIM(LevelSetVariableName) )
 
1604
     ELSE  
 
1605
       SurfSol => VariableGet( Solver % Mesh % Variables, 'Surface' )
 
1606
     END IF
 
1607
     Surface => Surfsol % Values
 
1608
     SurfPerm => SurfSol % Perm
 
1609
 
 
1610
     StiffMatrix => Solver % Matrix
 
1611
     ForceVector => StiffMatrix % RHS
 
1612
     Norm = Solver % Variable % Norm
 
1613
 
 
1614
     Diff = ListGetConstReal(Solver % Values,'Curvature Diffusion',GotIt)
 
1615
 
 
1616
!------------------------------------------------------------------------------
 
1617
!    Allocate some permanent storage, this is done first time only
 
1618
!------------------------------------------------------------------------------
 
1619
     IF ( .NOT. AllocationsDone ) THEN
 
1620
       N = Solver % Mesh % MaxElementNodes
 
1621
 
 
1622
       ALLOCATE( ElementNodes % x( N ),   &
 
1623
           ElementNodes % y( N ),   &
 
1624
           ElementNodes % z( N ),   &
 
1625
           ParentNodes % x( N ),   &
 
1626
           ParentNodes % y( N ),   &
 
1627
           ParentNodes % z( N ),   &
 
1628
           LocalForce( N ),       &
 
1629
           LocalStiffMatrix( N, N ), &
 
1630
           Surf( N ), &
 
1631
           STAT=istat )
 
1632
 
 
1633
       IF ( istat /= 0 ) THEN
 
1634
         CALL Fatal( 'CurvatureSolve', 'Memory allocation error.' )
 
1635
       END IF
 
1636
 
 
1637
       AllocationsDone = .TRUE.
 
1638
     END IF
 
1639
!------------------------------------------------------------------------------
 
1640
!    Do some additional initialization, and go for it
 
1641
!------------------------------------------------------------------------------
 
1642
 
 
1643
     totat = 0.0d0
 
1644
     totst = 0.0d0
 
1645
 
 
1646
     at = CPUTime()
 
1647
     
 
1648
     CALL Info( 'LevelSetCurvature','-------------------------------------', Level=4 )
 
1649
     CALL Info( 'LevelSetCurvature','Solving Level set curvature', Level=4 )
 
1650
     CALL Info( 'LevelSetCurvature','-------------------------------------', Level=4 )
 
1651
 
 
1652
     CALL InitializeToZero( StiffMatrix, ForceVector )
 
1653
 
 
1654
!------------------------------------------------------------------------------
 
1655
!      Do the assembly for bulk elements
 
1656
!------------------------------------------------------------------------------
 
1657
 
 
1658
     DO t=1,Solver % NumberOfActiveElements
 
1659
       
 
1660
       Element => Solver % Mesh % Elements(Solver % ActiveElements(t))
 
1661
       n = Element % TYPE % NumberOfNodes
 
1662
       NodeIndexes => Element % NodeIndexes
 
1663
 
 
1664
       ElementNodes % x(1:n) = Solver % Mesh % Nodes % x(NodeIndexes)
 
1665
       ElementNodes % y(1:n) = Solver % Mesh % Nodes % y(NodeIndexes)
 
1666
       ElementNodes % z(1:n) = Solver % Mesh % Nodes % z(NodeIndexes)
 
1667
 
 
1668
       Surf = Surface( SurfPerm(NodeIndexes) )
 
1669
 
 
1670
       CALL LocalMatrix( LocalStiffMatrix, LocalForce, &
 
1671
           Surf, Element, n, ElementNodes )
 
1672
         
 
1673
       CALL UpdateGlobalEquations( StiffMatrix, LocalStiffMatrix, &
 
1674
           ForceVector, LocalForce, n, 1, CurvPerm(NodeIndexes) )
 
1675
 
 
1676
     END DO     
 
1677
 
 
1678
!------------------------------------------------------------------------------
 
1679
!      Do the assembly for boundary elements
 
1680
!------------------------------------------------------------------------------
 
1681
 
 
1682
     DO t=Solver % Mesh % NumberOfBulkElements + 1, &
 
1683
         Solver % Mesh % NumberOfBulkElements + &
 
1684
         Solver % Mesh % NumberOfBoundaryElements
 
1685
      
 
1686
       Element => Solver % Mesh % Elements(t)
 
1687
!------------------------------------------------------------------------------
 
1688
       DO i=1,Model % NumberOfBCs
 
1689
         IF ( Element % BoundaryInfo % Constraint == Model % BCs(i) % Tag ) THEN
 
1690
 
 
1691
           n = Element % TYPE % NumberOfNodes
 
1692
           NodeIndexes => Element % NodeIndexes
 
1693
           IF ( ANY( CurvPerm(NodeIndexes) <= 0 ) ) CYCLE         
 
1694
 
 
1695
           IF ( Element % TYPE % ElementCode == 101 ) CYCLE
 
1696
           
 
1697
           IF ( .NOT. ListGetLogical(Model % BCs(i) % Values, &
 
1698
               'Levelset Curvature BC',gotIt) ) CYCLE
 
1699
 
 
1700
           ElementNodes % x(1:n) = Solver % Mesh % Nodes % x(NodeIndexes)
 
1701
           ElementNodes % y(1:n) = Solver % Mesh % Nodes % y(NodeIndexes)
 
1702
           ElementNodes % z(1:n) = Solver % Mesh % Nodes % z(NodeIndexes)
 
1703
           
 
1704
           Parent => Element % BoundaryInfo % Left
 
1705
          
 
1706
           stat = ASSOCIATED( Parent )
 
1707
           IF ( stat ) stat = stat .AND. ALL(CurvPerm(Parent % NodeIndexes) > 0)
 
1708
          
 
1709
           IF ( .NOT. stat ) THEN
 
1710
             Parent => Element % BoundaryInfo % Right
 
1711
             
 
1712
            stat = ASSOCIATED( Parent )
 
1713
            IF ( stat ) stat = ALL(CurvPerm(Parent % NodeIndexes(1:pn)) > 0)
 
1714
            
 
1715
            IF ( .NOT. stat )  THEN
 
1716
              CALL Warn( 'LevelSetCurvature', &
 
1717
                  'No curvature solution available for specified boundary' )
 
1718
              CYCLE
 
1719
            END IF
 
1720
          END IF
 
1721
          
 
1722
          pn = Parent % TYPE % NumberOfNodes           
 
1723
          ParentNodes % x(1:pn) = Solver % Mesh % Nodes % x(Parent % NodeIndexes)
 
1724
          ParentNodes % y(1:pn) = Solver % Mesh % Nodes % y(Parent % NodeIndexes)
 
1725
          ParentNodes % z(1:pn) = Solver % Mesh % Nodes % z(Parent % NodeIndexes)
 
1726
 
 
1727
          Surf(1:pn) = Surface(SurfPerm(Parent % NodeIndexes))
 
1728
          
 
1729
!------------------------------------------------------------------------------
 
1730
!             Get element matrix and rhs due to boundary conditions ...
 
1731
!------------------------------------------------------------------------------
 
1732
          
 
1733
          CALL LocalBoundary( LocalStiffMatrix, LocalForce,  &
 
1734
              Surf, Element, Parent, n, pn, ElementNodes, ParentNodes )
 
1735
!------------------------------------------------------------------------------
 
1736
!             Update global matrices from local matrices
 
1737
!------------------------------------------------------------------------------
 
1738
          CALL UpdateGlobalEquations( StiffMatrix, LocalStiffMatrix, &
 
1739
              ForceVector, LocalForce, n, 1, CurvPerm(NodeIndexes) )
 
1740
!------------------------------------------------------------------------------
 
1741
           END IF   
 
1742
        END DO   
 
1743
      END DO   
 
1744
!------------------------------------------------------------------------------
 
1745
 
 
1746
     CALL FinishAssembly( Solver,ForceVector )
 
1747
     
 
1748
     at = CPUTime() - at
 
1749
     WRITE(Message,'(a,F8.2)') 'Assembly done in time (s):',at
 
1750
     CALL Info( 'LevelSetCurvature',Message, Level=4 )
 
1751
 
 
1752
!------------------------------------------------------------------------------
 
1753
!    Solve the system and we are done.
 
1754
!------------------------------------------------------------------------------
 
1755
     st = CPUTime()
 
1756
     CALL SolveSystem( StiffMatrix, ParMatrix, ForceVector, &
 
1757
         Curvature, Norm, 1, Solver )
 
1758
     
 
1759
     st = CPUTIme()-st
 
1760
     WRITE(Message,'(a,F8.2)') 'Solution done in time (s):',st
 
1761
     CALL Info( 'LevelSetCurvature',Message, Level=4 )
 
1762
 
 
1763
!------------------------------------------------------------------------------
 
1764
 
 
1765
     Coeff = ListGetConstReal(Solver % Values,'Curvature Coefficient',GotIt) 
 
1766
     IF(GotIt) THEN
 
1767
       Curvature = Coeff * Curvature
 
1768
     END IF
 
1769
 
 
1770
     Alpha = ListGetConstReal(Model % Simulation,'Levelset Bandwidth',GotIt) 
 
1771
     IF(.NOT. GotIt) Alpha = ListGetConstReal(Solver % Values,'Levelset Bandwidth',GotIt)      
 
1772
     IF(GotIt) THEN
 
1773
       DO i=1,SIZE(CurvPerm) 
 
1774
         j = CurvPerm(i)
 
1775
         k = SurfPerm(i)
 
1776
         IF(j == 0 .OR. k == 0) CYCLE
 
1777
         Val = Surface(k)
 
1778
 
 
1779
         IF( Val < -Alpha) THEN
 
1780
           Delta = 0.0d0
 
1781
         ELSE IF(Val > Alpha) THEN
 
1782
           Delta = 0.0d0
 
1783
         ELSE
 
1784
           Delta = (1.0d0 + COS( (Val/Alpha) * PI ) ) / (2.0d0 * Alpha)
 
1785
         END IF
 
1786
         
 
1787
         Curvature(j) = Delta * Curvature(j)
 
1788
       END DO
 
1789
     END IF
 
1790
 
 
1791
 
 
1792
CONTAINS
 
1793
 
 
1794
   SUBROUTINE LocalMatrix( StiffMatrix,ForceVector, Surf, &
 
1795
       Element, n, Nodes )
 
1796
!------------------------------------------------------------------------------
 
1797
 
 
1798
     REAL(KIND=dp), DIMENSION(:)   :: ForceVector,Surf
 
1799
     REAL(KIND=dp), DIMENSION(:,:) :: StiffMatrix
 
1800
     INTEGER :: n
 
1801
     TYPE(Nodes_t) :: Nodes
 
1802
     TYPE(Element_t), POINTER :: Element
 
1803
 
 
1804
!------------------------------------------------------------------------------
 
1805
!    Local variables
 
1806
!------------------------------------------------------------------------------
 
1807
!
 
1808
     REAL(KIND=dp) :: Basis(n)
 
1809
     REAL(KIND=dp) :: dBasisdx(n,3),detJ
 
1810
     REAL(KIND=dp) :: Grad(3),GradAbs,A,B,xpos
 
1811
     INTEGER :: i,j,k,c,p,q,t,dim,N_Integ
 
1812
     REAL(KIND=dp) :: s,u,v,w
 
1813
     TYPE(GaussIntegrationPoints_t), TARGET :: IntegStuff
 
1814
     REAL(KIND=dp), DIMENSION(:), POINTER :: U_Integ,V_Integ,W_Integ,S_Integ
 
1815
 
 
1816
     LOGICAL :: stat
 
1817
 
 
1818
!------------------------------------------------------------------------------
 
1819
 
 
1820
     dim = CoordinateSystemDimension()
 
1821
     Grad = 0.0d0
 
1822
     ForceVector = 0.0D0
 
1823
     StiffMatrix = 0.0D0
 
1824
 
 
1825
!------------------------------------------------------------------------------
 
1826
!    Integration stuff
 
1827
!------------------------------------------------------------------------------
 
1828
     IntegStuff = GaussPoints( element )
 
1829
     U_Integ => IntegStuff % u
 
1830
     V_Integ => IntegStuff % v
 
1831
     W_Integ => IntegStuff % w
 
1832
     S_Integ => IntegStuff % s
 
1833
     N_Integ =  IntegStuff % n
 
1834
 
 
1835
!------------------------------------------------------------------------------
 
1836
!    Now we start integrating
 
1837
!------------------------------------------------------------------------------
 
1838
     DO t=1,N_Integ
 
1839
 
 
1840
       u = U_Integ(t)
 
1841
       v = V_Integ(t)
 
1842
       w = W_Integ(t)
 
1843
       
 
1844
!------------------------------------------------------------------------------
 
1845
!      Basis function values & derivatives at the integration point
 
1846
!------------------------------------------------------------------------------
 
1847
       stat = ElementInfo( Element,Nodes,u,v,w,detJ,Basis,dBasisdx)
 
1848
 
 
1849
       s = detJ * S_Integ(t)
 
1850
 
 
1851
!------------------------------------------------------------------------------
 
1852
!      Coordinatesystem dependent info
 
1853
!------------------------------------------------------------------------------
 
1854
       IF ( Coordinates == AxisSymmetric .OR. Coordinates == CylindricSymmetric ) THEN
 
1855
         xpos = SUM( Nodes % x(1:n) * Basis(1:n) )
 
1856
         s = xpos * s
 
1857
       END IF
 
1858
 
 
1859
!------------------------------------------------------------------------------
 
1860
 
 
1861
       DO i=1,dim
 
1862
         Grad(i) = SUM( dBasisdx(1:n,i) * Surf(1:n) )
 
1863
       END DO
 
1864
       GradAbs = SQRT( SUM( Grad * Grad ) )
 
1865
       IF ( GradAbs > 10*AEPS ) THEN
 
1866
         Grad = Grad / GradAbs
 
1867
       END IF
 
1868
 
 
1869
!------------------------------------------------------------------------------
 
1870
!      Loop over basis functions of both unknowns and weights
 
1871
!------------------------------------------------------------------------------
 
1872
       DO p=1,n
 
1873
 
 
1874
         DO q=1,n
 
1875
           A = Basis(q) * Basis(p)
 
1876
           A = A + Diff * SUM( dBasisdx(q,1:dim) * dBasisdx(p,1:dim) )
 
1877
 
 
1878
           StiffMatrix(p,q) = StiffMatrix(p,q) + s * A
 
1879
         END DO
 
1880
 
 
1881
         B = - SUM( dBasisdx(p, 1:dim) * Grad(1:dim))  
 
1882
         ForceVector(p) = ForceVector(p) + s * B        
 
1883
 
 
1884
       END DO
 
1885
     END DO
 
1886
 
 
1887
!------------------------------------------------------------------------------
 
1888
   END SUBROUTINE LocalMatrix
 
1889
 
 
1890
 
 
1891
!------------------------------------------------------------------------------
 
1892
   SUBROUTINE LocalBoundary( BoundaryMatrix, BoundaryVector, &
 
1893
        Surf, Element, Parent, n, pn, ElementNodes, ParentNodes )
 
1894
!------------------------------------------------------------------------------
 
1895
 
 
1896
     REAL(KIND=dp) :: BoundaryMatrix(:,:), BoundaryVector(:), Surf(:)
 
1897
     TYPE(Nodes_t)   :: ElementNodes, ParentNodes
 
1898
     TYPE(Element_t), POINTER :: Element, Parent
 
1899
     INTEGER :: n, pn
 
1900
 
 
1901
     REAL(KIND=dp) :: Basis(n), dBasisdx(n,3), Normal(3), Grad(3), GradAbs, detJ, A
 
1902
     REAL(KIND=dp) :: ParentBasis(pn), ParentdBasisdx(pn,3)
 
1903
     REAL(KIND=dp) :: u,v,w,s,x(n),y(n),z(n),xpos, Force
 
1904
     REAL(KIND=dp), POINTER :: U_Integ(:),V_Integ(:),W_Integ(:),S_Integ(:)
 
1905
 
 
1906
     INTEGER :: i,j,k,t,p,q,N_Integ,dim
 
1907
 
 
1908
     TYPE(GaussIntegrationPoints_t), TARGET :: IntegStuff
 
1909
 
 
1910
     LOGICAL :: stat
 
1911
!------------------------------------------------------------------------------
 
1912
 
 
1913
     BoundaryVector = 0.0d0
 
1914
     BoundaryMatrix = 0.0d0
 
1915
     dim = CoordinateSystemDimension()
 
1916
     Grad = 0.0d0
 
1917
 
 
1918
!------------------------------------------------------------------------------
 
1919
!    Integration stuff
 
1920
!------------------------------------------------------------------------------
 
1921
     IntegStuff = GaussPoints( Element )
 
1922
     U_Integ => IntegStuff % u
 
1923
     V_Integ => IntegStuff % v
 
1924
     W_Integ => IntegStuff % w
 
1925
     S_Integ => IntegStuff % s
 
1926
     N_Integ =  IntegStuff % n
 
1927
 
 
1928
!------------------------------------------------------------------------------
 
1929
!   Now we start integrating
 
1930
!------------------------------------------------------------------------------
 
1931
     DO t=1,N_Integ
 
1932
       u = U_Integ(t)
 
1933
       v = V_Integ(t)
 
1934
       w = W_Integ(t)
 
1935
!------------------------------------------------------------------------------
 
1936
!     Basis function values & derivates at the integration point
 
1937
!------------------------------------------------------------------------------
 
1938
       stat = ElementInfo( Element,ElementNodes,u,v,w,detJ,Basis,dBasisdx )
 
1939
       s = S_Integ(t) * detJ
 
1940
 
 
1941
!------------------------------------------------------------------------------
 
1942
!      Coordinatesystem dependent info
 
1943
!------------------------------------------------------------------------------
 
1944
      IF ( CurrentCoordinateSystem() /= Cartesian ) THEN
 
1945
        xpos = SUM( ElementNodes % x(1:n)*Basis(1:n) )
 
1946
        s = xpos * s
 
1947
      END IF
 
1948
 
 
1949
      Normal = Normalvector( Element, ElementNodes, u, v, .TRUE. )
 
1950
      
 
1951
      !------------------------------------------------------------------------------
 
1952
      ! Need parent element basis etc., for computing normal derivatives on boundary.
 
1953
      !------------------------------------------------------------------------------
 
1954
      
 
1955
      DO i = 1,n
 
1956
        DO j = 1,pn
 
1957
          IF ( Element % NodeIndexes(i) == Parent % NodeIndexes(j) ) THEN
 
1958
            x(i) = Parent % TYPE % NodeU(j)
 
1959
            y(i) = Parent % TYPE % NodeV(j)
 
1960
            z(i) = Parent % TYPE % NodeW(j)
 
1961
            EXIT
 
1962
          END IF
 
1963
        END DO
 
1964
      END DO
 
1965
      
 
1966
      u = SUM( Basis(1:n) * x(1:n) )
 
1967
      v = SUM( Basis(1:n) * y(1:n) )
 
1968
      w = SUM( Basis(1:n) * z(1:n) )
 
1969
      
 
1970
      stat = ElementInfo( Parent, ParentNodes, u, v, w, detJ, ParentBasis, ParentdBasisdx )
 
1971
          
 
1972
      DO j = 1,DIM
 
1973
        Grad(j) = SUM( ParentdBasisdx(1:pn,j) * Surf(1:pn) )
 
1974
      END DO
 
1975
      GradAbs = SQRT( SUM(Grad(1:dim) * Grad(1:dim)) )
 
1976
 
 
1977
      IF ( GradAbs > 10*AEPS ) THEN         
 
1978
        Grad = Grad / GradAbs
 
1979
      END IF
 
1980
 
 
1981
!------------------------------------------------------------------------------
 
1982
      
 
1983
      DO p=1,n
 
1984
        DO q=1,n
 
1985
          A = Diff * SUM( Normal(1:dim) * dBasisdx(q,1:dim)) * Basis(p)
 
1986
          BoundaryMatrix(p,q) = BoundaryMatrix(p,q) + s * A
 
1987
        END DO
 
1988
      END DO
 
1989
 
 
1990
      Force = SUM( Normal(1:dim) * Grad(1:dim) )
 
1991
      DO q=1,N
 
1992
        BoundaryVector(q) = BoundaryVector(q) + s * Basis(q) * Force
 
1993
      END DO
 
1994
 
 
1995
    END DO
 
1996
 
 
1997
  END SUBROUTINE LocalBoundary
 
1998
!------------------------------------------------------------------------------
 
1999
 
 
2000
!------------------------------------------------------------------------------
 
2001
END SUBROUTINE LevelSetCurvature
 
2002
!------------------------------------------------------------------------------
 
2003
 
 
2004
 
 
2005
 
 
2006
 
 
2007
FUNCTION LevelSetTimestep( Model ) RESULT( dt )
 
2008
  USE Types
 
2009
  USE Lists
 
2010
  USE Integration
 
2011
  USE ElementDescription
 
2012
 
 
2013
  IMPLICIT NONE
 
2014
 
 
2015
  TYPE(Model_t) :: Model
 
2016
  REAL(KIND=dp) :: dt
 
2017
 
 
2018
  TYPE(Solver_t), POINTER :: Solver
 
2019
  TYPE(Variable_t), POINTER ::  TimeVariable, SurfSol
 
2020
  REAL(KIND=dp), POINTER :: Surface(:), Surf(:)
 
2021
  INTEGER, POINTER :: SurfPerm(:), NodeIndexes(:)
 
2022
  TYPE(Nodes_t)   :: ElementNodes
 
2023
  TYPE(Element_t),POINTER :: CurrentElement
 
2024
  TYPE(ValueList_t), POINTER :: Material
 
2025
  REAL(KIND=dp), POINTER :: TimestepSizes(:,:)
 
2026
 
 
2027
  INTEGER :: i,j,k,n,t,elem,N_Integ,body_id,mat_id, dim, TimeIntervals
 
2028
  REAL(KIND=dp), ALLOCATABLE :: Basis(:), dBasisdx(:,:), NodalVelo(:,:)
 
2029
  REAL(KIND=dp) :: Val,Grad(3),Velo(3),NormalVelo,AbsVelo,GradAbs,detJ,&
 
2030
      MaxNormVelo,MaxAbsVelo
 
2031
  REAL(KIND=dp) :: s,u,v,w,dt0,prevdt,cumtime,dsmax
 
2032
  TYPE(GaussIntegrationPoints_t), TARGET :: IntegStuff
 
2033
  REAL(KIND=dp), DIMENSION(:), POINTER :: U_Integ,V_Integ,W_Integ,S_Integ
 
2034
  LOGICAL :: stat, GotIt, AllocationsDone = .FALSE.
 
2035
  CHARACTER(LEN=MAX_NAME_LEN) :: LevelSetVariableName
 
2036
 
 
2037
  SAVE AllocationsDone, Basis, dBasisdx, NodalVelo, Surf, ElementNodes, prevdt
 
2038
 
 
2039
 
 
2040
  TimestepSizes => ListGetConstRealArray( CurrentModel % Simulation,&
 
2041
      'Timestep Sizes', GotIt )
 
2042
  TimeIntervals = SIZE(TimestepSizes)
 
2043
 
 
2044
  IF(TimeIntervals > 1) THEN
 
2045
    CALL Warn('LevelSetTimestep','Implemented only for one Time Interval')
 
2046
  END IF
 
2047
 
 
2048
  dt0 = TimestepSizes(1,1)
 
2049
 
 
2050
  TimeVariable => VariableGet(Model % Variables, 'Time')
 
2051
  cumtime = TimeVariable % Values(1)
 
2052
 
 
2053
  dim = CoordinateSystemDimension()
 
2054
 
 
2055
  ! The variable that should be reinitialized
 
2056
  LevelSetVariableName = 'Surface'
 
2057
  SurfSol => VariableGet( Model % Variables, TRIM(LevelSetVariableName) )
 
2058
  IF(ASSOCIATED(SurfSol)) THEN
 
2059
    Surface  => SurfSol % Values
 
2060
    SurfPerm => SurfSol % Perm
 
2061
  ELSE
 
2062
    CALL Warn('LevelSetTimeStep','SurfSol does not exist: '//TRIM(LevelSetVariableName))
 
2063
    RETURN
 
2064
  END IF
 
2065
 
 
2066
  Solver => SurfSol % Solver
 
2067
  DsMax = ListGetConstReal(Model % Simulation,'LevelSet Courant Number',GotIt)
 
2068
  IF(.NOT. GotIt) DsMax = 1.0d0
 
2069
 
 
2070
  IF(.NOT. AllocationsDone) THEN
 
2071
    N = Solver % Mesh % MaxElementNodes
 
2072
    ALLOCATE( Basis(n), dBasisdx(n,3), NodalVelo(3,n), Surf(n), &
 
2073
        ElementNodes % x(n), ElementNodes % y(n), ElementNodes % z(n) )
 
2074
    prevdt = dt0
 
2075
    AllocationsDone = .TRUE. 
 
2076
  END IF
 
2077
 
 
2078
  MaxNormVelo = 0.0d0
 
2079
  MaxAbsVelo = 0.0d0
 
2080
 
 
2081
  DO elem=1,Solver % Mesh % NumberOfBulkElements
 
2082
    
 
2083
    CurrentElement => Solver % Mesh % Elements(elem)
 
2084
    n = CurrentElement % TYPE % NumberOfNodes
 
2085
    NodeIndexes => CurrentElement % NodeIndexes
 
2086
    
 
2087
    IF(ANY(SurfPerm(NodeIndexes) == 0)) CYCLE
 
2088
    
 
2089
    Surf(1:n) = Surface( SurfPerm(NodeIndexes) )
 
2090
    IF(ALL(Surf(1:n) < 0.0d0) .OR. ALL(Surf(1:n) > 0.0d0) ) CYCLE
 
2091
 
 
2092
    ElementNodes % x(1:n) = Solver % Mesh % Nodes % x(NodeIndexes)
 
2093
    ElementNodes % y(1:n) = Solver % Mesh % Nodes % y(NodeIndexes)
 
2094
    ElementNodes % z(1:n) = Solver % Mesh % Nodes % z(NodeIndexes)
 
2095
    
 
2096
    Model % CurrentElement => CurrentElement
 
2097
    body_id = CurrentElement % Bodyid    
 
2098
    mat_id = ListGetInteger( Model % Bodies( body_id ) % Values, 'Material' )
 
2099
    Material => Model % Materials(mat_id) % Values
 
2100
 
 
2101
!------------------------------------------------------------------------------
 
2102
!         Computed or given velocity field
 
2103
!------------------------------------------------------------------------------
 
2104
       
 
2105
    NodalVelo(1,1:n) = ListGetReal( Material,'Levelset Velocity 1',n,NodeIndexes,GotIt)
 
2106
    NodalVelo(2,1:n) = ListGetReal( Material,'Levelset Velocity 2',n,NodeIndexes,GotIt)
 
2107
    IF(dim == 3) NodalVelo(3,1:n) = ListGetReal( Material,'Levelset Velocity 3',n,NodeIndexes,GotIt)
 
2108
 
 
2109
!------------------------------------------------------------------------------
 
2110
!    Integration stuff
 
2111
!------------------------------------------------------------------------------
 
2112
    IntegStuff = GaussPoints( CurrentElement )
 
2113
    
 
2114
    U_Integ => IntegStuff % u
 
2115
    V_Integ => IntegStuff % v
 
2116
    W_Integ => IntegStuff % w
 
2117
    S_Integ => IntegStuff % s
 
2118
    N_Integ =  IntegStuff % n
 
2119
 
 
2120
!------------------------------------------------------------------------------
 
2121
!    Maximum at any integration point
 
2122
!------------------------------------------------------------------------------
 
2123
    
 
2124
    DO t=1,N_Integ
 
2125
 
 
2126
      u = U_Integ(t)
 
2127
      v = V_Integ(t)
 
2128
      w = W_Integ(t)
 
2129
         
 
2130
!------------------------------------------------------------------------------
 
2131
!      Basis function values & derivatives at the integration point
 
2132
!------------------------------------------------------------------------------
 
2133
 
 
2134
      stat = ElementInfo( CurrentElement,ElementNodes,u,v,w,detJ, Basis,dBasisdx)
 
2135
      
 
2136
      DO i=1,dim
 
2137
        Grad(i) = SUM( dBasisdx(1:n,i) * Surf(1:n) )
 
2138
        Velo(i) = SUM( Basis(1:n) * NodalVelo(i,1:n) )
 
2139
      END DO
 
2140
      
 
2141
      GradAbs = SQRT( SUM( Grad(1:dim) * Grad(1:dim) ) )         
 
2142
      NormalVelo = SUM( Grad(1:dim) * Velo(1:dim) ) / GradAbs
 
2143
      NormalVelo = NormalVelo / SQRT(detJ)
 
2144
      MaxNormVelo = MAX(MaxNormVelo, ABS(NormalVelo))
 
2145
      
 
2146
      AbsVelo = SQRT(SUM (Velo(1:dim) * Velo(1:dim)) )
 
2147
      AbsVelo = AbsVelo / SQRT(detJ)
 
2148
      MaxAbsVelo = MAX(MaxAbsVelo, ABS(AbsVelo))
 
2149
    END DO
 
2150
  END DO
 
2151
 
 
2152
  dt = dt0
 
2153
  IF( ListGetLogical(Model % Simulation,'LevelSet Timestep Directional',GotIt) ) THEN
 
2154
    IF( MaxNormVelo * dt0 > dsMax) dt = dsMax / MaxNormVelo
 
2155
  ELSE
 
2156
    IF( MaxAbsVelo * dt0 > dsMax) dt = dsMax / MaxAbsVelo
 
2157
  END IF
 
2158
  
 
2159
  IF( dt < dt0) THEN
 
2160
    IF(dt > prevdt) dt = 0.5d0 * (dt + prevdt)
 
2161
  END IF
 
2162
  prevdt = dt
 
2163
  
 
2164
  WRITE(Message,'(a,ES12.3)') 'Levelset timestep',dt
 
2165
  CALL Info( 'LevelSetTimestep',Message, Level=4 )
 
2166
  CALL ListAddConstReal(Model % Simulation,'res: Levelset timestep',dt)
 
2167
 
 
2168
END FUNCTION LevelSetTimestep
 
2169