CoCalc -- AdjointThickness_GradientSolver.F90

GitHub Repository: ElmerCSC/elmerfem
Path: blob/devel/elmerice/Solvers/AdjointThickness/AdjointThickness_GradientSolver.F90
³²⁰⁶ views
1
!/*****************************************************************************/
2
! *
3
! *  Elmer/Ice, a glaciological add-on to Elmer
4
! *  http://elmerice.elmerfem.org
5
! *
6
! * 
7
! *  This program is free software; you can redistribute it and/or
8
! *  modify it under the terms of the GNU General Public License
9
! *  as published by the Free Software Foundation; either version 2
10
! *  of the License, or (at your option) any later version.
11
! * 
12
! *  This program is distributed in the hope that it will be useful,
13
! *  but WITHOUT ANY WARRANTY; without even the implied warranty of
14
! *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
15
! *  GNU General Public License for more details.
16
! *
17
! *  You should have received a copy of the GNU General Public License
18
! *  along with this program (in file fem/GPL-2); if not, write to the 
19
! *  Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, 
20
! *  Boston, MA 02110-1301, USA.
21
! *
22
! *****************************************************************************/
23
! ******************************************************************************
24
! *
25
! *  Authors: F. Gillet-Chaulet
26
! *  Web:     http://elmerice.elmerfem.org
27
! *
28
! *  Original Date: Dec. 2020
29
! * 
30
! *****************************************************************************
31
!
32
SUBROUTINE AdjointThickness_GradientSolver_init0(Model,Solver,dt,TransientSimulation )
33
!------------------------------------------------------------------------------
34
  USE DefUtils
35
  IMPLICIT NONE
36
!------------------------------------------------------------------------------
37
  TYPE(Solver_t), TARGET :: Solver
38
  TYPE(Model_t) :: Model
39
  REAL(KIND=dp) :: dt
40
  LOGICAL :: TransientSimulation
41
!------------------------------------------------------------------------------
42
! Local variables
43
!------------------------------------------------------------------------------
44
  CHARACTER(LEN=MAX_NAME_LEN) :: Name
45

46
  Name = ListGetString( Solver % Values, 'Equation',UnFoundFatal=.TRUE.)
