FrontISTR 5.2.0
Large-scale structural analysis program with finit element method
Loading...
Searching...
No Matches
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
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
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
1249end 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