FrontISTR  5.2.0
Large-scale structural analysis program with finit element method
static_LIB_3d_vp.f90
Go to the documentation of this file.
1 !-------------------------------------------------------------------------------
2 ! Copyright (c) 2019 FrontISTR Commons
3 ! This software is released under the MIT License, see LICENSE.txt
4 !-------------------------------------------------------------------------------
5 
7 
8  use hecmw, only : kint, kreal
9  use elementinfo
10 
11  implicit none
12 
13 contains
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
25  integer(kind=kint), intent(in) :: nn
26  real(kind=kreal), intent(in) :: ecoord(3, nn)
27  type(tgaussstatus), intent(in) :: gausses(:)
28  real(kind=kreal), intent(out) :: stiff(:,:)
29  real(kind=kreal), intent(in) :: tincr
30  real(kind=kreal), intent(in), optional :: v(:, :)
31  real(kind=kreal), intent(in), optional :: temperature(nn)
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
502  integer(kind=kint), intent(in) :: nn
503  real(kind=kreal), intent(in) :: ecoord(3, nn)
504  real(kind=kreal), intent(in) :: v(4, nn)
505  real(kind=kreal), intent(in) :: dv(4, nn)
506  type(tgaussstatus), intent(inout) :: gausses(:)
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
590  integer(kind=kint), intent(in) :: nn
591  real(kind=kreal), intent(in) :: ecoord(3, nn)
592  real(kind=kreal), intent(in) :: v(4, nn)
593  real(kind=kreal), intent(in) :: dv(4, nn)
594  real(kind=kreal), intent(out) :: r(4*nn)
595  type(tgaussstatus), intent(inout) :: gausses(:)
596  real(kind=kreal), intent(in) :: tincr
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 
1249 end module m_static_lib_3d_vp
This module encapsulate the basic functions of all elements provide by this software.
Definition: element.f90:43
subroutine getshapefunc(fetype, localcoord, func)
Calculate the shape function in natural coodinate system.
Definition: element.f90:640
subroutine getquadpoint(fetype, np, pos)
Fetch the coordinate of gauss point.
Definition: element.f90:482
subroutine getglobalderiv(fetype, nn, localcoord, elecoord, det, gderiv)
Calculate shape derivative in global coordinate system.
Definition: element.f90:734
real(kind=kreal) function getweight(fetype, np)
Fetch the weight value in given gauss point.
Definition: element.f90:528
integer function numofquadpoints(fetype)
Obtains the number of quadrature points of the element.
Definition: element.f90:445
Definition: hecmw.f90:6
subroutine stf_c3_vp(etype, nn, ecoord, gausses, stiff, tincr, v, temperature)
subroutine update_c3_vp(etype, nn, ecoord, v, dv, gausses)
subroutine load_c3_vp(etype, nn, ecoord, v, dv, r, gausses, tincr)
This modules defines a structure to record history dependent parameter in static analysis.
Definition: mechgauss.f90:6
All data should be recorded in every quadrature points.
Definition: mechgauss.f90:13