40
41
42
43 use matparam_def_mod
44
45
46
47#include "implicit_f.inc"
48
49
50
51#include "mvsiz_p.inc"
52#include "param_c.inc"
53
54
55
56 INTEGER, INTENT(IN) :: NEL
57 INTEGER, INTENT(IN) :: NUMNOD
58 INTEGER NC1(MVSIZ),NC2(MVSIZ),NC3(MVSIZ),NC4(MVSIZ),
59 . NC5(MVSIZ),NC6(MVSIZ),NC7(MVSIZ),NC8(MVSIZ)
60 INTEGER IPARTS(*)
62 my_real ,
INTENT(IN) :: theaccfact
63 my_real ,
INTENT(IN) :: heat(mvsiz)
70 my_real :: px1(*), px2(*), px3(*), px4(*),
71 . py1(*), py2(*), py3(*), py4(*),
72 . pz1(*), pz2(*), pz3(*), pz4(*)
75 type (matparam_struct_) ,intent(in) :: mat_param
76
77
78
79 INTEGER I,M
80 my_real as, bs, kc, phix, phiy, phiz, a, b, rhocp, t0
81
82 as = mat_param%THERM%AS
83 bs = mat_param%THERM%BS
84 rhocp = mat_param%THERM%RHOCP
85 t0 = mat_param%THERM%TINI
86
87
88
89 DO i=1,nel
90 IF(off(i)==zero.OR.offg(i)<=zero) cycle
91 phix = tempnc(nc1(i))*px1(i) + tempnc(nc2(i))*px2(i) +
92 . tempnc(nc3(i))*px3(i) + tempnc(nc4(i))*px4(i) -
93 . tempnc(nc5(i))*px3(i) - tempnc(nc6(i))*px4(i) -
94 . tempnc(nc7(i))*px1(i) - tempnc(nc8(i))*px2(i)
95
96 phiy = tempnc(nc1(i))*py1(i) + tempnc(nc2(i))*py2(i) +
97 . tempnc(nc3(i))*py3(i) + tempnc(nc4(i))*py4(i) -
98 . tempnc(nc5(i))*py3(i) - tempnc(nc6(i))*py4(i) -
99 . tempnc(nc7(i))*py1(i) - tempnc(nc8(i))*py2(i)
100
101 phiz = tempnc(nc1(i))*pz1(i) + tempnc(nc2(i))*pz2(i) +
102 . tempnc(nc3(i))*pz3(i) + tempnc(nc4(i))*pz4(i) -
103 . tempnc(nc5(i))*pz3(i) - tempnc(nc6(i))*pz4(i) -
104 . tempnc(nc7(i))*pz1(i) - tempnc(nc8(i))*pz2(i)
105
106 kc = (as + bs*tempel(i))*vol(i)*dt1*theaccfact
107 phix = kc*phix
108 phiy = kc*phiy
109 phiz = kc*phiz
110
111
112
113 a = one_over_8 * heat(i)
114 b = phix*px1(i) + phiy*py1(i) + pz1(i)*phiz
115 fphi(i,1) = a - b
116 fphi(i,7) = a + b
117 b = phix*px2(i) + phiy*py2(i) + pz2(i)*phiz
118 fphi(i,2) = a - b
119 fphi(i,8) = a + b
120 b = phix*px3(i) + phiy*py3(i) + pz3(i)*phiz
121 fphi(i,3) = a - b
122 fphi(i,5) = a + b
123 b = phix*px4(i) + phiy*py4(i) + pz4(i)*phiz
124 fphi(i,4) = a - b
125 fphi(i,6) = a + b
126 ENDDO
127
128 DO i=1,nel
129 IF (off(i)==zero.OR.offg(i)<=zero) cycle
130 m = iparts(i)
131 partsav(28,m) = partsav(28,m) + rhocp*vol0(i)*(tempel(i)-t0) + heat(i)
132 ENDDO
133
134 RETURN