FrontISTR  5.7.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 
6 module m_static_lib_3d_vp
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  do nb = 1,nn
250  do na = 1,nn
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  enddo
299  enddo
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  do i = 1,4
640  do na = 1,nn
641  velo_new( 4*(na-1)+i ) = v(i, na)+dv(i, na)
642  enddo
643  enddo
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  do nb = 1,nn
825  do na = 1,nn
826 
827  mm(na, nb) = spfunc(na)*spfunc(nb)
828  aa(na, nb) = vx*( spfunc(na)*gderiv(nb, 1) ) &
829  +vy*( spfunc(na)*gderiv(nb, 2) ) &
830  +vz*( spfunc(na)*gderiv(nb, 3) )
831  dd(na, nb, 1, 1) = gderiv(na, 1)*gderiv(nb, 1)
832  dd(na, nb, 1, 2) = gderiv(na, 1)*gderiv(nb, 2)
833  dd(na, nb, 1, 3) = gderiv(na, 1)*gderiv(nb, 3)
834  dd(na, nb, 2, 1) = gderiv(na, 2)*gderiv(nb, 1)
835  dd(na, nb, 2, 2) = gderiv(na, 2)*gderiv(nb, 2)
836  dd(na, nb, 2, 3) = gderiv(na, 2)*gderiv(nb, 3)
837  dd(na, nb, 3, 1) = gderiv(na, 3)*gderiv(nb, 1)
838  dd(na, nb, 3, 2) = gderiv(na, 3)*gderiv(nb, 2)
839  dd(na, nb, 3, 3) = gderiv(na, 3)*gderiv(nb, 3)
840  trd(na, nb) = dd(na, nb, 1, 1) &
841  +dd(na, nb, 2, 2) &
842  +dd(na, nb, 3, 3)
843  bb(na, nb) = ( vx*vx )*dd(na, nb, 1, 1) &
844  +( vx*vy )*dd(na, nb, 1, 2) &
845  +( vx*vz )*dd(na, nb, 1, 3) &
846  +( vy*vx )*dd(na, nb, 2, 1) &
847  +( vy*vy )*dd(na, nb, 2, 2) &
848  +( vy*vz )*dd(na, nb, 2, 3) &
849  +( vz*vx )*dd(na, nb, 3, 1) &
850  +( vz*vy )*dd(na, nb, 3, 2) &
851  +( vz*vz )*dd(na, nb, 3, 3)
852  cc(na, nb, 1) = gderiv(na, 1)*spfunc(nb)
853  cc(na, nb, 2) = gderiv(na, 2)*spfunc(nb)
854  cc(na, nb, 3) = gderiv(na, 3)*spfunc(nb)
855 
856  ms(nb, na) = aa(na, nb)
857  as(na, nb) = bb(na, nb)
858  cs(na, nb, 1) = vx*dd(na, nb, 1, 1) &
859  +vy*dd(na, nb, 2, 1) &
860  +vz*dd(na, nb, 3, 1)
861  cs(na, nb, 2) = vx*dd(na, nb, 1, 2) &
862  +vy*dd(na, nb, 2, 2) &
863  +vz*dd(na, nb, 3, 2)
864  cs(na, nb, 3) = vx*dd(na, nb, 1, 3) &
865  +vy*dd(na, nb, 2, 3) &
866  +vz*dd(na, nb, 3, 3)
867  mp(na, nb, 1) = spfunc(nb)*gderiv(na, 1)
868  mp(na, nb, 2) = spfunc(nb)*gderiv(na, 2)
869  mp(na, nb, 3) = spfunc(nb)*gderiv(na, 3)
870  ap(nb, na, 1) = cs(na, nb, 1)
871  ap(nb, na, 2) = cs(na, nb, 2)
872  ap(nb, na, 3) = cs(na, nb, 3)
873  cp(na, nb) = trd(na, nb)
874  enddo
875  enddo
876 
877  !----------------------------------------------------------
878 
879  do nb = 1, nn
880 
881  do na = 1, nn
882 
883  i = 1
884  j = 1
885  isize = 4*(na-1)+i
886  jsize = 4*(nb-1)+j
887 
888  stiff(isize, jsize) &
889  = stiff(isize, jsize) &
890  +wg &
891  *( tincr_inv*rho*( mm(na, nb)+tau*ms(na, nb) ) &
892  +gamma*rho*( aa(na, nb)+tau*as(na, nb) ) &
893  +gamma*mu*( dd(na, nb, i, j)+trd(na, nb) ) )
894 
895  i = 1
896  j = 2
897  isize = 4*(na-1)+i
898  jsize = 4*(nb-1)+j
899 
900  stiff(isize, jsize) &
901  = stiff(isize, jsize) &
902  +wg &
903  *( gamma*mu*dd(na, nb, 2, 1) )
904 
905  i = 1
906  j = 3
907  isize = 4*(na-1)+i
908  jsize = 4*(nb-1)+j
909 
910  stiff(isize, jsize) &
911  = stiff(isize, jsize) &
912  +wg &
913  *( gamma*mu*dd(na, nb, 3, 1) )
914 
915  i = 1
916  j = 4
917  isize = 4*(na-1)+i
918  jsize = 4*(nb-1)+j
919 
920  stiff(isize, jsize) &
921  = stiff(isize, jsize) &
922  +wg &
923  *( -cc(na, nb, 1)+tau*cs(na, nb, 1) )
924 
925  i = 2
926  j = 1
927  isize = 4*(na-1)+i
928  jsize = 4*(nb-1)+j
929 
930  stiff(isize, jsize) &
931  = stiff(isize, jsize) &
932  +wg &
933  *( gamma*mu*dd(na, nb, 1, 2) )
934 
935  i = 2
936  j = 2
937  isize = 4*(na-1)+i
938  jsize = 4*(nb-1)+j
939 
940  stiff(isize, jsize) &
941  = stiff(isize, jsize) &
942  +wg &
943  *( tincr_inv*rho*( mm(na, nb)+tau*ms(na, nb) ) &
944  +gamma*rho*( aa(na, nb)+tau*as(na, nb) ) &
945  +gamma*mu*( dd(na, nb, i, j)+trd(na, nb) ) )
946 
947  i = 2
948  j = 3
949  isize = 4*(na-1)+i
950  jsize = 4*(nb-1)+j
951 
952  stiff(isize, jsize) &
953  = stiff(isize, jsize) &
954  +wg &
955  *( gamma*mu*dd(na, nb, 3, 2) )
956 
957  i = 2
958  j = 4
959  isize = 4*(na-1)+i
960  jsize = 4*(nb-1)+j
961 
962  stiff(isize, jsize) &
963  = stiff(isize, jsize) &
964  +wg &
965  *( -cc(na, nb, 2)+tau*cs(na, nb, 2) )
966 
967  i = 3
968  j = 1
969  isize = 4*(na-1)+i
970  jsize = 4*(nb-1)+j
971 
972  stiff(isize, jsize) &
973  = stiff(isize, jsize) &
974  +wg &
975  *( gamma*mu*dd(na, nb, 1, 3) )
976 
977  i = 3
978  j = 2
979  isize = 4*(na-1)+i
980  jsize = 4*(nb-1)+j
981 
982  stiff(isize, jsize) &
983  = stiff(isize, jsize) &
984  +wg &
985  *( gamma*mu*dd(na, nb, 2, 3) )
986 
987  i = 3
988  j = 3
989  isize = 4*(na-1)+i
990  jsize = 4*(nb-1)+j
991 
992  stiff(isize, jsize) &
993  = stiff(isize, jsize) &
994  +wg &
995  *( tincr_inv*rho*( mm(na, nb)+tau*ms(na, nb) ) &
996  +gamma*rho*( aa(na, nb)+tau*as(na, nb) ) &
997  +gamma*mu*( dd(na, nb, i, j)+trd(na, nb) ) )
998 
999  i = 3
1000  j = 4
1001  isize = 4*(na-1)+i
1002  jsize = 4*(nb-1)+j
1003 
1004  stiff(isize, jsize) &
1005  = stiff(isize, jsize) &
1006  +wg &
1007  *( -cc(na, nb, 3)+tau*cs(na, nb, 3) )
1008 
1009  i = 4
1010  j = 1
1011  isize = 4*(na-1)+i
1012  jsize = 4*(nb-1)+j
1013 
1014  stiff(isize, jsize) &
1015  = stiff(isize, jsize) &
1016  +wg &
1017  *( cc(nb, na, j) &
1018  +tincr_inv*tau*mp(na, nb, j) &
1019  +gamma*tau*ap(na, nb, j) )
1020 
1021  i = 4
1022  j = 2
1023  isize = 4*(na-1)+i
1024  jsize = 4*(nb-1)+j
1025 
1026  stiff(isize, jsize) &
1027  = stiff(isize, jsize) &
1028  +wg &
1029  *( cc(nb, na, j) &
1030  +tincr_inv*tau*mp(na, nb, j) &
1031  +gamma*tau*ap(na, nb, j) )
1032 
1033  i = 4
1034  j = 3
1035  isize = 4*(na-1)+i
1036  jsize = 4*(nb-1)+j
1037 
1038  stiff(isize, jsize) &
1039  = stiff(isize, jsize) &
1040  +wg &
1041  *( cc(nb, na, j) &
1042  +tincr_inv*tau*mp(na, nb, j) &
1043  +gamma*tau*ap(na, nb, j) )
1044 
1045  i = 4
1046  j = 4
1047  isize = 4*(na-1)+i
1048  jsize = 4*(nb-1)+j
1049 
1050  stiff(isize, jsize) &
1051  = stiff(isize, jsize) &
1052  +wg &
1053  *( rho_inv*tau*trd(na, nb) )
1054 
1055  end do
1056 
1057  end do
1058 
1059  !----------------------------------------------------------
1060 
1061  end do loopmatrix
1062 
1063  !--------------------------------------------------------------------
1064 
1065  b(:) = 0.0d0
1066 
1067  loopvector: do lx = 1, numofquadpoints(etype)
1068 
1069  !----------------------------------------------------------
1070 
1071  mu = gausses(lx)%pMaterial%variables(m_viscocity)
1072 
1073  rho = gausses(lx)%pMaterial%variables(m_density)
1074  rho_inv = 1.0d0/rho
1075 
1076  !----------------------------------------------------------
1077 
1078  call getquadpoint( etype, lx, naturalcoord(:) )
1079  call getshapefunc(etype, naturalcoord, spfunc)
1080  call getglobalderiv(etype, nn, naturalcoord, elem, det, gderiv)
1081 
1082  !----------------------------------------------------------
1083 
1084  wg = getweight(etype, lx)*det
1085 
1086  !----------------------------------------------------------
1087 
1088  vx = 0.0d0
1089  vy = 0.0d0
1090  vz = 0.0d0
1091 
1092  do na = 1, nn
1093 
1094  vx = vx+spfunc(na)*v(1, na)
1095  vy = vy+spfunc(na)*v(2, na)
1096  vz = vz+spfunc(na)*v(3, na)
1097 
1098  end do
1099 
1100  !----------------------------------------------------------
1101 
1102  do nb = 1,nn
1103  do na = 1,nn
1104  mm(na, nb) = spfunc(na)*spfunc(nb)
1105  aa(na, nb) = vx*( spfunc(na)*gderiv(nb, 1) ) &
1106  +vy*( spfunc(na)*gderiv(nb, 2) ) &
1107  +vz*( spfunc(na)*gderiv(nb, 3) )
1108  dd(na, nb, 1, 1) = gderiv(na, 1)*gderiv(nb, 1)
1109  dd(na, nb, 1, 2) = gderiv(na, 1)*gderiv(nb, 2)
1110  dd(na, nb, 1, 3) = gderiv(na, 1)*gderiv(nb, 3)
1111  dd(na, nb, 2, 1) = gderiv(na, 2)*gderiv(nb, 1)
1112  dd(na, nb, 2, 2) = gderiv(na, 2)*gderiv(nb, 2)
1113  dd(na, nb, 2, 3) = gderiv(na, 2)*gderiv(nb, 3)
1114  dd(na, nb, 3, 1) = gderiv(na, 3)*gderiv(nb, 1)
1115  dd(na, nb, 3, 2) = gderiv(na, 3)*gderiv(nb, 2)
1116  dd(na, nb, 3, 3) = gderiv(na, 3)*gderiv(nb, 3)
1117  bb(na, nb) = ( vx*vx )*dd(na, nb, 1, 1) &
1118  +( vx*vy )*dd(na, nb, 1, 2) &
1119  +( vx*vz )*dd(na, nb, 1, 3) &
1120  +( vy*vx )*dd(na, nb, 2, 1) &
1121  +( vy*vy )*dd(na, nb, 2, 2) &
1122  +( vy*vz )*dd(na, nb, 2, 3) &
1123  +( vz*vx )*dd(na, nb, 3, 1) &
1124  +( vz*vy )*dd(na, nb, 3, 2) &
1125  +( vz*vz )*dd(na, nb, 3, 3)
1126 
1127  ms(nb, na) = aa(na, nb)
1128  as(na, nb) = bb(na, nb)
1129  cs(na, nb, 1) = vx*dd(na, nb, 1, 1) &
1130  +vy*dd(na, nb, 2, 1) &
1131  +vz*dd(na, nb, 3, 1)
1132  cs(na, nb, 2) = vx*dd(na, nb, 1, 2) &
1133  +vy*dd(na, nb, 2, 2) &
1134  +vz*dd(na, nb, 3, 2)
1135  cs(na, nb, 3) = vx*dd(na, nb, 1, 3) &
1136  +vy*dd(na, nb, 2, 3) &
1137  +vz*dd(na, nb, 3, 3)
1138  mp(na, nb, 1) = spfunc(nb)*gderiv(na, 1)
1139  mp(na, nb, 2) = spfunc(nb)*gderiv(na, 2)
1140  mp(na, nb, 3) = spfunc(nb)*gderiv(na, 3)
1141  ap(nb, na, 1) = cs(na, nb, 1)
1142  ap(nb, na, 2) = cs(na, nb, 2)
1143  ap(nb, na, 3) = cs(na, nb, 3)
1144  enddo
1145  enddo
1146 
1147  !----------------------------------------------------------
1148 
1149  do na = 1, nn
1150 
1151  !----------------------------------------------------
1152 
1153  do i = 1, 3
1154 
1155  m_v(i) = 0.0d0
1156  a_v(i) = 0.0d0
1157  do j = 1, 3
1158  do k = 1, 3
1159  d_v(j, k, i) = 0.0d0
1160  end do
1161  end do
1162  ms_v(i) = 0.0d0
1163  as_v(i) = 0.0d0
1164  mp_dot_v = 0.0d0
1165  ap_dot_v = 0.0d0
1166 
1167  do nb = 1, nn
1168 
1169  ! Unsteady term
1170  m_v(i) = m_v(i)+mm(na, nb)*v(i, nb)
1171  ! Advection term
1172  a_v(i) = a_v(i)+aa(na, nb)*v(i, nb)
1173  ! Diffusion term
1174  do j = 1, 3
1175  do k = 1, 3
1176  d_v(j, k, i) = d_v(j, k, i)+dd(na, nb, j, k)*v(i, nb)
1177  end do
1178  end do
1179  ! Unsteady term (SUPG)
1180  ms_v(i) = ms_v(i)+ms(na, nb)*v(i, nb)
1181  ! Advection term (SUPG)
1182  as_v(i) = as_v(i)+as(na, nb)*v(i, nb)
1183  ! Unsteady term (PSPG)
1184  mp_dot_v = mp_dot_v+( mp(na, nb, 1)*v(1, nb) &
1185  +mp(na, nb, 2)*v(2, nb) &
1186  +mp(na, nb, 3)*v(3, nb) )
1187  ! Advection term (PSPG)
1188  ap_dot_v = ap_dot_v+( ap(na, nb, 1)*v(1, nb) &
1189  +ap(na, nb, 2)*v(2, nb) &
1190  +ap(na, nb, 3)*v(3, nb) )
1191  end do
1192 
1193  end do
1194 
1195  !----------------------------------------------------
1196 
1197  do i = 1, 3
1198 
1199  isize = 4*(na-1)+i
1200 
1201  b(isize) &
1202  = b(isize) &
1203  +wg &
1204  *( tincr_inv*rho*( m_v(i)+tau*ms_v(i) ) &
1205  -( 1.0d0-gamma )*rho*( a_v(i)+tau*as_v(i) ) &
1206  -( 1.0d0-gamma ) &
1207  *mu*( ( d_v(1, 1, i)+d_v(2, 2, i)+d_v(3, 3, i) ) &
1208  +( d_v(1, i, 1)+d_v(2, i, 2)+d_v(3, i, 3) ) ) )
1209 
1210  end do
1211 
1212  i = 4
1213  isize = 4*(na-1)+i
1214 
1215  b(isize) &
1216  = b(isize) &
1217  +wg &
1218  *( tincr_inv*tau*( mp_dot_v ) &
1219  -( 1.0d0-gamma )*tau*( ap_dot_v ) )
1220 
1221  !----------------------------------------------------
1222 
1223  end do
1224 
1225  !----------------------------------------------------------
1226 
1227  end do loopvector
1228 
1229  !----------------------------------------------------------------
1230 
1231  do isize = 1, 4*nn
1232 
1233  stiff_velo = 0.0d0
1234 
1235  do jsize = 1, 4*nn
1236 
1237  stiff_velo = stiff_velo+stiff(isize, jsize)*velo_new(jsize)
1238 
1239  end do
1240 
1241  r(isize) = b(isize)-stiff_velo
1242 
1243  end do
1244 
1245  !--------------------------------------------------------------------
1246  end subroutine load_c3_vp
1247  !--------------------------------------------------------------------
1248 
1249 
1250 end module m_static_lib_3d_vp
elementinfo::getweight
real(kind=kreal) function getweight(fetype, np)
Fetch the weight value in given gauss point.
Definition: element.f90:535
mcontact::mu
real(kind=kreal), save mu
penalty, default value
Definition: fstr_contact.f90:17
elementinfo::numofquadpoints
integer function numofquadpoints(fetype)
Obtains the number of quadrature points of the element.
Definition: element.f90:450
m_static_lib_3d_vp::stf_c3_vp
subroutine stf_c3_vp(etype, nn, ecoord, gausses, stiff, tincr, v, temperature)
Definition: static_LIB_3d_vp.f90:18
m_static_lib_3d_vp
Definition: static_LIB_3d_vp.f90:6
mmechgauss
This modules defines a structure to record history dependent parameter in static analysis.
Definition: mechgauss.f90:6
elementinfo::getglobalderiv
subroutine getglobalderiv(fetype, nn, localcoord, elecoord, det, gderiv)
Calculate shape derivative in global coordinate system.
Definition: element.f90:741
m_static_lib_3d_vp::update_c3_vp
subroutine update_c3_vp(etype, nn, ecoord, v, dv, gausses)
Definition: static_LIB_3d_vp.f90:495
elementinfo
This module encapsulate the basic functions of all elements provide by this software.
Definition: element.f90:43
elementinfo::getquadpoint
subroutine getquadpoint(fetype, np, pos)
Fetch the coordinate of gauss point.
Definition: element.f90:489
hecmw
Definition: hecmw.f90:6
elementinfo::getshapefunc
subroutine getshapefunc(fetype, localcoord, func)
Calculate the shape function in natural coordinate system.
Definition: element.f90:647
mmechgauss::tgaussstatus
All data should be recorded in every quadrature points.
Definition: mechgauss.f90:13
m_static_lib_3d_vp::load_c3_vp
subroutine load_c3_vp(etype, nn, ecoord, v, dv, r, gausses, tincr)
Definition: static_LIB_3d_vp.f90:583