34
35
36
38
39
40
41
42
43
44#include "implicit_f.inc"
45
46
47
48#include "mvsiz_p.inc"
49#include "scr17_c.inc"
50#include "com08_c.inc"
51#include "impl1_c.inc"
52
53
54
55 INTEGER IXC(NIXC,*),JFT,JLT,NNOD,NPLAT,IPLAT(*),NPT,
56 . ISTACK(MVSIZ,NPT),INOD(*)
58 . offg(*),off(*)
59 parameter(nnod = 4)
61 . vcore(mvsiz,3,nnod),vxyz(mvsiz,3*nnod,npt),
62 . vq(mvsiz,3,3),zl(*),di(mvsiz,6),
63 . y24_t(*),x13_t(*),x24_t(*),y13_t(*),
area(*),
64 . vni(4,4), del_ply(mvsiz,12,npt), vr(3,*)
65
66
67
68 INTEGER J,I,K,EP,NN(4), IPLY,IP,NPLAT0
69 INTEGER L,M,J1,J2
71 . pg,z1,z2,mx13,my13,mx23,my23,mx34,my34,gama1,gama2, x21,
72 . x34,y21,y34 ,z21,z34,l12,l34,x41,x32,y41,y32,z41,z32,xx1
73 . yy,xy,xz1,yz ,zz,y24,x24,y13,x13,corel(3,4),xx1,yy,off_l,
74 . d1,d2,dt05,dt025,exz,eyz,ddry,v13x,v24x,vhix,ddrz1,ddrz2,
75 . ddrz,ddrx,a1,a2,d3
77 . vg13(3),vg24(3),vghi(3),v13(mvsiz,3), v24(mvsiz,3),
78 . vhi(mvsiz,3), ar(3),d(6),alr(3),rr(3,nnod),
79 . area_i(mvsiz), del_iply(mvsiz,3,npt-1),dn_iply(mvsiz,12,npt-1),
80 . dn_ply(mvsiz,12,npt)
81 DATA pg/.577350269189626/
82
83
84 nplat0 = nplat
85 a1 = one - pg
86 a2 = one + pg
87
88 vni(1,1)= fourth*a2*a2
89 vni(2,1)= fourth*a1*a2
90 vni(3,1)= fourth*a1*a1
91 vni(4,1)= vni(2,1)
92
93 vni(1,2)= vni(2,1)
94 vni(2,2)= vni(1,1)
95 vni(3,2)= vni(2,1)
96 vni(4,2)= vni(3,1)
97
98 vni(1,3)= vni(3,1)
99 vni(2,3)= vni(2,1)
100 vni(3,3)= vni(1,1)
101 vni(4,3)= vni(2,1)
102
103 vni(1,4)= vni(2,1)
104 vni(2,4)= vni(3,1)
105 vni(3,4)= vni(2,1)
106 vni(4,4)= vni(1,1)
107
108 DO ep=jft,jlt
109 area_i(ep)=one/
max(em20,
area(ep))
110 ENDDO
111
112 DO j=1, npt
113 DO ep=jft,jlt
114 ip = istack(ep,j)
115
116 nn(1) = inod(ixc(2,ep))
117 nn(2) = inod(ixc(3,ep))
118 nn(3) = inod(ixc(4,ep))
119 nn(4) = inod(ixc(5,ep))
120
121 vg13(1)=
ply(ip)%V(1,nn(1)) -
ply(ip)%V(1,nn(3))
122 vg24(1)=
ply(ip)%V(1,nn(2)) -
ply(ip
123 vghi(1)=
ply(ip)%V(1,nn(1)) -
ply(ip)%V(1,nn(2))
124 . +
ply(ip)%V(1,nn(3)) -
ply(ip)%V(1,nn(4))
125
126 vg13(2)=
ply(ip)%V(2,nn(1)) -
ply(ip)%V(2,nn(3))
127 vg24(2)=
ply(ip)%V(2,nn(2)) -
ply(ip)%V(2,nn(4))
128 vghi(2)=
ply(ip)%V(2,nn(1)) -
ply(ip)%V(2,nn(2))
129 . +
ply(ip)%V(2,nn(3)) -
ply(ip)%V(2,nn(4))
130
131 vg13(3)=
ply(ip)%V(3,nn(1)) -
ply(ip)%V(3,nn(3))
132 vg24(3)=
ply(ip)%V(3,nn(2)) -
ply(ip)%V(3,nn(4))
133 vghi(3)=
ply(ip)%V(3,nn(1)) -
ply(ip)%V(3,nn(2))
134 . +
ply(ip)%V(3,nn(3)) -
ply(ip)%V(3,nn(4))
135
136 v13(ep,1) =(vq(ep,1,1)*vg13(1)+vq(ep,2,1)*vg13(2)
137 1 +vq(ep,3,1)*vg13(3))
138 v24(ep,1)=(vq(ep,1,1)*vg24(1)+vq(ep,2,1)*vg24(2)
139 1 +vq(ep,3,1)*vg24(3))
140 vhi(ep,1)=(vq(ep,1,1)*vghi(1)+vq(ep,2,1)*vghi(2)
141 1 +vq(ep,3,1)*vghi(3))
142 v13(ep,2)=(vq(ep,1,2)*vg13(1)+vq(ep,2,2)*vg13(2)
143 1 +vq(ep,3,2)*vg13(3))
144 v24(ep,2)=(vq(ep,1,2)*vg24(1)+vq(ep,2,2)*vg24(2)
145 1 +vq(ep,3,2)*vg24(3))
146 vhi(ep,2)=(vq(ep,1,2)*vghi(1)+vq(ep,2,2)*vghi(2)
147 1 +vq(ep,3,2)*vghi(3))
148 v13(ep,3)=(vq(ep,1,3)*vg13(1)+vq(ep,2,3)*vg13(2)
149 1 +vq(ep,3,3)*vg13(3))
150 v24(ep,3)=(vq(ep,1,3)*vg24(1)+vq(ep,2,3)*vg24(2)
151 1 +vq(ep,3,3)*vg24(3))
152 vhi(ep,3)=(vq(ep,1,3)*vghi(1)+vq(ep,2,3)*vghi(2)
153 1 +vq(ep,3,3)*vghi(3))
154
155 ENDDO
156
157 IF (impl_s==0) THEN
158 dt05 = half*dt1
159 dt025 =fourth*dt1
160 DO i=jft,jlt
161 exz = y24_t(i)*v13(i,3)-y13_t(i)*v24(i,3)
162 eyz = -x24_t(i)*v13(i,3)+x13_t(i)*v24(i,3)
163 ddry=dt05*exz*area_i(i)
164 ddrx=dt05*eyz*area_i(i)
165 v13x = v13(i,1)
166 v24x = v24(i,1)
167 vhix = vhi(i,1)
168 ddrz1=dt025*(v13(i,2)-v24(i,2))/(x13_t(i)-x24_t(i))
169 IF (abs(x13_t(i)-x24_t(i))<em10) ddrz1 = zero
170 v13(i,1) = v13(i,1)-ddry*v13(i,3)-ddrz1*v13
171 v24(i,1) = v24(i,1)-ddry*v24(i,3)-ddrz1*v24(i,2)
172 vhi(i,1) = vhi(i,1)-ddry*vhi(i,3)-ddrz1*vhi(i,2)
173 ddrz2=dt025*(v13x+v24x)/(y13_t(i)+y24_t(i))
174 IF (abs(y13_t(i)+y24_t(i))<em10) ddrz2 = zero
175 v13(i,2) = v13(i,2)-ddrx*v13(i,3)-ddrz2*v13x
176 v24(i,2) = v24(i,2)-ddrx*v24(i,3)-ddrz2*v24x
177 vhi(i,2) = vhi(i,2)-ddrx*vhi(i,3)-ddrz2*vhix
178 ENDDO
179 ENDIF
180
181#include "vectorize.inc"
182 DO i=jft,nplat0
183 ep =iplat(i)
184 vxyz(ep,1,j)=v13(ep,1)
185 vxyz(ep,2,j)=v13(ep,2)
186 vxyz(ep,3,j)=v13(ep,3)
187
188 vxyz(ep,4,j)=v24(ep,1)
189 vxyz(ep,5,j)=v24(ep,2)
190 vxyz(ep,6,j)=v24(ep,3)
191
192 vxyz(ep,7,j)=vhi(ep,1)
193 vxyz(ep,8,j)=vhi(ep,2)
194 vxyz(ep,9,j)=vhi(ep,3)
195
196 ENDDO
197#include "vectorize.inc"
198 DO i=nplat0+1,jlt
199 ep =iplat(i)
200 z1=zl(ep)
201 z2=z1*z1
202
203
204 x13 =(vcore(ep,1,1)-vcore(ep,1,3))*half
205 x24 =(vcore(ep,1,2)-vcore(ep,1,4))*half
206 y13 =(vcore(ep,2,1)-vcore(ep,2,3))*half
207 y24 =(vcore(ep,2,2)-vcore(ep,2,4))*half
208 mx13=(vcore(ep,1,1)+vcore(ep,1,3))*half
209 my13=(vcore(ep,2,1)+vcore(ep,2,3))*half
210
211
212 ar(1)=-z1*vhi(ep,2)+y13*v13(ep,3)+y24*v24(ep,3)+my13*vhi(ep,3)
213 ar(2)= z1*vhi(ep,1)-x13*v13(ep,3)-x24*v24(ep,3)-mx13*vhi(ep,3)
214 ar(3)= x13*v13(ep,2)+x24*v24(ep,2)+mx13*vhi(ep,2)
215 1 -y13*v13(ep,1)-y24*v24(ep,1)-my13*vhi(ep,1)
216
217 alr(1) =di(ep,1)*ar(1)+di(ep,4)*ar(2)+di(ep,5)*ar(3)
218 alr(2) =di(ep,4)*ar(1)+di(ep,2)*ar(2)+di(ep,6)*ar(3)
219 alr(3) =di(ep,5)*ar(1)+di(ep,6)*ar(2)+di(ep,3)*ar(3)
220
221 v13(ep,1)= half*v13(ep,1)+alr(3)*y13
222 v24(ep,1)= half*v24(ep,1)+alr(3)*y24
223 vhi(ep,1)= fourth*vhi(ep,1)+(alr(3)*my13-z1*alr(2))
224 v13(ep,2)= half*v13(ep,2)-alr(3)*x13
225 v24(ep,2)= half*v24(ep,2)-alr(3)*x24
226 vhi(ep,2)= fourth*vhi(ep,2)-(alr(3)*mx13-z1*alr(1))
227 v13(ep,3)= half*v13(ep,3)-(y13*alr(1)-x13*alr(2))
228 v24(ep,3)= half*v24(ep,3)-(y24*alr(1)-x24*alr(2))
229 vhi(ep,3)= fourth*vhi(ep,3)+(mx13*alr(2)-my13*alr(1))
230
231 vxyz(ep,1 ,j) = v13(ep,1) + vhi(ep,1)
232 vxyz(ep,4 ,j) = v24(ep,1) - vhi(ep,1)
233 vxyz(ep,7 ,j) = -v13(ep,1) + vhi(ep,1)
234 vxyz(ep,10,j) = -v24(ep,1) - vhi(ep,1)
235
236 vxyz(ep,2 ,j) = v13(ep,2) + vhi(ep,2)
237 vxyz(ep,5 ,j) = v24(ep,2) - vhi(ep,2)
238 vxyz(ep,8 ,j) = -v13(ep,2) + vhi(ep,2)
239 vxyz(ep,11,j) = -v24(ep,2) - vhi(ep,2)
240
241 vxyz(ep,3 ,j) = v13(ep,3) + vhi(ep,3)
242 vxyz(ep,6 ,j) = v24(ep,3) - vhi(ep,3)
243 vxyz(ep,9 ,j) = -v13(ep,3) + vhi(ep,3)
244 vxyz(ep,12,j) = -v24(ep,3) - vhi(ep,3)
245 ENDDO
246 ENDDO
247
248
249
250 off_l = 0.
251 DO ep=jft,jlt
252 off(ep) =
min(one,abs(offg(ep)))
253 off_l =
min(off_l,offg(ep))
254 ENDDO
255 IF(off_l < zero)THEN
256 DO j=1,npt
257 DO ep=jft,jlt
258 IF(offg(ep) < zero)THEN
259 vxyz(ep,1 ,j)= zero
260 vxyz(ep,2 ,j)= zero
261 vxyz(ep,3 ,j)= zero
262 vxyz(ep,4 ,j)= zero
263 vxyz(ep,5 ,j)= zero
264 vxyz(ep,6 ,j)= zero
265 vxyz(ep,7 ,j)= zero
266 vxyz(ep,8 ,j)= zero
267 vxyz(ep,9 ,j)= zero
268 vxyz(ep,10,j)= zero
269 vxyz(ep,11,j)= zero
270 vxyz(ep,12,j)= zero
271 ENDIF
272 ENDDO
273 ENDDO
274 ENDIF
275
276
277 DO j=1, npt
278 DO ep=jft,jlt
279 ip = istack(ep,j)
280 nn(1) = inod(ixc(2,ep))
281 nn(2) = inod(ixc(3,ep))
282 nn(3) = inod(ixc(4,ep))
283 nn(4) = inod(ixc(5,ep))
284
285 d1 =
ply(ip)%U(1,nn(1))
286 d2 =
ply(ip)%U(2,nn(1))
287 d3 =
ply(ip)%U(3,nn(1))
288 dn_ply(ep,1,j)=vq(ep,1,1)*d1 + vq(ep,2,1)*d2 + vq(ep,3,1)*d3
289 dn_ply(ep,2,j)=vq(ep,1,2)*d1 + vq(ep,2,2)*d2 + vq(ep,3,2)*d3
290 dn_ply(ep,3,j)=vq(ep,1,3)*d1 + vq(ep,2,3)*d2 + vq(ep,3,3)*d3
291
292 d1 =
ply(ip)%U(1,nn(2))
293 d2 =
ply(ip)%U(2,nn(2))
294 d3 =
ply(ip)%U(3,nn(2))
295 dn_ply(ep,4,j)=vq(ep,1,1)*d1 + vq(ep,2,1)*d2 + vq(ep,3,1)*d3
296 dn_ply(ep,5,j)=vq(ep,1,2)*d1 + vq(ep,2,2)*d2 + vq(ep,3,2)*d3
297 dn_ply(ep,6,j)=vq(ep,1,3)*d1 + vq(ep,2,3)*d2 + vq(ep,3,3)*d3
298
299 d1 =
ply(ip)%U(1,nn(3))
300 d2 =
ply(ip)%U(2,nn(3))
301 d3 =
ply(ip)%U(3,nn(3))
302 dn_ply(ep,7,j)=vq(ep,1,1)*d1 + vq(ep,2,1)*d2 + vq(ep,3,1)*d3
303 dn_ply(ep,8,j)=vq(ep,1,2)*d1 + vq(ep,2,2)*d2 + vq(ep,3,2)*d3
304 dn_ply(ep,9,j)=vq(ep,1,3)*d1 + vq(ep,2,3)*d2 + vq(ep,3,3)*d3
305
306 d1 =
ply(ip)%U(1,nn(4))
307 d2 =
ply(ip)%U(2,nn(4))
308 d3 =
ply(ip)%U(3,nn(4))
309 dn_ply(ep,10,j)=vq(ep,1,1)*d1 + vq(ep,2,1)*d2 + vq(ep,3,1)*d3
310 dn_ply(ep,11,j)=vq(ep,1,2)*d1 + vq(ep,2,2)*d2 + vq(ep,3,2)*d3
311 dn_ply(ep,12,j)=vq(ep,1,3)*d1 + vq(ep,2,3)*d2 + vq(ep,3,3)*d3
312 ENDDO
313 ENDDO
314
315 DO j=1, npt
316 DO ep=jft,jlt
317 del_ply(ep,1,j) = (
318 . dn_ply(ep,1,j)*vni(1,1) + dn_ply(ep,4,j)*vni(2,1) +
319 . dn_ply(ep,7,j)*vni(3,1) + dn_ply(ep,10,j)*vni(4,1) )
320 del_ply(ep,2,j) = (
321 . dn_ply(ep,2,j)*vni(1,1) + dn_ply(ep,5,j)*vni(2,1) +
322 . dn_ply(ep,8,j)*vni(3,1) + dn_ply(ep,11,j)*vni(4,1) )
323 del_ply(ep,3,j) = (
324 . dn_ply(ep,3,j)*vni(1,1) + dn_ply(ep,6,j)*vni(2,1) +
325 . dn_ply(ep,9,j)*vni(3,1) + dn_ply(ep,12,j)*vni(4,1) )
326
327 del_ply(ep,4,j) = (
328 . dn_ply(ep,1,j)*vni(1,2) + dn_ply(ep,4,j)*vni(2,2) +
329 . dn_ply(ep,7,j)*vni(3,2) + dn_ply(ep,10,j)*vni(4,2) )
330 del_ply(ep,5,j) = (
331 . dn_ply(ep,2,j)*vni(1,2) + dn_ply(ep,5,j)*vni(2,2) +
332 . dn_ply(ep,8,j)*vni(3,2) + dn_ply(ep,11,j)*vni(4,2) )
333 del_ply(ep,6,j) = (
334 . dn_ply(ep,3,j)*vni(1,2) + dn_ply(ep,6,j)*vni(2,2) +
335 . dn_ply(ep,9,j)*vni(3,2) + dn_ply(ep,12,j)*vni(4,2) )
336
337 del_ply(ep,7,j) = (
338 . dn_ply(ep,1,j)*vni(1,3) + dn_ply(ep,4,j)*vni(2,3) +
339 . dn_ply(ep,7,j)*vni(3,3) + dn_ply(ep,10,j)*vni(4,3) )
340 del_ply(ep,8,j) = (
341 . dn_ply(ep,2,j)*vni(1,3) + dn_ply(ep,5,j)*vni(2,3) +
342 . dn_ply(ep,8,j)*vni(3,3) + dn_ply(ep,11,j)*vni(4,3) )
343 del_ply(ep,9,j) = (
344 . dn_ply(ep,3,j)*vni(1,3) + dn_ply(ep,6,j)*vni(2,3) +
345 . dn_ply(ep,9,j)*vni(3,3) + dn_ply(ep,12,j)*vni(4,3) )
346
347 del_ply(ep,10,j) = (
348 . dn_ply(ep,1,j)*vni(1,4) + dn_ply(ep,4,j)*vni(2,4) +
349 . dn_ply(ep,7,j)*vni(3,4) + dn_ply(ep,10,j)*vni(4,4) )
350 del_ply(ep,11,j) = (
351 . dn_ply(ep,2,j)*vni(1,4) + dn_ply(ep,5,j)*vni(2,4) +
352 . dn_ply(ep,8,j)*vni(3,4) + dn_ply(ep,11,j)*vni(4,4) )
353 del_ply(ep,12,j) = (
354 . dn_ply(ep,3,j)*vni(1,4) + dn_ply(ep,6,j)*vni(2,4) +
355 . dn_ply(ep,9,j)*vni(3,4) + dn_ply(ep,12,j)*vni(4,4) )
356 ENDDO
357 ENDDO
358
359 RETURN
subroutine area(d1, x, x2, y, y2, eint, stif0)
type(ply_data), dimension(:), allocatable ply