8 use hecmw,
only : kint, kreal
16 (etype, nn, ecoord, gausses, stiff, tincr, &
24 integer(kind=kint),
intent(in) :: etype
25 integer(kind=kint),
intent(in) :: nn
26 real(kind=kreal),
intent(in) :: ecoord(3, nn)
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)
35 integer(kind=kint) :: i, j
36 integer(kind=kint) :: na, nb
37 integer(kind=kint) :: isize, jsize
38 integer(kind=kint) :: lx
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
56 real(kind=kreal) :: det, wg
57 real(kind=kreal) :: tau
58 real(kind=kreal),
parameter :: gamma = 0.5d0
62 tincr_inv = 1.0d0/tincr
66 elem(:, :) = ecoord(:, :)
96 volume_inv = 1.0d0/volume
100 naturalcoord(1) = 0.25d0
101 naturalcoord(2) = 0.25d0
102 naturalcoord(3) = 0.25d0
112 vx = vx+spfunc(na)*v(1, na)
113 vy = vy+spfunc(na)*v(2, na)
114 vz = vz+spfunc(na)*v(3, na)
118 v_dot_v = vx*vx+vy*vy+vz*vz
147 mu = mu+wg*gausses(lx)%pMaterial%variables(m_viscocity)
149 rho = rho+wg*gausses(lx)%pMaterial%variables(m_density)
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)
163 end do loopglobalderiv
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)
182 d = d+dabs( vx*dndx(na, 1)+vy*dndx(na, 2)+vz*dndx(na, 3) )
195 if( v_dot_v .LT. 1.0d-15 )
then
197 t3 = 4.0d0*mu/( rho*volume**(2.0d0/3.0d0) )
201 t3 = mu*d*d/( rho*v_dot_v )
207 tau = 1.0d0/dsqrt( t1*t1+t2*t2+t3*t3 )
217 mu = gausses(lx)%pMaterial%variables(m_viscocity)
219 rho = gausses(lx)%pMaterial%variables(m_density)
241 vx = vx+spfunc(na)*v(1, na)
242 vy = vy+spfunc(na)*v(2, na)
243 vz = vz+spfunc(na)*v(3, na)
249 forall(na = 1:nn, nb = 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) &
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)
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) &
285 cs(na, nb, 2) = vx*dd(na, nb, 1, 2) &
286 +vy*dd(na, nb, 2, 2) &
288 cs(na, nb, 3) = vx*dd(na, nb, 1, 3) &
289 +vy*dd(na, nb, 2, 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)
312 stiff(isize, jsize) &
313 = stiff(isize, jsize) &
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) ) )
324 stiff(isize, jsize) &
325 = stiff(isize, jsize) &
327 *( gamma*mu*dd(na, nb, 2, 1) )
334 stiff(isize, jsize) &
335 = stiff(isize, jsize) &
337 *( gamma*mu*dd(na, nb, 3, 1) )
344 stiff(isize, jsize) &
345 = stiff(isize, jsize) &
347 *( -cc(na, nb, 1)+tau*cs(na, nb, 1) )
354 stiff(isize, jsize) &
355 = stiff(isize, jsize) &
357 *( gamma*mu*dd(na, nb, 1, 2) )
364 stiff(isize, jsize) &
365 = stiff(isize, jsize) &
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) ) )
376 stiff(isize, jsize) &
377 = stiff(isize, jsize) &
379 *( gamma*mu*dd(na, nb, 3, 2) )
386 stiff(isize, jsize) &
387 = stiff(isize, jsize) &
389 *( -cc(na, nb, 2)+tau*cs(na, nb, 2) )
396 stiff(isize, jsize) &
397 = stiff(isize, jsize) &
399 *( gamma*mu*dd(na, nb, 1, 3) )
406 stiff(isize, jsize) &
407 = stiff(isize, jsize) &
409 *( gamma*mu*dd(na, nb, 2, 3) )
416 stiff(isize, jsize) &
417 = stiff(isize, jsize) &
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) ) )
428 stiff(isize, jsize) &
429 = stiff(isize, jsize) &
431 *( -cc(na, nb, 3)+tau*cs(na, nb, 3) )
438 stiff(isize, jsize) &
439 = stiff(isize, jsize) &
442 +tincr_inv*tau*mp(na, nb, j) &
443 +gamma*tau*ap(na, nb, j) )
450 stiff(isize, jsize) &
451 = stiff(isize, jsize) &
454 +tincr_inv*tau*mp(na, nb, j) &
455 +gamma*tau*ap(na, nb, j) )
462 stiff(isize, jsize) &
463 = stiff(isize, jsize) &
466 +tincr_inv*tau*mp(na, nb, j) &
467 +gamma*tau*ap(na, nb, j) )
474 stiff(isize, jsize) &
475 = stiff(isize, jsize) &
477 *( rho_inv*tau*trd(na, nb) )
494 (etype, nn, ecoord, v, dv, gausses)
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)
510 integer(kind=kint) :: lx
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
523 elem(:, :) = ecoord(:, :)
525 totalvelo(:, :) = v(:, :)+dv(:, :)
533 mu = gausses(lx)%pMaterial%variables(m_viscocity)
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) )
555 p = dot_product(totalvelo(4, 1:nn), spfunc(1:nn))
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)
568 gausses(lx)%strain_out(1:6) = gausses(lx)%strain(1:6)
569 gausses(lx)%stress_out(1:6) = gausses(lx)%stress(1:6)
582 (etype, nn, ecoord, v, dv, r, gausses, tincr)
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)
596 real(kind=kreal),
intent(in) :: tincr
600 integer(kind=kint) :: i, j, k
601 integer(kind=kint) :: na, nb
602 integer(kind=kint) :: isize, jsize
603 integer(kind=kint) :: lx
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
633 tincr_inv = 1.0d0/tincr
637 elem(:, :) = ecoord(:, :)
639 forall(na = 1:nn, i = 1:4)
641 velo_new( 4*(na-1)+i ) = v(i, na)+dv(i, na)
673 volume_inv = 1.0d0/volume
677 naturalcoord(1) = 0.25d0
678 naturalcoord(2) = 0.25d0
679 naturalcoord(3) = 0.25d0
689 vx = vx+spfunc(na)*v(1, na)
690 vy = vy+spfunc(na)*v(2, na)
691 vz = vz+spfunc(na)*v(3, na)
695 v_a_dot_v_a = vx*vx+vy*vy+vz*vz
724 mu = mu +wg*gausses(lx)%pMaterial%variables(m_viscocity)
725 rho = rho+wg*gausses(lx)%pMaterial%variables(m_density)
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)
739 end do loopglobalderiv
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)
758 d = d+dabs( vx*dndx(na, 1)+vy*dndx(na, 2)+vz*dndx(na, 3) )
771 if( v_a_dot_v_a .LT. 1.0d-15 )
then
773 t3 = 4.0d0*mu/( rho*volume**(2.0d0/3.0d0) )
777 t3 = mu*d*d/( rho*v_a_dot_v_a )
783 tau = 1.0d0/dsqrt( t1*t1+t2*t2+t3*t3 )
793 mu = gausses(lx)%pMaterial%variables(m_viscocity)
795 rho = gausses(lx)%pMaterial%variables(m_density)
816 vx = vx+spfunc(na)*v(1, na)
817 vy = vy+spfunc(na)*v(2, na)
818 vz = vz+spfunc(na)*v(3, na)
824 forall(na = 1:nn, nb = 1:nn)
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) &
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)
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) &
860 cs(na, nb, 2) = vx*dd(na, nb, 1, 2) &
861 +vy*dd(na, nb, 2, 2) &
863 cs(na, nb, 3) = vx*dd(na, nb, 1, 3) &
864 +vy*dd(na, nb, 2, 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)
887 stiff(isize, jsize) &
888 = stiff(isize, jsize) &
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) ) )
899 stiff(isize, jsize) &
900 = stiff(isize, jsize) &
902 *( gamma*mu*dd(na, nb, 2, 1) )
909 stiff(isize, jsize) &
910 = stiff(isize, jsize) &
912 *( gamma*mu*dd(na, nb, 3, 1) )
919 stiff(isize, jsize) &
920 = stiff(isize, jsize) &
922 *( -cc(na, nb, 1)+tau*cs(na, nb, 1) )
929 stiff(isize, jsize) &
930 = stiff(isize, jsize) &
932 *( gamma*mu*dd(na, nb, 1, 2) )
939 stiff(isize, jsize) &
940 = stiff(isize, jsize) &
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) ) )
951 stiff(isize, jsize) &
952 = stiff(isize, jsize) &
954 *( gamma*mu*dd(na, nb, 3, 2) )
961 stiff(isize, jsize) &
962 = stiff(isize, jsize) &
964 *( -cc(na, nb, 2)+tau*cs(na, nb, 2) )
971 stiff(isize, jsize) &
972 = stiff(isize, jsize) &
974 *( gamma*mu*dd(na, nb, 1, 3) )
981 stiff(isize, jsize) &
982 = stiff(isize, jsize) &
984 *( gamma*mu*dd(na, nb, 2, 3) )
991 stiff(isize, jsize) &
992 = stiff(isize, jsize) &
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) ) )
1003 stiff(isize, jsize) &
1004 = stiff(isize, jsize) &
1006 *( -cc(na, nb, 3)+tau*cs(na, nb, 3) )
1013 stiff(isize, jsize) &
1014 = stiff(isize, jsize) &
1017 +tincr_inv*tau*mp(na, nb, j) &
1018 +gamma*tau*ap(na, nb, j) )
1025 stiff(isize, jsize) &
1026 = stiff(isize, jsize) &
1029 +tincr_inv*tau*mp(na, nb, j) &
1030 +gamma*tau*ap(na, nb, j) )
1037 stiff(isize, jsize) &
1038 = stiff(isize, jsize) &
1041 +tincr_inv*tau*mp(na, nb, j) &
1042 +gamma*tau*ap(na, nb, j) )
1049 stiff(isize, jsize) &
1050 = stiff(isize, jsize) &
1052 *( rho_inv*tau*trd(na, nb) )
1070 mu = gausses(lx)%pMaterial%variables(m_viscocity)
1072 rho = gausses(lx)%pMaterial%variables(m_density)
1093 vx = vx+spfunc(na)*v(1, na)
1094 vy = vy+spfunc(na)*v(2, na)
1095 vz = vz+spfunc(na)*v(3, na)
1101 forall(na = 1:nn, nb = 1:nn)
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)
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)
1158 d_v(j, k, i) = 0.0d0
1169 m_v(i) = m_v(i)+mm(na, nb)*v(i, nb)
1171 a_v(i) = a_v(i)+aa(na, nb)*v(i, nb)
1175 d_v(j, k, i) = d_v(j, k, i)+dd(na, nb, j, k)*v(i, nb)
1179 ms_v(i) = ms_v(i)+ms(na, nb)*v(i, nb)
1181 as_v(i) = as_v(i)+as(na, nb)*v(i, nb)
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) )
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) )
1203 *( tincr_inv*rho*( m_v(i)+tau*ms_v(i) ) &
1204 -( 1.0d0-gamma )*rho*( a_v(i)+tau*as_v(i) ) &
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) ) ) )
1217 *( tincr_inv*tau*( mp_dot_v ) &
1218 -( 1.0d0-gamma )*tau*( ap_dot_v ) )
1236 stiff_velo = stiff_velo+stiff(isize, jsize)*velo_new(jsize)
1240 r(isize) = b(isize)-stiff_velo
This module encapsulate the basic functions of all elements provide by this software.
subroutine getshapefunc(fetype, localcoord, func)
Calculate the shape function in natural coodinate system.
subroutine getquadpoint(fetype, np, pos)
Fetch the coordinate of gauss point.
subroutine getglobalderiv(fetype, nn, localcoord, elecoord, det, gderiv)
Calculate shape derivative in global coordinate system.
real(kind=kreal) function getweight(fetype, np)
Fetch the weight value in given gauss point.
integer function numofquadpoints(fetype)
Obtains the number of quadrature points of the element.
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.
All data should be recorded in every quadrature points.