40
41
42
44
45
46
47#include "implicit_f.inc"
48
49
50
51 INTEGER NSN, NMN, I0, NIR, I2SIZE, IDEL2,
52 . IRECT(4,*), MSR(*), NSV(*), IRTL(*), WEIGHT(*),
53 . IADI2(NIR,*),INDXC(*),IADX(*)
54
56 . a(3,*),ar(3,*), x(3,*),fskyi2(i2size,*),mmass(*),
57 . dpara(7,*), ms(*), in(*),stifn(*),stifr(*),dmast,adm(*),
58 . smass(*), siner(*),fsav(*), crst(2,*),fncont(3,*),
59 . fncontp(3,*) ,ftcontp(3,*)
60 TYPE (H3D_DATABASE) :: H3D_DATA
61
62
63
64#include "com01_c.inc"
65#include "scr14_c.inc"
66#include "scr16_c.inc"
67
68
69
70 INTEGER I, J, K,J1,J2,J3,J4, II, L, JJ, NN,NISKY2,I0BASE
71
73 . h(4),
74 . s,t,ss, st, xmsi, fs(3),sp,sm,tp,tm,
75 . mx,my,mz,det,fx0,fy0,fz0,ins,
76 . x0,x1,x2,x3,x4,xs,y0,y1,y2,y3,y4,ys,z0,z1,z2,z3,z4,zs,
77 . x12,x22,x32,x42,y12,y22,y32,y42,z12,z22,z32,z42,
78 . xx,yy,zz,xxx,yyy,zzz,xy,yz,zx,xy2,yz2,zx2,
79 . a1,a2,a3,b1,b2,b3,c1,c2,c3,mr,mrx,mry,mrz,inx,iny,inz,stf,fact,
80 . fx(4),fy(4),fz(4)
81
82
83
84
85
86
87 i0base = i0
88
89
90
91
92 IF(anim_n(2)+outp_n(2)+h3d_data%N_SCAL_DMAS >0) THEN
93 DO ii=1,nmn
94 j=msr(ii)
95 adm(j) = adm(j)*mmass(ii)
96 ENDDO
97 ENDIF
98
99
100
101
102
103
104
105
106
107
108
109 DO ii=1,nsn
110 k = indxc(ii)
111 IF (k == 0) cycle
112 i = nsv(k)
113 IF(i>0)THEN
114 l=irtl(ii)
115
116 s = crst(1,ii)
117 t = crst(2,ii)
118 sp=one + s
119 sm=one - s
120 tp=fourth*(one + t)
121 tm=fourth*(one - t)
122
123 h(1)=one/nir
124 h(2)=one/nir
125 h(3)=one/nir
126 h(4)=one/nir
127
128 j1=irect(1,l)
129 j2=irect(2,l)
130 j3=irect(3,l)
131 j4=irect(4,l)
132 x1=x(1,j1)
133 y1=x(2,j1)
134 z1=x(3,j1)
135 x2=x(1,j2)
136 y2=x(2,j2)
137 z2=x(3,j2)
138 x3=x(1,j3)
139 y3=x(2,j3)
140 z3=x(3,j3)
141 x4=x(1,j4)
142 y4=x(2,j4)
143 z4=x(3,j4)
144 x0=fourth*(x1+x2+x3+x4)
145 y0=fourth*(y1+y2+y3+y4)
146 z0=fourth*(z1+z2+z3+z4)
147 x1=x1-x0
148 y1=y1-y0
149 z1=z1-z0
150 x2=x2-x0
151 y2=y2-y0
152 z2=z2-z0
153 x3=x3-x0
154 y3=y3-y0
155 z3=z3-z0
156 x4=x4-x0
157 y4=y4-y0
158 z4=z4-z0
159 xs=x(1,i)-x0
160 ys=x(2,i)-y0
161 zs=x(3,i)-z0
162
163 x12=x1*x1
164 x22=x2*x2
165 x32=x3*x3
166 x42=x4*x4
167 y12=y1*y1
168 y22=y2*y2
169 y32=y3*y3
170 y42=y4*y4
171 z12=z1*z1
172 z22=z2*z2
173 z32=z3*z3
174 z42=z4*z4
175 xx=x12 + x22 + x32 + x42
176 yy=y12 + y22 + y32 + y42
177 zz=z12 + z22 + z32 + z42
178 xy=x1*y1 + x2*y2 + x3*y3 + x4*y4
179 yz=y1*z1 + y2*z2 + y3*z3 + y4*z4
180 zx=z1*x1 + z2*x2 + z3*x3 + z4*x4
181 zzz=xx+yy
182 xxx=yy+zz
183 yyy=zz+xx
184 xy2=xy*xy
185 yz2=yz*yz
186 zx2=zx*zx
187 det= xxx*yyy*zzz - xxx*yz2 - yyy*zx2 - zzz*xy2
188 . - two*xy*yz*zx
189 det=one/det
190 b1=zzz*yyy-yz2
191 b2=xxx*zzz-zx2
192 b3=yyy*xxx-xy2
193 c3=zzz*xy+yz*zx
194 c1=xxx*yz+zx*xy
195 c2=yyy*zx+xy*yz
196
197 dpara(1,ii)=det
198 dpara(2,ii)=b1
199 dpara(3,ii)=b2
200 dpara(4,ii)=b3
201 dpara(5,ii)=c1
202 dpara(6,ii)=c2
203 dpara(7,ii)=c3
204
205 IF (weight(i)==1) THEN
206 xmsi=ms(i)
207 fs(1)=a(1,i)
208 fs(2)=a(2,i)
209 fs(3)=a(3,i)
210 IF (iroddl==1) THEN
211 ins=in(i)
212 mx=ar(1,i) + ys*fs(3) - zs*fs(2)
213 my=ar(2,i) + zs*fs(1) - xs*fs(3)
214 mz=ar(3,i) + xs*fs(2) - ys*fs(1)
215 ELSE
216 ins=zero
217 mx=ys*fs(3) - zs*fs(2)
218 my=zs*fs(1) - xs*fs(3)
219 mz=xs*fs(2) - ys*fs(1)
220 ENDIF
221
222 a1=det*(mx*b1+my*c3+mz*c2)
223 a2=det*(my*b2+mz*c1+mx*c3)
224 a3=det*(mz*b3+mx*c2+my*c1)
225
226 fx0=fs(1)*fourth
227 fy0=fs(2)*fourth
228 fz0=fs(3)*fourth
229
230
231
232
233
234 inx=ins + xmsi*(xs*xs+ys*ys+zs*zs)
235 mrx = (b1+c3+c2)
236 mry = (b2+c1+c3)
237 mrz = (b3+c2+c1)
238 mr=det*inx*
max(mrx,mry,mrz)
239
240
241
242
243
244 fact = one
245 IF (iroddl==1) THEN
246 IF ((in(j1)>zero.AND.in(j2)>zero.AND.in(j3)>zero.AND.in(j4)>zero)) THEN
247
248 fact = zero
249 ENDIF
250 ENDIF
251
252 xmsi=fourth*xmsi+mr*fact
253
254 IF (iroddl == 1) THEN
255 stf = fourth*stifn(i)+ det*
max(mrx,mry,mrz)*(stifr(i)+stifn(i)*(xs*xs+ys*ys+zs*zs))
256 ELSE
257 stf = fourth*stifn(i)+ det*
max(mrx,mry,mrz)*(stifn(i)*(xs*xs+ys*ys+zs*zs))
258 ENDIF
259
260 i0 = i0base + iadx(k)
261 nn = iadi2(1,i0)
262 fx(1) = fx0 + a2*z1 - a3*y1
263 fy(1) = fy0 + a3*x1 - a1*z1
264 fz(1) = fz0 + a1*y1 - a2*x1
265 fskyi2(1,nn) = fx(1)
266 fskyi2(2,nn) = fy(1)
267 fskyi2(3,nn) = fz(1)
268 fskyi2(4,nn) = xmsi
269 fskyi2(5,nn) = stf
270 IF (iroddl == 1) THEN
271 fskyi2(6,nn) = zero
272 fskyi2(7,nn) = zero
273 fskyi2(8,nn) = zero
274 fskyi2(9,nn) = inx*fourth*(one-fact)
275 fskyi2(10,nn)= zero
276 ENDIF
277
278 nn = iadi2(2,i0)
279 fx(2) = fx0 + a2*z2 - a3*y2
280 fy(2) = fy0 + a3*x2 - a1*z2
281 fz(2) = fz0 + a1*y2 - a2*x2
282 fskyi2(1,nn) = fx(2)
283 fskyi2(2,nn) = fy(2)
284 fskyi2(3,nn) = fz(2)
285 fskyi2(4,nn) = xmsi
286 fskyi2(5,nn) = stf
287 IF (iroddl == 1) THEN
288 fskyi2(6,nn) = zero
289 fskyi2(7,nn) = zero
290 fskyi2(8,nn) = zero
291 fskyi2(9,nn) = inx*fourth*(one-fact)
292 fskyi2(10,nn)= zero
293 ENDIF
294
295 nn = iadi2(3,i0)
296 fx(3) = fx0 + a2*z3 - a3*y3
297 fy(3) = fy0 + a3*x3 - a1*z3
298 fz(3) = fz0 + a1*y3 - a2*x3
299 fskyi2(1,nn) = fx(3)
300 fskyi2(2,nn) = fy(3)
301 fskyi2(3,nn) = fz(3)
302 fskyi2(4,nn) = xmsi
303 fskyi2(5,nn) = stf
304 IF (iroddl == 1) THEN
305 fskyi2(6,nn) = zero
306 fskyi2(7,nn) = zero
307 fskyi2(8,nn) = zero
308 fskyi2(9,nn) = inx*fourth*(one-fact)
309 fskyi2(10,nn)= zero
310 ENDIF
311
312 nn = iadi2(4,i0)
313 fx(4) = fx0 + a2*z4 - a3*y4
314 fy(4) = fy0 + a3*x4 - a1*z4
315 fz(4) = fz0 + a1*y4 - a2*x4
316 fskyi2(1,nn) = fx(4)
317 fskyi2(2,nn) = fy(4)
318 fskyi2(3,nn) = fz(4)
319 fskyi2(4,nn) = xmsi
320 fskyi2(5,nn) = stf
321 IF (iroddl == 1) THEN
322 fskyi2(6,nn) = zero
323 fskyi2(7,nn) = zero
324 fskyi2(8,nn) = zero
325 fskyi2(9,nn) = inx*fourth*(one-fact)
326 fskyi2(10,nn)= zero
327 ENDIF
328
329 dmast = dmast + 4.*xmsi - ms(i)
330
331 IF (anim_n(2)+outp_n(2)+h3d_data%N_SCAL_DMAS >0) THEN
332 adm(j1) = adm(j1) + xmsi - fourth*ms(i)
333 adm(j2) = adm(j2) + xmsi - fourth*ms(i)
334 adm(j3) = adm(j3) + xmsi - fourth*ms(i)
335 adm(j4) = adm(j4) + xmsi - fourth*ms(i)
336 ENDIF
337 ENDIF
338
339 IF(idel2/=0.AND.ms(i)/=0.)smass(ii)=ms(i)
340 ms(i)=zero
341 stifn(i)=em20
342
343 IF (iroddl==1) THEN
344 IF(idel2/=0.AND.in(i)/=0.)siner(ii)=in(i)
345 in(i)=zero
346 stifr(i)=em20
347 ENDIF
348
349
351 . irect(1,l),nir ,fsav ,fncont ,fncontp,
352 . ftcontp ,weight ,h3d_data,i ,h)
353
354
355 ELSEIF(weight(-i)==1) THEN
356 i0 = i0base + iadx(k)
357 nn = iadi2(1,i0)
358 fskyi2(1,nn) = zero
359 fskyi2(2,nn) = zero
360 fskyi2(3,nn) = zero
361 fskyi2(4,nn) = zero
362 fskyi2(5,nn) = zero
363 IF (iroddl == 1) THEN
364 fskyi2(6,nn) = zero
365 fskyi2(7,nn) = zero
366 fskyi2(8,nn) = zero
367 fskyi2(9,nn) = zero
368 fskyi2(10,nn)= zero
369 ENDIF
370 nn = iadi2(2,i0)
371 fskyi2(1,nn) = zero
372 fskyi2(2,nn) = zero
373 fskyi2(3,nn) = zero
374 fskyi2(4,nn) = zero
375 fskyi2(5,nn) = zero
376 IF (iroddl == 1) THEN
377 fskyi2(6,nn) = zero
378 fskyi2(7,nn) = zero
379 fskyi2(8,nn) = zero
380 fskyi2(9,nn) = zero
381 fskyi2(10,nn)= zero
382 ENDIF
383 fskyi2(10,nn)= zero
384 nn = iadi2(3,i0)
385 fskyi2(1,nn) = zero
386 fskyi2(2,nn) = zero
387 fskyi2(3,nn) = zero
388 fskyi2(4,nn) = zero
389 fskyi2(5,nn) = zero
390 IF (iroddl == 1) THEN
391 fskyi2(6,nn) = zero
392 fskyi2(7,nn) = zero
393 fskyi2(8,nn) = zero
394 fskyi2(9,nn) = zero
395 fskyi2(10,nn)= zero
396 ENDIF
397 nn = iadi2(4,i0)
398 fskyi2(1,nn) = zero
399 fskyi2(2,nn) = zero
400 fskyi2(3,nn) = zero
401 fskyi2(4,nn) = zero
402 fskyi2(5,nn) = zero
403 IF (iroddl == 1) THEN
404 fskyi2(6,nn) = zero
405 fskyi2(7,nn) = zero
406 fskyi2(8,nn) = zero
407 fskyi2(9,nn) = zero
408 fskyi2(10,nn)= zero
409 ENDIF
410
411 ENDIF
412 ENDDO
413
414
415
416 IF(anim_n(2)+outp_n(2)+h3d_data%N_SCAL_DMAS >0) THEN
417#include "vectorize.inc"
418 DO ii=1,nmn
419 j=msr(ii)
420 adm(j) = adm(j)/
max(mmass(ii),em20)
421 ENDDO
422 ENDIF
423
424 RETURN
subroutine i2forces(x, fs, fx, fy, fz, irect, nir, fsav, fncont, fncontp, ftcontp, weight, h3d_data, nsl, h)