35
36
37
38
39
40
41#include "implicit_f.inc"
42
43
44
45#include "mvsiz_p.inc"
46
47
48
49#include "com08_c.inc"
50
51
52
53 INTEGER , INTENT(IN) :: NEL
54
55 my_real,
DIMENSION(MVSIZ) ,
INTENT(IN) ::
56 . x1, x2, x3, x4, x5, x6, x7, x8,
57 . y1, y2, y3, y4, y5, y6, y7, y8,
58 . z1, z2, z3, z4, z5, z6, z7, z8,
59 . vz1,vz2,vz3,vz4,vz5,vz6,vz7,vz8,
60 . rho,vol,ssp,nu,off
61 my_real,
DIMENSION(NEL) ,
INTENT(IN) :: vol0
62 my_real,
DIMENSION(NEL) ,
INTENT(INOUT) :: eint
63 my_real,
DIMENSION(MVSIZ) ,
INTENT(INOUT) ::
64 . f31,f32,f33,f35,f36,f37
65 my_real,
DIMENSION(NEL,2) ,
INTENT(INOUT) :: fhour
66
67
68
69 INTEGER :: I,J
70
72 . dett ,
73 . jaci1, jaci2, jaci3,
74 . jaci4, jaci5, jaci6,
75 . jaci7, jaci8, jaci9,
76 . x17(mvsiz) , x28(mvsiz) , x35(mvsiz) , x46(mvsiz
77 . y17(mvsiz) , y28(mvsiz) , y35(mvsiz) , y46(mvsiz),
78 . z17(mvsiz) , z28(mvsiz) , z35(mvsiz) , z46(mvsiz),
79 . jac_59_68(mvsiz), jac_67_49(mvsiz), jac_48_57(mvsiz),jac_19_37(mvsiz),
80 . jaci12, jaci45, jaci78,
81 . x_17_46 , x_28_35 ,
82 . y_17_46 , y_28_35 ,
83 . z_17_46 , z_28_35 ,
84 . hx,hy,hz,h1x,h1y,h1z,h2x,h2y,h2z,h3y,h3z,h4x,h4y,h4z,
85 . px1(mvsiz), px2(mvsiz), px3(mvsiz), px4(mvsiz),
86 . py1(mvsiz), py2(mvsiz), py3(mvsiz), py4(mvsiz),
87 . pz1(mvsiz), pz2(mvsiz), pz3(mvsiz), pz4(mvsiz),
88 . px1h1(mvsiz),px2h1(mvsiz),px3h1(mvsiz),px4h1(mvsiz),
89 . px1h2(mvsiz),px2h2(mvsiz),px3h2(mvsiz),px4h2(mvsiz),
90 . jac1(mvsiz),jac2(mvsiz),jac3
91 . jac4(mvsiz),jac5(mvsiz),jac6(mvsiz),
92 . jac7(mvsiz),jac8(mvsiz),jac9(mvsiz),det(mvsiz)
94 . g11(mvsiz),g21(mvsiz),g31(mvsiz),g41(mvsiz),
95 . g51(mvsiz),g61(mvsiz),g71(mvsiz),g81(mvsiz),
96 . g12(mvsiz),g22(mvsiz),g32(mvsiz),g42(mvsiz),
97 . g52(mvsiz),g62(mvsiz),g72(mvsiz),g82(mvsiz),
98 . fcl(mvsiz),nfhz1(mvsiz),nfhz2(mvsiz),vz17,vz28,vz35,vz46,
99 . h1vz,h2vz,hgz1(mvsiz),hgz2(mvsiz),cc,
100 . gg ,e_dt(mvsiz),thk,thk_1(mvsiz)
101
102 DO i=1,nel
103 x17(i)=x7(i)-x1(i)
104 x28(i)=x8(i)-x2(i)
105 x35(i)=x5(i)-x3(i)
106 x46(i)=x6(i)-x4(i)
107 y17(i)=y7(i)-y1(i)
108 y28(i)=y8(i)-y2(i)
109 y35(i)=y5(i)-y3(i)
110 y46(i)=y6(i)-y4(i)
111 z17(i)=z7(i)-z1(i)
112 z28(i)=z8(i)-z2(i)
113 z35(i)=z5(i)-z3(i)
114 z46(i)=z6(i)-z4(i)
115 END DO
116
117
118 DO i=1,nel
119
120 jac3(i)=x17(i)+x28(i)-x35(i)-x46(i)
121 jac1(i)=y17(i)+y28(i)-y35(i)-y46(i)
122 jac2(i)=z17(i)+z28(i)-z35(i)-z46(i)
123
124 x_17_46=x17(i)+x46(i)
125 x_28_35=x28(i)+x35(i)
126 y_17_46=y17(i)+y46(i)
127 y_28_35=y28(i)+y35(i)
128 z_17_46=z17(i)+z46(i)
129 z_28_35=z28(i)+z35(i)
130
131 jac6(i)=x_17_46+x_28_35
132 jac4(i)=y_17_46+y_28_35
133 jac5(i)=z_17_46+z_28_35
134
135 jac9(i)=x_17_46-x_28_35
136 jac7(i)=y_17_46-y_28_35
137 jac8(i)=z_17_46-z_28_35
138 jac_59_68(i)=jac5(i)*jac9(i)-jac6(i)*jac8(i)
139 jac_67_49(i)=jac6(i)*jac7(i)-jac4(i)*jac9(i)
140 jac_19_37(i)=jac1(i)*jac9(i)-jac3(i)*jac7(i)
141 jac_48_57(i)=jac4(i)*jac8(i)-jac5(i)*jac7(i)
142 det(i)=one_over_64*(jac1(i)*jac_59_68(i)+jac2(i)*jac_67_49(i)+jac3(i
143 thk =fourth*jac5(i)
144 thk_1(i) =one/thk
145 ENDDO
146
147 DO i=1,nel
148 dett=one_over_64/det(i)
149 jaci1=dett*jac_59_68(i)
150 jaci4=dett*jac_67_49(i)
151 jaci7=dett*jac_48_57(i)
152 jaci2=dett*(-jac2(i)*jac9(i)+jac3(i)*jac8(i))
153 jaci5=dett*jac_19_37(i)
154 jaci8=dett*(-jac1(i)*jac8(i)+jac2(i)*jac7(i))
155 jaci3=dett*( jac2(i)*jac6(i)-jac3(i)*jac5(i))
156 jaci6=dett*(-jac1(i)*jac6(i)+jac3(i)*jac4(i))
157 jaci9=dett*( jac1(i)*jac5(i)-jac2(i)*jac4(i))
158
159 jaci12=jaci1-jaci2
160 jaci45=jaci4-jaci5
161 jaci78=jaci7-jaci8
162
163 py3(i)= jaci12+jaci3
164 pz3(i)= jaci45+jaci6
165 px3(i)= jaci78+jaci9
166 py4(i)= jaci12-jaci3
167 pz4(i)= jaci45-jaci6
168 px4(i)= jaci78-jaci9
169
170 jaci12=jaci1+jaci2
171 jaci45=jaci4+jaci5
172 jaci78=jaci7+jaci8
173
174 py1(i)=-jaci12-jaci3
175 pz1(i)=-jaci45-jaci6
176 px1(i)=-jaci78-jaci9
177 py2(i)=-jaci12+jaci3
178 pz2(i)=-jaci45+jaci6
179 px2(i)=-jaci78+jaci9
180 ENDDO
181
182
183 DO i=1,nel
184 h1x=x1(i)+x2(i)-x3(i)-x4(i)-x5(i)-x6(i)+x7(i)+x8(i)
185 h1y=y1(i)+y2(i)-y3(i)-y4(i)-y5(i)-y6(i)+y7(i)+y8(i)
186 h1z=z1(i)+z2(i)-z3(i)-z4(i)-z5(i)-z6(i)+z7(i)+z8(i)
187 hx=one_over_8*h1x
188 hy=one_over_8*h1y
189 hz=one_over_8*h1z
190 px1h1(i)=px1(i)*hx+ py1(i)*hy+pz1(i)*hz
191 px2h1(i)=px2(i)*hx+ py2(i)*hy+pz2(i)*hz
192 px3h1(i)=px3(i)*hx+ py3(i)*hy+pz3(i)*hz
193 px4h1(i)=px4(i)*hx+ py4(i)*hy+pz4(i)*hz
194 END DO
195 DO i=1,nel
196 g11(i)= one_over_8-px1h1(i)
197 g21(i)= one_over_8-px2h1(i)
198 g31(i)=-one_over_8-px3h1(i)
199 g41(i)=-one_over_8-px4h1(i)
200 g51(i)=-one_over_8+px3h1(i)
201 g61(i)=-one_over_8+px4h1(i)
202 g71(i)= one_over_8+px1h1(i)
203 g81(i)= one_over_8+px2h1(i)
204 ENDDO
205
206
207 DO i=1,nel
208 h2x=x1(i)-x2(i)-x3(i)+x4(i)-x5(i)+x6(i)+x7(i)-x8(i)
209 h2y=y1(i)-y2(i)-y3(i)+y4(i)-y5(i)+y6(i)+y7(i)-y8(i)
210 h2z=z1(i)-z2(i)-z3(i)+z4(i)-z5(i)+z6(i)+z7(i)-z8(i)
211 hx=one_over_8*h2x
212 hy=one_over_8*h2y
213 hz=one_over_8*h2z
214 px1h2(i)=px1(i)*hx+ py1(i)*hy+pz1(i)*hz
215 px2h2(i)=px2(i)*hx+ py2(i)*hy+pz2(i)*hz
216 px3h2(i)=px3(i)*hx+ py3(i)*hy+pz3(i)*hz
217 px4h2(i)=px4(i)*hx+ py4(i)*hy+pz4(i)*hz
218 END DO
219 DO i=1,nel
220 g12(i)= one_over_8-px1h2(i)
221 g22(i)=-one_over_8-px2h2(i)
222 g32(i)=-one_over_8-px3h2(i)
223 g42(i)= one_over_8-px4h2(i)
224 g52(i)=-one_over_8+px3h2(i)
225 g62(i)= one_over_8+px4h2(i)
226 g72(i)= one_over_8+px1h2(i)
227 g82(i)=-one_over_8+px2h2(i)
228 ENDDO
229 DO i=1,nel
230 vz17=vz1(i)-vz7(i)
231 vz28=vz2(i)-vz8(i)
232 vz35=vz3(i)-vz5(i)
233 vz46=vz4(i)-vz6(i)
234 h1vz=vz1(i)+vz2(i)-vz3(i)-vz4(i)
235 . -vz5(i)-vz6(i)+vz7(i)+vz8(i)
236 hgz1(i)=one_over_8*h1vz-
237 . (px1h1(i)*vz17+px2h1(i)*vz28+
238 . px3h1(i)*vz35+px4h1(i)*vz46)
239 h2vz=vz1(i)-vz2(i)-vz3(i)+vz4(i)
240 . -vz5(i)+vz6(i)+vz7(i)-vz8(i)
241 hgz2(i)=one_over_8*h2vz-
242 . (px1h2(i)*vz17+px2h2(i)*vz28+
243 . px3h2(i)*vz35+px4h2(i)*vz46)
244 ENDDO
245
246 DO i=1,nel
247 gg =rho(i)*ssp(i)*ssp(i)*(one-two*nu(i))/(one-nu(i))
248 e_dt(i)=gg*off(i)*dt1*(one +nu(i))
249 ENDDO
250 DO i=1,nel
251 cc =zep1*e_dt(i)*thk_1(i)
252 fhour(i,1) = off(i)*fhour(i,1) + cc*hgz1(i)
253 fhour(i,2) = off(i)*fhour(i,2) + cc*hgz2(i)
254 ENDDO
255 DO i=1,nel
256 cc = thk_1(i)*vol(i)
257 fcl(i)=em02*rho(i)*off(i)*ssp(i)*vol(i)**two_third
258 nfhz1(i) = cc*fhour(i,1) + fcl(i)*hgz1(i)
259 nfhz2(i) = cc*fhour(i,2) + fcl(i)*hgz2(i)
260 ENDDO
261
262
263 DO i=1,nel
264 f31(i) =f31(i)-g11(i)*nfhz1(i)-g12(i)*nfhz2(i)
265 f32(i) =f32(i)-g21(i)*nfhz1(i)-g22(i)*nfhz2(i)
266 f33(i) =f33(i)-g31(i)*nfhz1(i)-g32(i)*nfhz2(i)
267 . -g41(i)*nfhz1(i)-g42(i)*nfhz2(i)
268 f35(i) =f35(i)-g51(i)*nfhz1(i)-g52(i)*nfhz2(i)
269 f36(i) =f36(i)-g61(i)*nfhz1(i)-g62(i)*nfhz2(i)
270 f37(i) =f37(i)-g71(i)*nfhz1(i)-g72(i)*nfhz2(i)
271 . -g81(i)*nfhz1(i)-g82(i)*nfhz2(i)
272 ENDDO
273 DO i=1,nel
274 eint(i)= eint(i)+dt1*(
275 . nfhz1(i)*hgz1(i) + nfhz2(i)*hgz2(i) )
277 ENDDO
278
279 RETURN
280