1!-------------------------------------------------------------------------------
2! Copyright (c) 2019 FrontISTR Commons
3! This software is released under the MIT License, see LICENSE.txt
4!-------------------------------------------------------------------------------
5
6module m_static_LIB_3d_vp
7
8  use hecmw, only : kint, kreal
9  use elementInfo
10
11  implicit none
12
13contains
14  !--------------------------------------------------------------------
15  subroutine STF_C3_vp                                  &
16      (etype, nn, ecoord, gausses, stiff, tincr, &
17      v, temperature)
18    !--------------------------------------------------------------------
19
20    use mMechGauss
21
22    !--------------------------------------------------------------------
23
24    integer(kind=kint), intent(in) :: etype                     !< element type
25    integer(kind=kint), intent(in) :: nn                        !< the number of elemental nodes
26    real(kind=kreal),   intent(in) :: ecoord(3, nn)             !< coordinates of elemental nodes
27    type(tGaussStatus), intent(in) :: gausses(:)                !< status of qudrature points
28    real(kind=kreal),   intent(out) :: stiff(:,:)               !< stiff matrix
29    real(kind=kreal),   intent(in) :: tincr                     !< time increment
30    real(kind=kreal),   intent(in), optional :: v(:, :)         !< nodal velocity
31    real(kind=kreal),   intent(in), optional :: temperature(nn) !< temperature
32
33    !--------------------------------------------------------------------
34
35    integer(kind=kint) :: i, j
36    integer(kind=kint) :: na, nb
37    integer(kind=kint) :: isize, jsize
38    integer(kind=kint) :: LX
39
40    real(kind=kreal) :: MM(nn, nn), AA(nn, nn), DD(nn, nn, 3, 3), &
41      trD(nn, nn), BB(nn, nn), CC(nn, nn, 3),   &
42      MS(nn, nn), AS(nn, nn), CS(nn, nn, 3),    &
43      MP(nn, nn, 3), AP(nn, nn, 3), CP(nn, nn)
44    real(kind=kreal) :: spfunc(nn), gderiv(nn, 3)
45    real(kind=kreal) :: elem(3, nn)
46    real(kind=kreal) :: naturalCoord(3)
47    real(kind=kreal) :: dndx(nn, 3)
48    real(kind=kreal) :: tincr_inv
49    real(kind=kreal) :: volume, volume_inv
50    real(kind=kreal) :: mu
51    real(kind=kreal) :: rho, rho_inv
52    real(kind=kreal) :: vx, vy, vz
53    real(kind=kreal) :: t1, t2, t3
54    real(kind=kreal) :: v_dot_v
55    real(kind=kreal) :: d
56    real(kind=kreal) :: det, wg
57    real(kind=kreal) :: tau
58    real(kind=kreal), parameter :: gamma = 0.5D0
59
60    !--------------------------------------------------------------------
61
62    tincr_inv = 1.0D0/tincr
63
64    !--------------------------------------------------------------------
65
66    elem(:, :) = ecoord(:, :)
67
68    !--------------------------------------------------------------------
69
70    t1 = 2.0D0*tincr_inv
71
72    !---------------------------------------------------------------
73
74    volume = 0.0D0
75
76    loopVolume: do LX = 1, NumOfQuadPoints(etype)
77
78      !----------------------------------------------------------
79
80      call getQuadPoint( etype, LX, naturalCoord(:) )
81      call getShapeFunc(etype, naturalCoord, spfunc)
82      call getGlobalDeriv(etype, nn, naturalCoord, elem, det, gderiv)
83
84      !----------------------------------------------------------
85
86      wg = getWeight(etype, LX)*det
87
88      !----------------------------------------------------------
89
90      volume = volume+wg
91
92      !----------------------------------------------------------
93
94    end do loopVolume
95
96    volume_inv = 1.0D0/volume
97
98    !---------------------------------------------------------------
99
100    naturalCoord(1) = 0.25D0
101    naturalCoord(2) = 0.25D0
102    naturalCoord(3) = 0.25D0
103
104    call getShapeFunc(etype, naturalCoord, spfunc)
105
106    vx = 0.0D0
107    vy = 0.0D0
108    vz = 0.0D0
109
110    do na = 1, nn
111
112      vx = vx+spfunc(na)*v(1, na)
113      vy = vy+spfunc(na)*v(2, na)
114      vz = vz+spfunc(na)*v(3, na)
115
116    end do
117
118    v_dot_v = vx*vx+vy*vy+vz*vz
119
120    !---------------------------------------------------------------
121
122    mu  = 0.0D0
123    rho = 0.0D0
124
125    do na = 1, nn
126
127      dndx(na, 1) = 0.0D0
128      dndx(na, 2) = 0.0D0
129      dndx(na, 3) = 0.0D0
130
131    end do
132
133    loopGlobalDeriv: do LX = 1, NumOfQuadPoints(etype)
134
135      !----------------------------------------------------------
136
137      call getQuadPoint( etype, LX, naturalCoord(1:3) )
138      call getShapeFunc(etype, naturalCoord, spfunc)
139      call getGlobalDeriv(etype, nn, naturalCoord, elem, det, gderiv)
140
141      !----------------------------------------------------------
142
143      wg = getWeight(etype, LX)*det
144
145      !----------------------------------------------------------
146
147      mu = mu+wg*gausses(LX)%pMaterial%variables(M_VISCOCITY)
148
149      rho = rho+wg*gausses(LX)%pMaterial%variables(M_DENSITY)
150
151      !----------------------------------------------------------
152
153      do na = 1, nn
154
155        dndx(na, 1) = dndx(na, 1)+wg*gderiv(na, 1)
156        dndx(na, 2) = dndx(na, 2)+wg*gderiv(na, 2)
157        dndx(na, 3) = dndx(na, 3)+wg*gderiv(na, 3)
158
159      end do
160
161      !----------------------------------------------------------
162
163    end do loopGlobalDeriv
164
165    mu  = volume_inv*mu
166    rho = volume_inv*rho
167
168    do na = 1, nn
169
170      dndx(na, 1) = volume_inv*dndx(na, 1)
171      dndx(na, 2) = volume_inv*dndx(na, 2)
172      dndx(na, 3) = volume_inv*dndx(na, 3)
173
174    end do
175
176    !---------------------------------------------------------------
177
178    d = 0.0D0
179
180    do na = 1, nn
181
182      d = d+dabs( vx*dndx(na, 1)+vy*dndx(na, 2)+vz*dndx(na, 3) )
183
184    end do
185
186    ! h_es3d = 2.0D0/( d/DSQRT( v_dot_v ) )
187
188    !---------------------------------------------------------------
189
190    ! t2 = 2.0D0*DSQRT( v_dot_v )/h_es3d
191    t2 = d
192
193    !----------------------------------------------------------------
194
195    if( v_dot_v .LT. 1.0D-15 ) then
196
197      t3 = 4.0D0*mu/( rho*volume**(2.0D0/3.0D0) )
198
199    else
200
201      t3 = mu*d*d/( rho*v_dot_v )
202
203    end if
204
205    !----------------------------------------------------------------
206
207    tau = 1.0D0/dsqrt( t1*t1+t2*t2+t3*t3 )
208
209    !--------------------------------------------------------------------
210
211    stiff(:, :) = 0.0D0
212
213    loopGauss: do LX = 1, NumOfQuadPoints(etype)
214
215      !----------------------------------------------------------
216
217      mu = gausses(LX)%pMaterial%variables(M_VISCOCITY)
218
219      rho = gausses(LX)%pMaterial%variables(M_DENSITY)
220      rho_inv = 1.0D0/rho
221
222      !----------------------------------------------------------
223
224      call getQuadPoint( etype, LX, naturalCoord(1:3) )
225      call getShapeFunc( etype, naturalCoord(:), spfunc(:) )
226      call getGlobalDeriv( etype, nn, naturalCoord(:), elem(:,:), &
227        det, gderiv(:,:) )
228
229      !----------------------------------------------------------
230
231      wg = getWeight(etype, LX)*det
232
233      !----------------------------------------------------------
234
235      vx = 0.0D0
236      vy = 0.0D0
237      vz = 0.0D0
238
239      do na = 1, nn
240
241        vx = vx+spfunc(na)*v(1, na)
242        vy = vy+spfunc(na)*v(2, na)
243        vz = vz+spfunc(na)*v(3, na)
244
245      end do
246
247      !----------------------------------------------------------
248
249      forall(na = 1:nn, nb = 1:nn)
250
251        MM(na, nb) = spfunc(na)*spfunc(nb)
252        AA(na, nb) = vx*( spfunc(na)*gderiv(nb, 1) ) &
253          +vy*( spfunc(na)*gderiv(nb, 2) ) &
254          +vz*( spfunc(na)*gderiv(nb, 3) )
255        DD(na, nb, 1, 1) = gderiv(na, 1)*gderiv(nb, 1)
256        DD(na, nb, 1, 2) = gderiv(na, 1)*gderiv(nb, 2)
257        DD(na, nb, 1, 3) = gderiv(na, 1)*gderiv(nb, 3)
258        DD(na, nb, 2, 1) = gderiv(na, 2)*gderiv(nb, 1)
259        DD(na, nb, 2, 2) = gderiv(na, 2)*gderiv(nb, 2)
260        DD(na, nb, 2, 3) = gderiv(na, 2)*gderiv(nb, 3)
261        DD(na, nb, 3, 1) = gderiv(na, 3)*gderiv(nb, 1)
262        DD(na, nb, 3, 2) = gderiv(na, 3)*gderiv(nb, 2)
263        DD(na, nb, 3, 3) = gderiv(na, 3)*gderiv(nb, 3)
264        trD(na, nb) = DD(na, nb, 1, 1) &
265          +DD(na, nb, 2, 2) &
266          +DD(na, nb, 3, 3)
267        BB(na, nb) = ( vx*vx )*DD(na, nb, 1, 1) &
268          +( vx*vy )*DD(na, nb, 1, 2) &
269          +( vx*vz )*DD(na, nb, 1, 3) &
270          +( vy*vx )*DD(na, nb, 2, 1) &
271          +( vy*vy )*DD(na, nb, 2, 2) &
272          +( vy*vz )*DD(na, nb, 2, 3) &
273          +( vz*vx )*DD(na, nb, 3, 1) &
274          +( vz*vy )*DD(na, nb, 3, 2) &
275          +( vz*vz )*DD(na, nb, 3, 3)
276        CC(na, nb, 1) = gderiv(na, 1)*spfunc(nb)
277        CC(na, nb, 2) = gderiv(na, 2)*spfunc(nb)
278        CC(na, nb, 3) = gderiv(na, 3)*spfunc(nb)
279
280        MS(nb, na) = AA(na, nb)
281        AS(na, nb) = BB(na, nb)
282        CS(na, nb, 1) = vx*DD(na, nb, 1, 1) &
283          +vy*DD(na, nb, 2, 1) &
284          +vz*DD(na, nb, 3, 1)
285        CS(na, nb, 2) = vx*DD(na, nb, 1, 2) &
286          +vy*DD(na, nb, 2, 2) &
287          +vz*DD(na, nb, 3, 2)
288        CS(na, nb, 3) = vx*DD(na, nb, 1, 3) &
289          +vy*DD(na, nb, 2, 3) &
290          +vz*DD(na, nb, 3, 3)
291        MP(na, nb, 1) = spfunc(nb)*gderiv(na, 1)
292        MP(na, nb, 2) = spfunc(nb)*gderiv(na, 2)
293        MP(na, nb, 3) = spfunc(nb)*gderiv(na, 3)
294        AP(nb, na, 1) = CS(na, nb, 1)
295        AP(nb, na, 2) = CS(na, nb, 2)
296        AP(nb, na, 3) = CS(na, nb, 3)
297        CP(na, nb) = trD(na, nb)
298
299      end forall
300
301      !----------------------------------------------------------
302
303      do nb = 1, nn
304
305        do na = 1, nn
306
307          i = 1
308          j = 1
309          isize = 4*(na-1)+i
310          jsize = 4*(nb-1)+j
311
312          stiff(isize, jsize)                              &
313            = stiff(isize, jsize)                            &
314            +wg                                             &
315            *( tincr_inv*rho*( MM(na, nb)+tau*MS(na, nb) ) &
316            +gamma*rho*( AA(na, nb)+tau*AS(na, nb) )     &
317            +gamma*mu*( DD(na, nb, i, j)+trD(na, nb) ) )
318
319          i = 1
320          j = 2
321          isize = 4*(na-1)+i
322          jsize = 4*(nb-1)+j
323
324          stiff(isize, jsize)              &
325            = stiff(isize, jsize)            &
326            +wg                             &
327            *( gamma*mu*DD(na, nb, 2, 1) )
328
329          i = 1
330          j = 3
331          isize = 4*(na-1)+i
332          jsize = 4*(nb-1)+j
333
334          stiff(isize, jsize)              &
335            = stiff(isize, jsize)            &
336            +wg                             &
337            *( gamma*mu*DD(na, nb, 3, 1) )
338
339          i = 1
340          j = 4
341          isize = 4*(na-1)+i
342          jsize = 4*(nb-1)+j
343
344          stiff(isize, jsize)                     &
345            = stiff(isize, jsize)                   &
346            +wg                                    &
347            *( -CC(na, nb, 1)+tau*CS(na, nb, 1) )
348
349          i = 2
350          j = 1
351          isize = 4*(na-1)+i
352          jsize = 4*(nb-1)+j
353
354          stiff(isize, jsize)              &
355            = stiff(isize, jsize)            &
356            +wg                             &
357            *( gamma*mu*DD(na, nb, 1, 2) )
358
359          i = 2
360          j = 2
361          isize = 4*(na-1)+i
362          jsize = 4*(nb-1)+j
363
364          stiff(isize, jsize)                              &
365            = stiff(isize, jsize)                            &
366            +wg                                             &
367            *( tincr_inv*rho*( MM(na, nb)+tau*MS(na, nb) ) &
368            +gamma*rho*( AA(na, nb)+tau*AS(na, nb) )     &
369            +gamma*mu*( DD(na, nb, i, j)+trD(na, nb) ) )
370
371          i = 2
372          j = 3
373          isize = 4*(na-1)+i
374          jsize = 4*(nb-1)+j
375
376          stiff(isize, jsize)              &
377            = stiff(isize, jsize)            &
378            +wg                             &
379            *( gamma*mu*DD(na, nb, 3, 2) )
380
381          i = 2
382          j = 4
383          isize = 4*(na-1)+i
384          jsize = 4*(nb-1)+j
385
386          stiff(isize, jsize)                     &
387            = stiff(isize, jsize)                   &
388            +wg                                    &
389            *( -CC(na, nb, 2)+tau*CS(na, nb, 2) )
390
391          i = 3
392          j = 1
393          isize = 4*(na-1)+i
394          jsize = 4*(nb-1)+j
395
396          stiff(isize, jsize)              &
397            = stiff(isize, jsize)            &
398            +wg                             &
399            *( gamma*mu*DD(na, nb, 1, 3) )
400
401          i = 3
402          j = 2
403          isize = 4*(na-1)+i
404          jsize = 4*(nb-1)+j
405
406          stiff(isize, jsize)              &
407            = stiff(isize, jsize)            &
408            +wg                             &
409            *( gamma*mu*DD(na, nb, 2, 3) )
410
411          i = 3
412          j = 3
413          isize = 4*(na-1)+i
414          jsize = 4*(nb-1)+j
415
416          stiff(isize, jsize)                               &
417            = stiff(isize, jsize)                             &
418            +wg                                              &
419            *( tincr_inv*rho*( MM(na, nb)+tau*MS(na, nb) )  &
420            +gamma*rho*( AA(na, nb)+tau*AS(na, nb) )      &
421            +gamma*mu*( DD(na, nb, i, j)+trD(na, nb) ) )
422
423          i = 3
424          j = 4
425          isize = 4*(na-1)+i
426          jsize = 4*(nb-1)+j
427
428          stiff(isize, jsize)                     &
429            = stiff(isize, jsize)                   &
430            +wg                                    &
431            *( -CC(na, nb, 3)+tau*CS(na, nb, 3) )
432
433          i = 4
434          j = 1
435          isize = 4*(na-1)+i
436          jsize = 4*(nb-1)+j
437
438          stiff(isize, jsize)              &
439            = stiff(isize, jsize)            &
440            +wg                             &
441            *( CC(nb, na, j)               &
442            +tincr_inv*tau*MP(na, nb, j) &
443            +gamma*tau*AP(na, nb, j) )
444
445          i = 4
446          j = 2
447          isize = 4*(na-1)+i
448          jsize = 4*(nb-1)+j
449
450          stiff(isize, jsize)              &
451            = stiff(isize, jsize)            &
452            +wg                             &
453            *( CC(nb, na, j)               &
454            +tincr_inv*tau*MP(na, nb, j) &
455            +gamma*tau*AP(na, nb, j) )
456
457          i = 4
458          j = 3
459          isize = 4*(na-1)+i
460          jsize = 4*(nb-1)+j
461
462          stiff(isize, jsize)              &
463            = stiff(isize, jsize)            &
464            +wg                             &
465            *( CC(nb, na, j)               &
466            +tincr_inv*tau*MP(na, nb, j) &
467            +gamma*tau*AP(na, nb, j) )
468
469          i = 4
470          j = 4
471          isize = 4*(na-1)+i
472          jsize = 4*(nb-1)+j
473
474          stiff(isize, jsize)            &
475            = stiff(isize, jsize)          &
476            +wg                           &
477            *( rho_inv*tau*trD(na, nb) )
478
479        end do
480
481      end do
482
483      !----------------------------------------------------------
484
485    end do loopGauss
486
487    !--------------------------------------------------------------------
488  end subroutine STF_C3_vp
489  !--------------------------------------------------------------------
490
491
492  !--------------------------------------------------------------------
493  subroutine UPDATE_C3_vp                        &
494      (etype, nn, ecoord, v, dv, gausses)
495    !--------------------------------------------------------------------
496
497    use mMechGauss
498
499    !--------------------------------------------------------------------
500
501    integer(kind=kint), intent(in) :: etype            !< element type
502    integer(kind=kint), intent(in) :: nn               !< the number of elemental nodes
503    real(kind=kreal), intent(in) :: ecoord(3, nn)      !< coordinates of elemental nodes
504    real(kind=kreal), intent(in) :: v(4, nn)           !< nodal velcoity
505    real(kind=kreal), intent(in) :: dv(4, nn)          !< nodal velocity increment
506    type(tGaussStatus), intent(inout) :: gausses(:)    !< status of qudrature points
507
508    !--------------------------------------------------------------------
509
510    integer(kind=kint) :: LX
511
512    real(kind=kreal) :: elem(3, nn)
513    real(kind=kreal) :: totalvelo(4, nn)
514    real(kind=kreal) :: spfunc(nn), gderiv(nn, 3)
515    real(kind=kreal) :: gveloderiv(3, 3)
516    real(kind=kreal) :: naturalCoord(3)
517    real(kind=kreal) :: det
518    real(kind=kreal) :: mu
519    real(kind=kreal) :: p
520
521    !--------------------------------------------------------------------
522
523    elem(:, :) = ecoord(:, :)
524
525    totalvelo(:, :) = v(:, :)+dv(:, :)
526
527    !--------------------------------------------------------------------
528
529    loopMatrix: do LX = 1, NumOfQuadPoints(etype)
530
531      !----------------------------------------------------------
532
533      mu = gausses(LX)%pMaterial%variables(M_VISCOCITY)
534
535      !----------------------------------------------------------
536
537      call getQuadPoint( etype, LX, naturalCoord(:) )
538      call getShapeFunc(etype, naturalCoord, spfunc)
539      call getGlobalDeriv(etype, nn, naturalCoord, elem, det, gderiv)
540
541      !----------------------------------------------------------
542
543      ! Deformation rate tensor
544      gveloderiv(1:3, 1:3) = matmul( totalvelo(1:3, 1:nn), gderiv(1:nn, 1:3) )
545      gausses(LX)%strain(1) = gveloderiv(1, 1)
546      gausses(LX)%strain(2) = gveloderiv(2, 2)
547      gausses(LX)%strain(3) = gveloderiv(3, 3)
548      gausses(LX)%strain(4) = 0.5D0*( gveloderiv(1, 2)+gveloderiv(2, 1) )
549      gausses(LX)%strain(5) = 0.5D0*( gveloderiv(2, 3)+gveloderiv(3, 2) )
550      gausses(LX)%strain(6) = 0.5D0*( gveloderiv(3, 1)+gveloderiv(1, 3) )
551
552      !----------------------------------------------------------
553
554      ! Pressure
555      p = dot_product(totalvelo(4, 1:nn), spfunc(1:nn))
556
557      ! Cauchy stress tensor
558      gausses(LX)%stress(1) = -p+2.0D0*mu*gausses(LX)%strain(1)
559      gausses(LX)%stress(2) = -p+2.0D0*mu*gausses(LX)%strain(2)
560      gausses(LX)%stress(3) = -p+2.0D0*mu*gausses(LX)%strain(3)
561      gausses(LX)%stress(4) = 2.0D0*mu*gausses(LX)%strain(4)
562      gausses(LX)%stress(5) = 2.0D0*mu*gausses(LX)%strain(5)
563      gausses(LX)%stress(6) = 2.0D0*mu*gausses(LX)%strain(6)
564
565      !----------------------------------------------------------
566
567      !set stress and strain for output
568      gausses(LX)%strain_out(1:6) = gausses(LX)%strain(1:6)
569      gausses(LX)%stress_out(1:6) = gausses(LX)%stress(1:6)
570
571    end do loopMatrix
572
573    !----------------------------------------------------------------
574
575    !--------------------------------------------------------------------
576  end subroutine UPDATE_C3_vp
577  !--------------------------------------------------------------------
578
579
580  !--------------------------------------------------------------------
581  subroutine LOAD_C3_vp                                    &
582      (etype, nn, ecoord, v, dv, r, gausses, tincr)
583    !--------------------------------------------------------------------
584
585    use mMechGauss
586
587    !--------------------------------------------------------------------
588
589    integer(kind=kint), intent(in)    :: etype         !< element type
590    integer(kind=kint), intent(in)    :: nn            !< the number of elemental nodes
591    real(kind=kreal), intent(in)      :: ecoord(3, nn) !< coordinates of elemental nodes
592    real(kind=kreal), intent(in)      :: v(4, nn)      !< nodal dislplacements
593    real(kind=kreal), intent(in)      :: dv(4, nn)     !< nodal velocity increment
594    real(kind=kreal), intent(out)     :: r(4*nn)       !< elemental residual
595    type(tGaussStatus), intent(inout) :: gausses(:)    !< status of qudrature points
596    real(kind=kreal), intent(in)      :: tincr         !< time increment
597
598    !--------------------------------------------------------------------
599
600    integer(kind=kint) :: i, j, k
601    integer(kind=kint) :: na, nb
602    integer(kind=kint) :: isize, jsize
603    integer(kind=kint) :: LX
604
605    real(kind=kreal) :: elem(3, nn)
606    real(kind=kreal) :: velo_new(4*nn)
607    real(kind=kreal) :: stiff(4*nn, 4*nn)
608    real(kind=kreal) :: b(4*nn)
609    real(kind=kreal) :: MM(nn, nn), AA(nn, nn), DD(nn, nn, 3, 3), &
610      trD(nn, nn), BB(nn, nn), CC(nn, nn, 3),   &
611      MS(nn, nn), AS(nn, nn), CS(nn, nn, 3),    &
612      MP(nn, nn, 3), AP(nn, nn, 3), CP(nn, nn)
613    real(kind=kreal) :: spfunc(nn), gderiv(nn, 3)
614    real(kind=kreal) :: naturalCoord(3)
615    real(kind=kreal) :: dndx(nn, 3)
616    real(kind=kreal) :: tincr_inv
617    real(kind=kreal) :: volume, volume_inv
618    real(kind=kreal) :: mu
619    real(kind=kreal) :: rho, rho_inv
620    real(kind=kreal) :: vx, vy, vz
621    real(kind=kreal) :: t1, t2, t3
622    real(kind=kreal) :: v_a_dot_v_a
623    real(kind=kreal) :: d
624    real(kind=kreal) :: det, wg
625    real(kind=kreal) :: tau
626    real(kind=kreal) :: m_v(3), a_v(3), d_v(3, 3, 3),        &
627      ms_v(3), as_v(3), mp_dot_v, ap_dot_v
628    real(kind=kreal) :: stiff_velo
629    real(kind=kreal), parameter :: gamma = 0.5D0
630
631    !--------------------------------------------------------------------
632
633    tincr_inv = 1.0D0/tincr
634
635    !--------------------------------------------------------------------
636
637    elem(:, :) = ecoord(:, :)
638
639    forall(na = 1:nn, i = 1:4)
640
641      velo_new( 4*(na-1)+i ) = v(i, na)+dv(i, na)
642
643    end forall
644
645    !--------------------------------------------------------------------
646
647    t1 = 2.0D0*tincr_inv
648
649    !---------------------------------------------------------------
650
651    volume = 0.0D0
652
653    loopVolume: do LX = 1, NumOfQuadPoints(etype)
654
655      !----------------------------------------------------------
656
657      call getQuadPoint( etype, LX, naturalCoord(:) )
658      call getShapeFunc(etype, naturalCoord, spfunc)
659      call getGlobalDeriv(etype, nn, naturalCoord, elem, det, gderiv)
660
661      !----------------------------------------------------------
662
663      wg = getWeight(etype, LX)*det
664
665      !----------------------------------------------------------
666
667      volume = volume+wg
668
669      !----------------------------------------------------------
670
671    end do loopVolume
672
673    volume_inv = 1.0D0/volume
674
675    !---------------------------------------------------------------
676
677    naturalCoord(1) = 0.25D0
678    naturalCoord(2) = 0.25D0
679    naturalCoord(3) = 0.25D0
680
681    call getShapeFunc(etype, naturalCoord, spfunc)
682
683    vx = 0.0D0
684    vy = 0.0D0
685    vz = 0.0D0
686
687    do na = 1, nn
688
689      vx = vx+spfunc(na)*v(1, na)
690      vy = vy+spfunc(na)*v(2, na)
691      vz = vz+spfunc(na)*v(3, na)
692
693    end do
694
695    v_a_dot_v_a = vx*vx+vy*vy+vz*vz
696
697    !---------------------------------------------------------------
698
699    do na = 1, nn
700
701      dndx(na, 1) = 0.0D0
702      dndx(na, 2) = 0.0D0
703      dndx(na, 3) = 0.0D0
704
705    end do
706
707    mu = 0.0d0
708    rho = 0.0d0
709
710    loopGlobalDeriv: do LX = 1, NumOfQuadPoints(etype)
711
712      !----------------------------------------------------------
713
714      call getQuadPoint( etype, LX, naturalCoord(1:3) )
715      call getShapeFunc(etype, naturalCoord, spfunc)
716      call getGlobalDeriv(etype, nn, naturalCoord, elem, det, gderiv)
717
718      !----------------------------------------------------------
719
720      wg = getWeight(etype, LX)*det
721
722      !----------------------------------------------------------
723
724      mu  = mu +wg*gausses(LX)%pMaterial%variables(M_VISCOCITY)
725      rho = rho+wg*gausses(LX)%pMaterial%variables(M_DENSITY)
726
727      !----------------------------------------------------------
728
729      do na = 1, nn
730
731        dndx(na, 1) = dndx(na, 1)+wg*gderiv(na, 1)
732        dndx(na, 2) = dndx(na, 2)+wg*gderiv(na, 2)
733        dndx(na, 3) = dndx(na, 3)+wg*gderiv(na, 3)
734
735      end do
736
737      !----------------------------------------------------------
738
739    end do loopGlobalDeriv
740
741    mu  = volume_inv*mu
742    rho = volume_inv*rho
743
744    do na = 1, nn
745
746      dndx(na, 1) = volume_inv*dndx(na, 1)
747      dndx(na, 2) = volume_inv*dndx(na, 2)
748      dndx(na, 3) = volume_inv*dndx(na, 3)
749
750    end do
751
752    !---------------------------------------------------------------
753
754    d = 0.0D0
755
756    do na = 1, nn
757
758      d = d+dabs( vx*dndx(na, 1)+vy*dndx(na, 2)+vz*dndx(na, 3) )
759
760    end do
761
762    ! h_es3d = 2.0D0/( d/DSQRT( v_dot_v ) )
763
764    !---------------------------------------------------------------
765
766    ! t2 = 2.0D0*DSQRT( v_dot_v )/h_es3d
767    t2 = d
768
769    !----------------------------------------------------------------
770
771    if( v_a_dot_v_a .LT. 1.0D-15 ) then
772
773      t3 = 4.0D0*mu/( rho*volume**(2.0D0/3.0D0) )
774
775    else
776
777      t3 = mu*d*d/( rho*v_a_dot_v_a )
778
779    end if
780
781    !----------------------------------------------------------------
782
783    tau = 1.0D0/dsqrt( t1*t1+t2*t2+t3*t3 )
784
785    !--------------------------------------------------------------------
786
787    stiff(:, :) = 0.0D0
788
789    loopMatrix: do LX = 1, NumOfQuadPoints(etype)
790
791      !----------------------------------------------------------
792
793      mu = gausses(LX)%pMaterial%variables(M_VISCOCITY)
794
795      rho = gausses(LX)%pMaterial%variables(M_DENSITY)
796      rho_inv = 1.0D0/rho
797
798      !----------------------------------------------------------
799
800      call getQuadPoint( etype, LX, naturalCoord(:) )
801      call getShapeFunc(etype, naturalCoord, spfunc)
802      call getGlobalDeriv(etype, nn, naturalCoord, elem, det, gderiv)
803
804      !----------------------------------------------------------
805
806      wg = getWeight(etype, LX)*det
807
808      !----------------------------------------------------------
809
810      vx = 0.0D0
811      vy = 0.0D0
812      vz = 0.0D0
813
814      do na = 1, nn
815
816        vx = vx+spfunc(na)*v(1, na)
817        vy = vy+spfunc(na)*v(2, na)
818        vz = vz+spfunc(na)*v(3, na)
819
820      end do
821
822      !----------------------------------------------------------
823
824      forall(na = 1:nn, nb = 1:nn)
825
826        MM(na, nb) = spfunc(na)*spfunc(nb)
827        AA(na, nb) = vx*( spfunc(na)*gderiv(nb, 1) ) &
828          +vy*( spfunc(na)*gderiv(nb, 2) ) &
829          +vz*( spfunc(na)*gderiv(nb, 3) )
830        DD(na, nb, 1, 1) = gderiv(na, 1)*gderiv(nb, 1)
831        DD(na, nb, 1, 2) = gderiv(na, 1)*gderiv(nb, 2)
832        DD(na, nb, 1, 3) = gderiv(na, 1)*gderiv(nb, 3)
833        DD(na, nb, 2, 1) = gderiv(na, 2)*gderiv(nb, 1)
834        DD(na, nb, 2, 2) = gderiv(na, 2)*gderiv(nb, 2)
835        DD(na, nb, 2, 3) = gderiv(na, 2)*gderiv(nb, 3)
836        DD(na, nb, 3, 1) = gderiv(na, 3)*gderiv(nb, 1)
837        DD(na, nb, 3, 2) = gderiv(na, 3)*gderiv(nb, 2)
838        DD(na, nb, 3, 3) = gderiv(na, 3)*gderiv(nb, 3)
839        trD(na, nb) = DD(na, nb, 1, 1) &
840          +DD(na, nb, 2, 2) &
841          +DD(na, nb, 3, 3)
842        BB(na, nb) = ( vx*vx )*DD(na, nb, 1, 1) &
843          +( vx*vy )*DD(na, nb, 1, 2) &
844          +( vx*vz )*DD(na, nb, 1, 3) &
845          +( vy*vx )*DD(na, nb, 2, 1) &
846          +( vy*vy )*DD(na, nb, 2, 2) &
847          +( vy*vz )*DD(na, nb, 2, 3) &
848          +( vz*vx )*DD(na, nb, 3, 1) &
849          +( vz*vy )*DD(na, nb, 3, 2) &
850          +( vz*vz )*DD(na, nb, 3, 3)
851        CC(na, nb, 1) = gderiv(na, 1)*spfunc(nb)
852        CC(na, nb, 2) = gderiv(na, 2)*spfunc(nb)
853        CC(na, nb, 3) = gderiv(na, 3)*spfunc(nb)
854
855        MS(nb, na) = AA(na, nb)
856        AS(na, nb) = BB(na, nb)
857        CS(na, nb, 1) = vx*DD(na, nb, 1, 1) &
858          +vy*DD(na, nb, 2, 1) &
859          +vz*DD(na, nb, 3, 1)
860        CS(na, nb, 2) = vx*DD(na, nb, 1, 2) &
861          +vy*DD(na, nb, 2, 2) &
862          +vz*DD(na, nb, 3, 2)
863        CS(na, nb, 3) = vx*DD(na, nb, 1, 3) &
864          +vy*DD(na, nb, 2, 3) &
865          +vz*DD(na, nb, 3, 3)
866        MP(na, nb, 1) = spfunc(nb)*gderiv(na, 1)
867        MP(na, nb, 2) = spfunc(nb)*gderiv(na, 2)
868        MP(na, nb, 3) = spfunc(nb)*gderiv(na, 3)
869        AP(nb, na, 1) = CS(na, nb, 1)
870        AP(nb, na, 2) = CS(na, nb, 2)
871        AP(nb, na, 3) = CS(na, nb, 3)
872        CP(na, nb) = trD(na, nb)
873
874      end forall
875
876      !----------------------------------------------------------
877
878      do nb = 1, nn
879
880        do na = 1, nn
881
882          i = 1
883          j = 1
884          isize = 4*(na-1)+i
885          jsize = 4*(nb-1)+j
886
887          stiff(isize, jsize)                              &
888            = stiff(isize, jsize)                            &
889            +wg                                             &
890            *( tincr_inv*rho*( MM(na, nb)+tau*MS(na, nb) ) &
891            +gamma*rho*( AA(na, nb)+tau*AS(na, nb) )     &
892            +gamma*mu*( DD(na, nb, i, j)+trD(na, nb) ) )
893
894          i = 1
895          j = 2
896          isize = 4*(na-1)+i
897          jsize = 4*(nb-1)+j
898
899          stiff(isize, jsize)              &
900            = stiff(isize, jsize)            &
901            +wg                             &
902            *( gamma*mu*DD(na, nb, 2, 1) )
903
904          i = 1
905          j = 3
906          isize = 4*(na-1)+i
907          jsize = 4*(nb-1)+j
908
909          stiff(isize, jsize)              &
910            = stiff(isize, jsize)            &
911            +wg                             &
912            *( gamma*mu*DD(na, nb, 3, 1) )
913
914          i = 1
915          j = 4
916          isize = 4*(na-1)+i
917          jsize = 4*(nb-1)+j
918
919          stiff(isize, jsize)                     &
920            = stiff(isize, jsize)                   &
921            +wg                                    &
922            *( -CC(na, nb, 1)+tau*CS(na, nb, 1) )
923
924          i = 2
925          j = 1
926          isize = 4*(na-1)+i
927          jsize = 4*(nb-1)+j
928
929          stiff(isize, jsize)              &
930            = stiff(isize, jsize)            &
931            +wg                             &
932            *( gamma*mu*DD(na, nb, 1, 2) )
933
934          i = 2
935          j = 2
936          isize = 4*(na-1)+i
937          jsize = 4*(nb-1)+j
938
939          stiff(isize, jsize)                              &
940            = stiff(isize, jsize)                            &
941            +wg                                             &
942            *( tincr_inv*rho*( MM(na, nb)+tau*MS(na, nb) ) &
943            +gamma*rho*( AA(na, nb)+tau*AS(na, nb) )     &
944            +gamma*mu*( DD(na, nb, i, j)+trD(na, nb) ) )
945
946          i = 2
947          j = 3
948          isize = 4*(na-1)+i
949          jsize = 4*(nb-1)+j
950
951          stiff(isize, jsize)              &
952            = stiff(isize, jsize)            &
953            +wg                             &
954            *( gamma*mu*DD(na, nb, 3, 2) )
955
956          i = 2
957          j = 4
958          isize = 4*(na-1)+i
959          jsize = 4*(nb-1)+j
960
961          stiff(isize, jsize)                     &
962            = stiff(isize, jsize)                   &
963            +wg                                    &
964            *( -CC(na, nb, 2)+tau*CS(na, nb, 2) )
965
966          i = 3
967          j = 1
968          isize = 4*(na-1)+i
969          jsize = 4*(nb-1)+j
970
971          stiff(isize, jsize)              &
972            = stiff(isize, jsize)            &
973            +wg                             &
974            *( gamma*mu*DD(na, nb, 1, 3) )
975
976          i = 3
977          j = 2
978          isize = 4*(na-1)+i
979          jsize = 4*(nb-1)+j
980
981          stiff(isize, jsize)              &
982            = stiff(isize, jsize)            &
983            +wg                             &
984            *( gamma*mu*DD(na, nb, 2, 3) )
985
986          i = 3
987          j = 3
988          isize = 4*(na-1)+i
989          jsize = 4*(nb-1)+j
990
991          stiff(isize, jsize)                              &
992            = stiff(isize, jsize)                            &
993            +wg                                             &
994            *( tincr_inv*rho*( MM(na, nb)+tau*MS(na, nb) ) &
995            +gamma*rho*( AA(na, nb)+tau*AS(na, nb) )     &
996            +gamma*mu*( DD(na, nb, i, j)+trD(na, nb) ) )
997
998          i = 3
999          j = 4
1000          isize = 4*(na-1)+i
1001          jsize = 4*(nb-1)+j
1002
1003          stiff(isize, jsize)                     &
1004            = stiff(isize, jsize)                   &
1005            +wg                                    &
1006            *( -CC(na, nb, 3)+tau*CS(na, nb, 3) )
1007
1008          i = 4
1009          j = 1
1010          isize = 4*(na-1)+i
1011          jsize = 4*(nb-1)+j
1012
1013          stiff(isize, jsize)              &
1014            = stiff(isize, jsize)            &
1015            +wg                             &
1016            *( CC(nb, na, j)               &
1017            +tincr_inv*tau*MP(na, nb, j) &
1018            +gamma*tau*AP(na, nb, j) )
1019
1020          i = 4
1021          j = 2
1022          isize = 4*(na-1)+i
1023          jsize = 4*(nb-1)+j
1024
1025          stiff(isize, jsize)              &
1026            = stiff(isize, jsize)            &
1027            +wg                             &
1028            *( CC(nb, na, j)               &
1029            +tincr_inv*tau*MP(na, nb, j) &
1030            +gamma*tau*AP(na, nb, j) )
1031
1032          i = 4
1033          j = 3
1034          isize = 4*(na-1)+i
1035          jsize = 4*(nb-1)+j
1036
1037          stiff(isize, jsize)              &
1038            = stiff(isize, jsize)            &
1039            +wg                             &
1040            *( CC(nb, na, j)               &
1041            +tincr_inv*tau*MP(na, nb, j) &
1042            +gamma*tau*AP(na, nb, j) )
1043
1044          i = 4
1045          j = 4
1046          isize = 4*(na-1)+i
1047          jsize = 4*(nb-1)+j
1048
1049          stiff(isize, jsize)            &
1050            = stiff(isize, jsize)          &
1051            +wg                           &
1052            *( rho_inv*tau*trD(na, nb) )
1053
1054        end do
1055
1056      end do
1057
1058      !----------------------------------------------------------
1059
1060    end do loopMatrix
1061
1062    !--------------------------------------------------------------------
1063
1064    b(:) = 0.0D0
1065
1066    loopVector: do LX = 1, NumOfQuadPoints(etype)
1067
1068      !----------------------------------------------------------
1069
1070      mu = gausses(LX)%pMaterial%variables(M_VISCOCITY)
1071
1072      rho = gausses(LX)%pMaterial%variables(M_DENSITY)
1073      rho_inv = 1.0D0/rho
1074
1075      !----------------------------------------------------------
1076
1077      call getQuadPoint( etype, LX, naturalCoord(:) )
1078      call getShapeFunc(etype, naturalCoord, spfunc)
1079      call getGlobalDeriv(etype, nn, naturalCoord, elem, det, gderiv)
1080
1081      !----------------------------------------------------------
1082
1083      wg = getWeight(etype, LX)*det
1084
1085      !----------------------------------------------------------
1086
1087      vx = 0.0D0
1088      vy = 0.0D0
1089      vz = 0.0D0
1090
1091      do na = 1, nn
1092
1093        vx = vx+spfunc(na)*v(1, na)
1094        vy = vy+spfunc(na)*v(2, na)
1095        vz = vz+spfunc(na)*v(3, na)
1096
1097      end do
1098
1099      !----------------------------------------------------------
1100
1101      forall(na = 1:nn, nb = 1:nn)
1102
1103        MM(na, nb) = spfunc(na)*spfunc(nb)
1104        AA(na, nb) = vx*( spfunc(na)*gderiv(nb, 1) ) &
1105          +vy*( spfunc(na)*gderiv(nb, 2) ) &
1106          +vz*( spfunc(na)*gderiv(nb, 3) )
1107        DD(na, nb, 1, 1) = gderiv(na, 1)*gderiv(nb, 1)
1108        DD(na, nb, 1, 2) = gderiv(na, 1)*gderiv(nb, 2)
1109        DD(na, nb, 1, 3) = gderiv(na, 1)*gderiv(nb, 3)
1110        DD(na, nb, 2, 1) = gderiv(na, 2)*gderiv(nb, 1)
1111        DD(na, nb, 2, 2) = gderiv(na, 2)*gderiv(nb, 2)
1112        DD(na, nb, 2, 3) = gderiv(na, 2)*gderiv(nb, 3)
1113        DD(na, nb, 3, 1) = gderiv(na, 3)*gderiv(nb, 1)
1114        DD(na, nb, 3, 2) = gderiv(na, 3)*gderiv(nb, 2)
1115        DD(na, nb, 3, 3) = gderiv(na, 3)*gderiv(nb, 3)
1116        BB(na, nb) = ( vx*vx )*DD(na, nb, 1, 1) &
1117          +( vx*vy )*DD(na, nb, 1, 2) &
1118          +( vx*vz )*DD(na, nb, 1, 3) &
1119          +( vy*vx )*DD(na, nb, 2, 1) &
1120          +( vy*vy )*DD(na, nb, 2, 2) &
1121          +( vy*vz )*DD(na, nb, 2, 3) &
1122          +( vz*vx )*DD(na, nb, 3, 1) &
1123          +( vz*vy )*DD(na, nb, 3, 2) &
1124          +( vz*vz )*DD(na, nb, 3, 3)
1125
1126        MS(nb, na) = AA(na, nb)
1127        AS(na, nb) = BB(na, nb)
1128        CS(na, nb, 1) = vx*DD(na, nb, 1, 1) &
1129          +vy*DD(na, nb, 2, 1) &
1130          +vz*DD(na, nb, 3, 1)
1131        CS(na, nb, 2) = vx*DD(na, nb, 1, 2) &
1132          +vy*DD(na, nb, 2, 2) &
1133          +vz*DD(na, nb, 3, 2)
1134        CS(na, nb, 3) = vx*DD(na, nb, 1, 3) &
1135          +vy*DD(na, nb, 2, 3) &
1136          +vz*DD(na, nb, 3, 3)
1137        MP(na, nb, 1) = spfunc(nb)*gderiv(na, 1)
1138        MP(na, nb, 2) = spfunc(nb)*gderiv(na, 2)
1139        MP(na, nb, 3) = spfunc(nb)*gderiv(na, 3)
1140        AP(nb, na, 1) = CS(na, nb, 1)
1141        AP(nb, na, 2) = CS(na, nb, 2)
1142        AP(nb, na, 3) = CS(na, nb, 3)
1143
1144      end forall
1145
1146      !----------------------------------------------------------
1147
1148      do na = 1, nn
1149
1150        !----------------------------------------------------
1151
1152        do i = 1, 3
1153
1154          m_v(i) = 0.0D0
1155          a_v(i) = 0.0D0
1156          do j = 1, 3
1157            do k = 1, 3
1158              d_v(j, k, i) = 0.0D0
1159            end do
1160          end do
1161          ms_v(i) = 0.0D0
1162          as_v(i) = 0.0D0
1163          mp_dot_v = 0.0D0
1164          ap_dot_v = 0.0D0
1165
1166          do nb = 1, nn
1167
1168            ! Unsteady term
1169            m_v(i) = m_v(i)+MM(na, nb)*v(i, nb)
1170            ! Advection term
1171            a_v(i) = a_v(i)+AA(na, nb)*v(i, nb)
1172            ! Diffusion term
1173            do j = 1, 3
1174              do k = 1, 3
1175                d_v(j, k, i) = d_v(j, k, i)+DD(na, nb, j, k)*v(i, nb)
1176              end do
1177            end do
1178            ! Unsteady term (SUPG)
1179            ms_v(i) = ms_v(i)+MS(na, nb)*v(i, nb)
1180            ! Advection term (SUPG)
1181            as_v(i) = as_v(i)+AS(na, nb)*v(i, nb)
1182            ! Unsteady term (PSPG)
1183            mp_dot_v = mp_dot_v+( MP(na, nb, 1)*v(1, nb)   &
1184              +MP(na, nb, 2)*v(2, nb)   &
1185              +MP(na, nb, 3)*v(3, nb) )
1186            ! Advection term (PSPG)
1187            ap_dot_v = ap_dot_v+( AP(na, nb, 1)*v(1, nb)   &
1188              +AP(na, nb, 2)*v(2, nb)   &
1189              +AP(na, nb, 3)*v(3, nb) )
1190          end do
1191
1192        end do
1193
1194        !----------------------------------------------------
1195
1196        do i = 1, 3
1197
1198          isize = 4*(na-1)+i
1199
1200          b(isize)                                                  &
1201            = b(isize)                                                &
1202            +wg                                                      &
1203            *( tincr_inv*rho*( m_v(i)+tau*ms_v(i) )                 &
1204            -( 1.0D0-gamma )*rho*( a_v(i)+tau*as_v(i) )           &
1205            -( 1.0D0-gamma )                                      &
1206            *mu*( ( d_v(1, 1, i)+d_v(2, 2, i)+d_v(3, 3, i) )     &
1207            +( d_v(1, i, 1)+d_v(2, i, 2)+d_v(3, i, 3) ) ) )
1208
1209        end do
1210
1211        i = 4
1212        isize = 4*(na-1)+i
1213
1214        b(isize)                                &
1215          = b(isize)                              &
1216          +wg                                    &
1217          *( tincr_inv*tau*( mp_dot_v )         &
1218          -( 1.0D0-gamma )*tau*( ap_dot_v ) )
1219
1220        !----------------------------------------------------
1221
1222      end do
1223
1224      !----------------------------------------------------------
1225
1226    end do loopVector
1227
1228    !----------------------------------------------------------------
1229
1230    do isize = 1, 4*nn
1231
1232      stiff_velo = 0.0D0
1233
1234      do jsize = 1, 4*nn
1235
1236        stiff_velo = stiff_velo+stiff(isize, jsize)*velo_new(jsize)
1237
1238      end do
1239
1240      r(isize) = b(isize)-stiff_velo
1241
1242    end do
1243
1244    !--------------------------------------------------------------------
1245  end subroutine LOAD_C3_vp
1246  !--------------------------------------------------------------------
1247
1248
1249end module m_static_LIB_3d_vp
1250