47
  CALL ListAddNewString( Solver % Values,'Variable',&
48
          '-nooutput '//TRIM(Name)//'_var')
49
  CALL ListAddLogical(Solver % Values, 'Optimize Bandwidth',.FALSE.)
50

51
END SUBROUTINE AdjointThickness_GradientSolver_init0
52
!*****************************************************************************
53
!-----------------------------------------------------------------------------
54
SUBROUTINE AdjointThickness_GradientSolver( Model,Solver,dt,TransientSimulation )
55
  USE DefUtils
56
  IMPLICIT NONE
57

58
  !------------------------------------------------------------------------------
59
  !    external variables
60
  !------------------------------------------------------------------------------
61
  TYPE(Model_t) :: Model
62
  TYPE(Solver_t):: Solver
63
  REAL(KIND=dp) :: dt
64
  LOGICAL :: TransientSimulation
65
  !------------------------------------------------------------------------------
66
  !    Local variables
67
  !------------------------------------------------------------------------------
68
  TYPE(ValueList_t), POINTER :: BodyForce, SolverParams
69

70
  INTEGER :: NMAX,NSDOFs
71
  INTEGER :: istat
72
  INTEGER :: t,i,j,n
73
  LOGICAL :: ConvectionVar, Found
74

75
  REAL(KIND=dp) :: Norm
76

77
  CHARACTER(LEN=MAX_NAME_LEN)  :: FlowSolName
78
  TYPE(Variable_t), POINTER :: FlowSol
79
  REAL(KIND=dp), POINTER :: FlowSolution(:)
80
  INTEGER, POINTER :: FlowPerm(:)
81

82
  REAL(KIND=dp), ALLOCATABLE,SAVE :: STIFF(:,:),FORCE(:),LOAD(:),Velo(:,:)
83

84
  CHARACTER(LEN=MAX_NAME_LEN)  :: SolverName= 'AdjointThickness_ThicknessSolver'
85

86
  TYPE(Nodes_t),SAVE   :: ElementNodes
87
  TYPE(Element_t),POINTER :: CurrentElement
88
  INTEGER, POINTER :: NodeIndexes(:)
89

90
  LOGICAL, SAVE :: AllocationsDone = .FALSE.
91

92
  CHARACTER(LEN=MAX_NAME_LEN) :: ThickSolName,ThickbSolName
93
  TYPE(Variable_t), POINTER :: ThickSol,ThickbSol
94
  REAL(KIND=dp), POINTER :: Thick(:),Thickb(:)
95
  INTEGER, POINTER :: ThickPerm(:),ThickbPerm(:)
96

97
  LOGICAL :: Reset, ComputeDJDsmbTop,ComputeDJDsmbBot,ComputeDJDUV
98
  TYPE(Variable_t), POINTER :: DJDsmbTopSol,DJDsmbBotSol,DJDUVSol
99
  INTEGER,POINTER :: DJDsmbTopPerm(:),DJDsmbBotPerm(:),DJDUVPerm(:)
100
  REAL(KIND=dp), POINTER :: DJDsmbTop(:),DJDsmbBot(:),DJDUV(:)
101
  REAL(KIND=dp), ALLOCATABLE ,SAVE:: LOADb(:),Velob(:,:)
102

103
  !------------------------------------------------------------------------------
104
  !    Go
105
  !------------------------------------------------------------------------------
106
  SolverParams => GetSolverParams()
107

108
  !------------------------------------------------------------------------------
109
  !    check incompatibilities
110
  !------------------------------------------------------------------------------
111
  IF ( CurrentCoordinateSystem() /= Cartesian ) THEN
112
          CALL FATAL(SolverName,'sorry only for cartesian coordinate system')
113
  END IF
114
  IF (TransientSimulation) THEN
115
          CALL FATAL(SolverName,'sorry works only in steady state')
116
  ENDIF
117

118
  !------------------------------------------------------------------------------
119
  !    Allocate some permanent storage, this is done first time only
120
  !------------------------------------------------------------------------------
121
  IF ( .NOT. AllocationsDone .OR. Solver % Mesh % Changed) THEN
122
     NMAX = Model % MaxElementNodes
123

124
     IF ( AllocationsDone ) THEN
125
        DEALLOCATE( ElementNodes % x,    &
126
             ElementNodes % y,    &
127
             ElementNodes % z,    &
128
             FORCE,    &
129
             STIFF, &
130
             LOAD, LOADb,&
131
             Velo,Velob)
132
     END IF
133

134
     ALLOCATE( ElementNodes % x( NMAX ),    &
135
          ElementNodes % y( NMAX ),    &
136
          ElementNodes % z( NMAX ),    &
137
          FORCE( NMAX ),    &
138
          STIFF( NMAX, NMAX ), &
139
          LOAD(NMAX) , LOADb(NMAX),&
140
          Velo( 2, NMAX ), Velob( 2, NMAX ),&
141
          STAT=istat )
142
     IF ( istat /= 0 ) THEN
143
        CALL Fatal(SolverName,'Memory allocation error 1, Aborting.')
144
     END IF
145

146
     CALL Info(SolverName,'Memory allocations done' )
147
     AllocationsDone = .TRUE.
148
  END IF
149

150

151
  !------------------------------------------------------------------------------
152
  !    Get variables 
153
  !------------------------------------------------------------------------------
154
   ThickSolName =  GetString(SolverParams ,'Thickness Solution Name', Found)
155
   IF(.NOT.Found) THEN
156
     CALL WARN(SolverName,'<Thickness Solution Name> not Found; assume default H')
157
     ThickSolName = 'H'
158
   ENDIF
159
   ThickSol => VariableGet( Solver % Mesh % Variables, ThickSolName, UnFoundFatal=.TRUE.)
160
   Thick => ThickSol % Values
161
   ThickPerm => ThickSol % Perm
162

163
   ThickbSolName = GetString(SolverParams ,'Adjoint Solution Name', Found)
164
   IF(.NOT.Found) THEN
165
     CALL WARN(SolverName,'<Adjoint Solution Name> not Found; assume default Adjoint')
166
     ThickSolName = 'Adjoint'
167
   ENDIF
168
   ThickbSol => VariableGet( Solver % Mesh % Variables, ThickbSolName,UnFoundFatal=.TRUE.)
169
   Thickb => ThickbSol % Values
170
   ThickbPerm => ThickbSol % Perm
171

172
   ComputeDJDsmbTop = GetLogical(SolverParams ,'ComputeDJDsmbTop' , Found)
173
   IF (ComputeDJDsmbTop) THEN
174
     DJDsmbTopSol => VariableGet( Solver % Mesh % Variables, 'DJDsmbTop', UnFoundFatal=.TRUE.)
175
     DJDsmbTop => DJDsmbTopSol % Values
176
     DJDsmbTopPerm => DJDsmbTopSol % Perm
177
     Reset =  GetLogical( SolverParams,'Reset DJDsmbTop', Found)
178
     if (Reset.OR.(.NOT.Found)) DJDsmbTop = 0.0
179
   ENDIF
180

181
   ComputeDJDsmbBot = GetLogical(SolverParams ,'ComputeDJDsmbBot' , Found)
182
   IF (ComputeDJDsmbBot) THEN
183
     DJDsmbBotSol => VariableGet( Solver % Mesh % Variables, 'DJDsmbBot',UnFoundFatal=.TRUE. )
184
     DJDsmbBot => DJDsmbBotSol % Values
185
     DJDsmbBotPerm => DJDsmbBotSol % Perm
186
     Reset =  GetLogical( SolverParams,'Reset DJDsmbBot', Found)
187
     if (Reset.OR.(.NOT.Found)) DJDsmbBot = 0.0
188
   ENDIF
189

190
  !------------------------------------------------------------------------------
191
  !    Get Flow solution
192
  !------------------------------------------------------------------------------
193
   ConvectionVar=.True.
194
   FlowSolName =  GetString(SolverParams ,'Flow Solution Name', Found)
195
   IF(.NOT.Found) THEN        
196
        WRITE(Message,'(A)') &
197
           '<Flow Solution Name> Not Found; will look for <convection velocity> in body forces'
198
        CALL Info(SolverName,Message,level=10)
199
        ConvectionVar=.False.
200
        NSDOFS=GetInteger(SolverParams ,'Convection Dimension',Found)
201
        IF(.NOT.Found) &
202
            CALL Fatal(SolverName,'if <Flow Solution Name> not given prescribe <Convection Dimension>')
203
   ELSE
204
        FlowSol => VariableGet( Solver % Mesh % Variables, FlowSolName,UnFoundFatal=.TRUE.)
205
        FlowPerm     => FlowSol % Perm
206
        NSDOFs     =  FlowSol % DOFs
207
        FlowSolution => FlowSol % Values
208
   END IF
209

210
   ComputeDJDUV = GetLogical(SolverParams ,'ComputeDJDUV' , Found)
211
   IF (ComputeDJDUV) THEN
212
     DJDUVSol => VariableGet( Solver % Mesh % Variables, 'DJDUV',UnFoundFatal=.TRUE. )
213
     DJDUV => DJDUVSol % Values
214
     DJDUVPerm => DJDUVSol % Perm
215
     IF (DJDUVSol % DOFs.NE.NSDOFs) &
216
          CALL FATAL(SolverName,'DJDUV DOFs is different from the velocity DOFs')
217
     Reset =  GetLogical( SolverParams,'Reset DJDUV', Found)
218
     if (Reset.OR.(.NOT.Found)) DJDUV = 0.0
219
   ENDIF
220

221
     !------------------------------------------------------------------------------
222
     !    Do the assembly
223
     !------------------------------------------------------------------------------
224
     DO t=1,Solver % NumberOfActiveElements
225

226
        CurrentElement => GetActiveElement(t)
227
        n = GetElementNOFNodes(CurrentElement)
228
        NodeIndexes => CurrentElement % NodeIndexes
229

230
        ! set coords of highest occurring dimension to zero (to get correct path element)
231
        !-------------------------------------------------------------------------------
232
        ElementNodes % x(1:n) = Solver % Mesh % Nodes % x(NodeIndexes)
233
        IF (NSDOFs == 1) THEN
234
           ElementNodes % y(1:n) = 0.0
235
           ElementNodes % z(1:n) = 0.0
236
        ELSE IF (NSDOFs == 2) THEN
237
           ElementNodes % y(1:n) = Solver % Mesh % Nodes % y(NodeIndexes)
238
           ElementNodes % z(1:n) = 0.0_dp
239
        ELSE
240
           WRITE(Message,'(a,i0,a)')&
241
                'It is not possible to compute Thickness evolution if Flow Sol DOFs=',&
242
                NSDOFs, ' . Aborting'
243
           CALL Fatal( SolverName, Message) 
244
           STOP   
245
        END IF
246

247
        ! get pointers on BF
248
        !----------------------------------------------------------------
249
        BodyForce => GetBodyForce(CurrentElement)
250

251
        ! get flow soulution and velocity field from it
252
        !----------------------------------------------
253
        Velo = 0.0_dp
254
        !----------------------------------------------------
255
        ! get velocity profile
256
        IF (ConvectionVar) Then
257
          DO i=1,n
258
             j = NSDOFs*FlowPerm(NodeIndexes(i))
259
              !2D problem - 1D Thickness evolution
260
              IF(NSDOFs == 1) THEN 
261
                 Velo(1,i) = FlowSolution( j ) 
262
                 Velo(2,i) = 0.0_dp
263
              !2D problem - 
264
              ELSE IF (NSDOFs == 2) THEN
265
                 Velo(1,i) = FlowSolution( j-1 ) 
266
                 Velo(2,i) = FlowSolution( j ) 
267
              ELSE
268
                 WRITE(Message,'(a)') 'Velocity should be size 1 or 2'
269
                 CALL Fatal( SolverName, Message)
270
              END IF
271
           END DO
272
       ELSE
273
        IF (ASSOCIATED( BodyForce ) ) THEN
274
          Velo(1,1:n) = ListGetReal( BodyForce, 'Convection Velocity 1',n, NodeIndexes,UnFoundFatal=.TRUE. )
275
          if (NSDOFs.eq.2) &
276
           Velo(2,1:n) = ListGetReal( BodyForce, 'Convection Velocity 2',n, NodeIndexes,UnFoundFatal=.TRUE. )
277
        END IF
278
       END IF
279
        !------------------------------------------------------------------------------
280
        !      get the accumulation/ablation rate (i.e. normal surface flux)
281
        !      from the body force section
282
        !------------------------------------------------------------------------------
283
        LOAD=0.0_dp
284
        IF (ASSOCIATED( BodyForce ) ) THEN
285
              LOAD(1:n) = LOAD(1:n) +   &
286
                      GetReal( BodyForce, 'Top Surface Accumulation', Found )
287
              LOAD(1:n) = LOAD(1:n) +   &
288
                      GetReal( BodyForce, 'Bottom Surface Accumulation', Found )
289
        END IF
290

291

292
        !------------------------------------------------------------------------------
293
        !      Get element local matrix, and rhs vector
294
        !------------------------------------------------------------------------------
295
        CALL LocalMatrix( STIFF, FORCE,LOAD,LOADb,&
296
                Velo, Velob,NSDOFs, &
297
             CurrentElement, n, ElementNodes, NodeIndexes)
298

299
        If (ComputeDJDsmbTop) &
300
           DJDsmbTop(DJDsmbTopPerm(NodeIndexes(1:n)))=DJDsmbTop(DJDsmbTopPerm(NodeIndexes(1:n)))+LOADb(1:n)
301
        If (ComputeDJDsmbBot) &
302
           DJDsmbBot(DJDsmbBotPerm(NodeIndexes(1:n)))=DJDsmbBot(DJDsmbBotPerm(NodeIndexes(1:n)))+LOADb(1:n)
303
        If (ComputeDJDUV) then
304

305
        Do j=1,n
306
          Do i=1,NSDOFs
307
             DJDUV(NSDOFs*(DJDUVPerm(NodeIndexes(j))-1)+i)=DJDUV(NSDOFs*(DJDUVPerm(NodeIndexes(j))-1)+i)+&
308
               Velob(i,j)
309
           End do
310
          End Do
311
        End if
312

313
     END DO ! End loop bulk elements
314

315
     !------------------------------------------------------------------------------
316
   CONTAINS
317

318
     !------------------------------------------------------------------------------
319
     !==============================================================================
320
     SUBROUTINE LocalMatrix( STIFF, FORCE,&
321
          LOAD,LOADb,  Velo, Velob,  NSDOFs, &
322
          Element, n, Nodes, NodeIndexes)
323
       !------------------------------------------------------------------------------
324
       !------------------------------------------------------------------------------
325
       !      external variables:
326
       !      ------------------------------------------------------------------------
327
       REAL(KIND=dp) :: STIFF(:,:),FORCE(:), LOAD(:), Velo(:,:)
328
       REAL(KIND=dp) :: LOADb(:), Velob(:,:)
329

330
       INTEGER :: n, NodeIndexes(:), NSDOFs
331
       TYPE(Nodes_t) :: Nodes
332
       TYPE(Element_t), POINTER :: Element
333

334
       !------------------------------------------------------------------------------
335
       !      internal variables:
336
       !      ------------------------------------------------------------------------
337
       TYPE(GaussIntegrationPoints_t) :: IntegStuff
338
       REAL(KIND=dp) :: Basis(n),dBasisdx(n,3)
339
       REAL(KIND=dp) :: SU(n),SW(n)
340
       REAL(KIND=dp) :: Vgauss(2), UNorm,divu,Source
341
       REAL(KIND=dp) :: U,V,W,S,SqrtElementMetric
342
       REAL(KIND=dp) :: Tau,hk
343
       INTEGER :: i,j,t,p,q
344
       LOGICAL :: stat
345

346
       REAL(KIND=dp) :: Sourceb,divub,Vgaussb(2),SUb(n),SWb(n),Taub,Unormb
347
       !------------------------------------------------------------------------------
348

349
       LOADb = 0.0_dp
350
       Velob = 0.0_dp
351

352

353
       hK = ElementDiameter( Element, Nodes )
354

355
       !
356
       !      Numerical integration:
357
       !      ----------------------
358
       IntegStuff = GaussPoints( Element )
359

360
       SU = 0.0_dp
361
       SW = 0.0_dp
362

363
       DO t = 1,IntegStuff % n
364
          U = IntegStuff % u(t)
365
          V = IntegStuff % v(t)
366
          W = IntegStuff % w(t)
367
          S = IntegStuff % s(t)
368
          !
369
          !        Basis function values & derivatives at the integration point:
370
          !        -------------------------------------------------------------
371
          stat = ElementInfo( Element,Nodes,U,V,W,SqrtElementMetric, &
372
               Basis,dBasisdx, Bubbles=.False.)
373

374
          !        Correction from metric
375
          !        ----------------------
376
          S = S * SqrtElementMetric
377

378
          !
379
          !        Velocities and (norm of) gradient of free surface and source function 
380
          !        at Gauss point
381
          !        ---------------------------------------------------------------------
382

383
          Vgauss=0.0_dp
384
          DO i=1,NSDOFs
385
             Vgauss(i) = SUM( Basis(1:n)*(Velo(i,1:n)))
386
          END DO
387

388
          divu = 0.0_dp
389
          DO i=1,NSDOFs
390
             divu = divu +  SUM( dBasisdx(1:n,i)*(Velo(i,1:n)))
391
          END DO
392

393
          UNorm = SQRT( SUM( Vgauss(1:NSDOFs)**2 ) )
394
          Tau = 0.0_dp
395
          IF (UNorm .GT. 0.0_dp) Tau = hK / ( 2*Unorm )
396

397
          DO p=1,n
398
             SU(p) = 0.0_dp
399
             SW(p) = 0.0_dp
400
             DO i=1,NSDOFs
401
                SU(p) = SU(p) + Vgauss(i) * dBasisdx(p,i)
402
                SW(p) = SW(p) + Vgauss(i) * dBasisdx(p,i)
403
             END DO
404
          END DO
405

406
          !        Stiffness matrix:
407
          !        -----------------
408
          !DO p=1,n
409
          !   DO q=1,n
410
          !      DO i=1,NSDOFs
411
          !         STIFF(p,q) = STIFF(p,q) + &
412
          !              s * Vgauss(i) * dBasisdx(q,i) * Basis(p)
413
          !      END DO
414
          !      STIFF(p,q) =  STIFF(p,q) + s * Tau * SU(q) * SW(p)
415
          !      STIFF(p,q) =  STIFF(p,q) + s * divu * Basis(q) * (Basis(p) + Tau*SW(p))
416
          !   END DO
417
          !END DO
418

419
          !! Adjoint part of the stiffness matrix
420
          divub = 0._dp
421
          Vgaussb = 0._dp
422
          SUb = 0._dp
423
          SWb = 0._dp
424
          Taub =0._dp
425
          Do p=1,n
426
            Do q=1,n
427
               SUb(q) = SUb(q) + s * Tau * SW(p) * &
428
                      (-Thick(ThickPerm(NodeIndexes(q)))*Thickb(ThickbPerm(NodeIndexes(p))))
429
               SWb(p) = SWb(p) + s * Tau * SU(q) * &
430
                      (-Thick(ThickPerm(NodeIndexes(q)))*Thickb(ThickbPerm(NodeIndexes(p))))
431
               SWb(p) = SWb(p) + s * divu * Basis(q) * Tau * &
432
                      (-Thick(ThickPerm(NodeIndexes(q)))*Thickb(ThickbPerm(NodeIndexes(p))))
433
               Taub = Taub + s *  SU(q) * SW(p) * &
434
                      (-Thick(ThickPerm(NodeIndexes(q)))*Thickb(ThickbPerm(NodeIndexes(p))))
435
               Taub = Taub + s * divu * Basis(q) * SW(p) * &
436
                     (-Thick(ThickPerm(NodeIndexes(q)))*Thickb(ThickbPerm(NodeIndexes(p))))
437
               divub = divub + s  * Basis(q) * (Basis(p) + Tau*SW(p)) * &
438
                      (-Thick(ThickPerm(NodeIndexes(q)))*Thickb(ThickbPerm(NodeIndexes(p))))
439
               Do i=1,NSDOFs
440
                 Vgaussb(i)=Vgaussb(i) + s  * dBasisdx(q,i) * Basis(p) * &
441
                     (-Thick(ThickPerm(NodeIndexes(q)))*Thickb(ThickbPerm(NodeIndexes(p))))
442
               End do
443
            End do
444
          End do
445

446
          !        Get accumulation/ablation function 
447
          !        --------------------------------------------------------- 
448
          Source = 0.0_dp
449
          Source=SUM(Basis(1:n)*LOAD(1:n))
450

451
          !        Assemble force vector:
452
          !        ---------------------
453
          !FORCE(1:n) = FORCE(1:n) &
454
          !     + Source * (Basis(1:n) + Tau*SW(1:n)) * s
455
          
456
          !! Adjoint Part of the accumulation/ablation function
457
          Sourceb=0._dp
458
          Do p=1,n
459
            Sourceb = Sourceb + (Basis(p) + Tau*SW(p)) * s * Thickb(ThickbPerm(NodeIndexes(p)))
460
            SWb(p) = SWb(p) + Source * Tau * s * Thickb(ThickbPerm(NodeIndexes(p)))
461
            Taub = Taub + Source * SW(p) * s * Thickb(ThickbPerm(NodeIndexes(p)))
462
          End Do
463
          LOADb(1:n) = LOADb(1:n) + Sourceb * Basis(1:n)
464
          !!
465

466
          !DO p=1,n
467
          !   SU(p) = 0.0_dp
468
          !   SW(p) = 0.0_dp
469
          !   DO i=1,NSDOFs
470
          !      SU(p) = SU(p) + Vgauss(i) * dBasisdx(p,i)
471
          !      SW(p) = SW(p) + Vgauss(i) * dBasisdx(p,i)
472
          !   END DO
473
          !END IF
474

475
          DO p=1,n
476
            DO i=1,NSDOFs
477
              Vgaussb(i)=Vgaussb(i)+SUb(p)*dBasisdx(p,i)
478
              Vgaussb(i)=Vgaussb(i)+SWb(p)*dBasisdx(p,i)
479
            END DO
480
          END DO
481

482
          !UNorm = SQRT( SUM( Vgauss(1:NSDOFs)**2 ) )
483
          !Tau=0._dp
484
          !IF (UNorm .NE. 0.0_dp) Tau = hK / ( 2*Unorm )
485

486
          IF (UNorm .GT. 0.0_dp) THEN
487
             Unormb = -0.5*hK*(Unorm**(-2)) * Taub
488
          ELSE
489
             Unormb = 0._dp
490
          END IF
491

492
          Do i=1,NSDOFs
493
            IF (UNorm .GT. 0.0_dp) &
494
               Vgaussb(i) = Vgaussb(i) + Unormb * 0.5*(SUM( Vgauss(1:NSDOFs)**2)**(-0.5))*2.0*Vgauss(i)
495
          END DO
496

497
          Do p=1,n
498
            Do i=1,NSDOFs
499
              Velob(i,p)=Velob(i,p)+Vgaussb(i)*Basis(p)+divub*dBasisdx(p,i)
500
            End do
501
          End do
502

503

504
       END DO
505

506
       !------------------------------------------------------------------------------
507
     END SUBROUTINE LocalMatrix
508

509
     !------------------------------------------------------------------------------
510
   END SUBROUTINE AdjointThickness_GradientSolver
511
!------------------------------------------------------------------------------
512

513
Product

Resources

Company
