41
42
43
45
46
47
48#include "implicit_f.inc"
49
50
51
52 INTEGER NSN, NMN, IDEL2,
53 . IRECT(4,*), MSR(*), NSV(*), IRTL(*), WEIGHT(*),INDXC(NSN),MSEGTYP2(*)
54
56 . x(3,*),v(3,*),a(3,*),ms(*),crst(2,*),stifn(*),mmass(*),smass(*),
57 . fsav(*),fncont(3,*),in(*),siner(*),dpara(7,*),ar(3,*),stifr(*),csts_bis(2,*),
58 . t2fac_sms(*),fncontp(3,*) ,ftcontp(3,*)
59 TYPE (H3D_DATABASE) :: H3D_DATA
60
61
62
63#include "com01_c.inc"
64#include "impl1_c.inc"
65#include "sms_c.inc"
66
67
68
69 INTEGER NIR, I, J, K, II, L, JJ,
70 . IX1,IX2,IX3,IX4
71
73 . h(4),xmsj, ss, st, xmsi,fs(3),sp,sm,tp,tm,
74 . e1x,e1y,e1z,e2x,e2y,e2z,e3x,e3y,e3z,
75 . x1,x2,x3,x4,y1,y2,y3,y4,z1,z2,z3,z4,x0,y0,z0,xs(3),xm(3),
76 . stifm,fmx(4),fmy(4),fmz(4),fx(4),fy(4),fz(4),
77 . rx(4),ry(4),rz(4),rs(3),flx,fly,flz,fac_triang,stbrk,
78 . mxs,mys,mzs,stifmr,dwdu,rm(3),betax,betay,h2(4),mxi,myi,mzi
79
80 nir=2
81 IF(n2d==0)nir=4
82
83 IF(impl_s>0) THEN
84 DO ii=1,nsn
85 k = indxc(ii)
86 IF (k == 0) cycle
87 i = nsv(k)
88 IF(i>0)THEN
89 l=irtl(ii)
90
91 xmsi=ms(i)*weight(i)
92 fs(1)=a(1,i)*weight(i)
93 fs(2)=a(2,i)*weight(i)
94 fs(3)=a(3,i)*weight(i)
95
96 IF (iroddl == 1) THEN
97 mxi=ar(1,i)*weight(i)
98 myi=ar(2,i)*weight(i)
99 mzi=ar(3,i)*weight(i)
100 ENDIF
101
102 ix1 = irect(1,l)
103 ix2 = irect(2,l)
104 ix3 = irect(3,l)
105 ix4 = irect(4,l)
106
107 IF (ix3 == ix4) THEN
108
109 nir = 3
110 h(1) = crst(1,ii)
111 h(2) = crst(2,ii)
112 h(3) = one-crst(1,ii)-crst(2,ii)
113 h(4) = zero
114 h2(1) = csts_bis(1,ii)
115 h2(2) = csts_bis(2,ii)
116 h2(3) = one-csts_bis(1,ii)-csts_bis(2,ii)
117 h2(4) = zero
118 ELSE
119
120 nir = 4
121 ss=crst(1,ii)
122 st=crst(2,ii)
123 sp=one + ss
124 sm=one - ss
125 tp=fourth*(one + st)
126 tm=fourth*(one - st)
127 h(1)=tm*sm
128 h(2)=tm*sp
129 h(3)=tp*sp
130 h(4)=tp*sm
131
132
133 ss=csts_bis(1,ii)
134 st=csts_bis(2,ii)
135 sp=one + ss
136 sm=one - ss
137 tp=fourth*(one + st)
138 tm=fourth*(one - st)
139 h2(1)=tm*sm
140 h2(2)=tm*sp
141 h2(3)=tp*sp
142 h2(4)=tp*sm
143 ENDIF
144
145 IF (msegtyp2(l)==0) THEN
146
147
148
149
150
151
152
153 x1 = x(1,ix1)
154 y1 = x(2,ix1)
155 z1 = x(3,ix1)
156 x2 = x(1,ix2)
157 y2 = x(2,ix2)
158 z2 = x(3,ix2)
159 x3 = x(1,ix3)
160 y3 = x(2,ix3)
161 z3 = x(3,ix3)
162 x4 = x(1,ix4)
163 y4 = x(2,ix4)
164 z4 = x(3,ix4)
165 xs(1) = x(1,i)
166 xs(2) = x(2,i)
167 xs(3) = x(3,i)
168
169 CALL i2rep(x1 ,x2 ,x3 ,x4 ,
170 . y1 ,y2 ,y3 ,y4 ,
171 . z1 ,z2 ,z3 ,z4 ,
172 . e1x ,e1y ,e1z ,
173 . e2x ,e2y ,e2z ,
174 . e3x ,e3y ,e3z ,nir)
175
176
177 IF (nir == 4) THEN
178 fac_triang = one
179 x0 = fourth*(x1 + x2 + x3 + x4)
180 y0 = fourth*(y1 + y2 + y3 + y4)
181 z0 = fourth*(z1 + z2 + z3 + z4)
182 ELSE
183 fac_triang = zero
184 x0 = third*(x1 + x2 + x3)
185 y0 = third*(y1 + y2 + y3)
186 z0 = third*(z1 + z2 + z3)
187 ENDIF
188
189 xs(1) = xs(1) - x0
190 xs(2) = xs(2) - y0
191 xs(3) = xs(3) - z0
192
193 x1 = x1 - x0
194 y1 = y1 - y0
195 z1 = z1 - z0
196 x2 = x2 - x0
197 y2 = y2 - y0
198 z2 = z2 - z0
199 x3 = x3 - x0
200 y3 = y3 - y0
201 z3 = z3 - z0
202 x4 = x4 - x0
203 y4 = y4 - y0
204 z4 = z4 - z0
205 IF (nir==3) THEN
206 x4 = zero
207 y4 = zero
208 z4 = zero
209 END IF
210
211 xm(1) = x1*h(1) + x2*h(2) + x3*h(3) + x4*h(4)
212 xm(2) = y1*h(1) + y2*h(2) + y3*h(3) + y4*h(4)
213 xm(3) = z1*h(1) + z2*h(2) + z3*h(3) + z4*h(4)
214
215
216
217 rs(1) = xs(1)*e1x + xs(2)*e1y + xs(3)*e1z
218 rs(2) = xs(1)*e2x + xs(2)*e2y + xs(3)*e2z
219 rs(3) = xs(1)*e3x + xs(2)*e3y + xs(3)*e3z
220
221 rm(1) = xm(1)*e1x + xm(2)*e1y + xm(3)*e1z
222 rm(2) = xm(1)*e2x + xm(2)*e2y + xm(3)*e2z
223 rm(3) = xm(1)*e3x + xm(2)*e3y + xm(3)*e3z
224
225 rx(1) = e1x*x1 + e1y*y1 + e1z*z1
226 ry(1) = e2x*x1 + e2y*y1 + e2z*z1
227 rz(1) = e3x*x1 + e3y*y1 + e3z*z1
228 rx(2) = e1x*x2 + e1y*y2 + e1z*z2
229 ry(2) = e2x*x2 + e2y*y2 + e2z*z2
230 rz(2) = e3x*x2 + e3y*y2 + e3z*z2
231 rx(3) = e1x*x3 + e1y*y3 + e1z*z3
232 ry(3) = e2x*x3 + e2y*y3 + e2z*z3
233 rz(3) = e3x*x3 + e3y*y3 + e3z*z3
234 rx(4) = e1x*x4 + e1y*y4 + e1z*z4
235 ry(4) = e2x*x4 + e2y*y4 + e2z*z4
236 rz(4) = e3x*x4 + e3y*y4 + e3z*z4
237
238
239 CALL i2cin_rot27(stbrk,rs,rm,rx(1),ry(1),rz(1),rx(2),ry(2),rz(2),rx(3),ry(3),rz(3),
240 . rx(4),ry(4),rz(4),dpara(1,ii),dwdu,e1x,e1y,e1z,e2x,e2y,e2z,e3x,e3y,e3z,
241 . nir,betax,betay)
242
243
244
245 flx = fs(1)*e1x + fs(2)*e1y + fs(3)*e1z
246 fly = fs(1)*e2x + fs(2)*e2y + fs(3)*e2z
247 flz = fs(1)*e3x + fs(2)*e3y + fs(3)*e3z
248
249 DO j=1,4
250 fmx(j) = h(j)*flx
251 fmy(j) = h(j)*fly
252 fmz(j) = h(j)*flz
253 ENDDO
254
255
256 IF (iroddl==1) THEN
257 mxs = mxi*e1x + myi*e1y + mzi*e1z
258 mys = mxi*e2x + myi*e2y + mzi*e2z
259 mzs = mxi*e3x + myi*e3y + mzi*e3z
260
261
263 . fmx ,fmy ,fmz ,h ,stifm ,
264 . mxs ,mys ,mzs ,stifmr ,betax ,
265 . betay)
266 ELSE
267 mxs = zero
268 mys = zero
269 mzs = zero
270
271
273 . fmx ,fmy ,fmz ,h ,stifm ,
274 . mxs ,mys ,mzs ,stifmr ,betax ,
275 . betay)
276 stifmr = zero
277 ENDIF
278
279
280
281 DO j=1,4
282 fx(j) = e1x*fmx(j) + e2x*fmy(j) + e3x*fmz(j)
283 fy(j) = e1y*fmx(j) + e2y*fmy(j) + e3y*fmz(j)
284 fz(j) = e1z*fmx(j) + e2z*fmy(j) + e3z*fmz(j)
285 ENDDO
286
287 ELSE
288
289
290
291
292 stifm=zero
293 stbrk=zero
294 stifmr=zero
295 dwdu=zero
296
297 DO j=1,4
298 fx(j) = fs(1)*h(j)
299 fy(j) = fs(2)*h(j)
300 fz(j) = fs(3)*h(j)
301 ENDDO
302
303 ENDIF
304
305 IF(iroddl/=0)THEN
306 DO jj=1,nir
307 j=irect(jj,l)
308 a(1,j)=a(1,j)+fx(jj)
309 a(2,j)=a(2,j)+fy(jj)
310 a(3,j)=a(3,j)+fz(jj)
311 ms(j)=ms(j)+xmsi*h2(jj)
312 stifn(j)=stifn(j)+weight(i)*(stifn(i)*(one+stbrk)*(abs(h(jj))+stifm)+stifr(i)*stifmr*dwdu)
313 ENDDO
314 ELSE
315 DO jj=1,nir
316 j=irect(jj,l)
317 a(1,j)=a(1,j)+fx(jj)
318 a(2,j)=a(2,j)+fy(jj)
319 a(3,j)=a(3,j)+fz(jj)
320 ms(j)=ms(j)+xmsi*h2(jj)
321 stifn(j)=stifn(j)+weight(i)*(stifn(i)*(one+stbrk)*(abs(h(jj))+stifm))
322 ENDDO
323 END IF
324
325
327 . irect(1,l),nir ,fsav ,fncont ,fncontp,
328 . ftcontp ,weight ,h3d_data,i ,h)
329
330 IF(iroddl==0)THEN
331 IF(idel2/=0.AND.ms(i)/=0.)smass(ii)=ms(i)
332 ms(i)=zero
333 stifn(i)=em20
334 a(1,i)=zero
335 a(2,i)=zero
336 a(3,i)=zero
337 ENDIF
338
339 ENDIF
340
341 ENDDO
342
343
344 ELSE
345
346 DO ii=1,nsn
347 k = indxc(ii)
348 IF (k == 0) cycle
349 i = nsv(k)
350 IF(i>0)THEN
351 l=irtl(ii)
352
353 xmsi=ms(i)*weight(i)
354 fs(1)=a(1,i)*weight(i)
355 fs(2)=a(2,i)*weight(i)
356 fs(3)=a(3,i)*weight(i)
357
358 IF (iroddl == 1) THEN
359 mxi=ar(1,i)*weight(i)
360 myi=ar(2,i)*weight(i)
361 mzi=ar(3,i)*weight(i)
362 ENDIF
363
364 ix1 = irect(1,l)
365 ix2 = irect(2,l)
366 ix3 = irect(3,l)
367 ix4 = irect(4,l)
368
369 IF (ix3 == ix4) THEN
370
371 nir = 3
372 h(1) = crst(1,ii)
373 h(2) = crst(2,ii)
374 h(3) = one-crst(1,ii)-crst(2,ii)
375 h(4) = zero
376 h2(1) = csts_bis(1,ii)
377 h2(2) = csts_bis(2,ii)
378 h2(3) = one-csts_bis(1,ii)-csts_bis(2,ii)
379 h2(4) = zero
380 ELSE
381
382 nir = 4
383 ss=crst(1,ii)
384 st=crst(2,ii)
385 sp=one + ss
386 sm=one - ss
387 tp=fourth*(one + st)
388 tm=fourth*(one - st)
389 h(1)=tm*sm
390 h(2)=tm*sp
391 h(3)=tp*sp
392 h(4)=tp*sm
393
394
395 ss=csts_bis(1,ii)
396 st=csts_bis(2,ii)
397 sp=one + ss
398 sm=one - ss
399 tp=fourth*(one + st)
400 tm=fourth*(one - st)
401 h2(1)=tm*sm
402 h2(2)=tm*sp
403 h2(3)=tp*sp
404 h2(4)=tp*sm
405 ENDIF
406
407 IF (msegtyp2(l)==0) THEN
408
409
410
411
412
413
414
415 x1 = x(1,ix1)
416 y1 = x(2,ix1)
417 z1 = x(3,ix1)
418 x2 = x(1,ix2)
419 y2 = x(2,ix2)
420 z2 = x(3,ix2)
421 x3 = x(1,ix3)
422 y3 = x(2,ix3)
423 z3 = x(3,ix3)
424 x4 = x(1,ix4)
425 y4 = x(2,ix4)
426 z4 = x(3,ix4)
427 xs(1) = x(1,i)
428 xs(2) = x(2,i)
429 xs(3) = x(3,i)
430
431 CALL i2rep(x1 ,x2 ,x3 ,x4 ,
432 . y1 ,y2 ,y3 ,y4 ,
433 . z1 ,z2 ,z3 ,z4 ,
434 . e1x ,e1y ,e1z ,
435 . e2x ,e2y ,e2z ,
436 . e3x ,e3y ,e3z ,nir)
437
438
439 IF (nir == 4) THEN
440 fac_triang = one
441 x0 = fourth*(x1 + x2 + x3 + x4)
442 y0 = fourth*(y1 + y2 + y3 + y4)
443 z0 = fourth*(z1 + z2 + z3 + z4)
444 ELSE
445 fac_triang = zero
446 x0 = third*(x1 + x2 + x3)
447 y0 = third*(y1 + y2 + y3)
448 z0 = third*(z1 + z2 + z3)
449 ENDIF
450
451 xs(1) = xs(1) - x0
452 xs(2) = xs(2) - y0
453 xs(3) = xs(3) - z0
454
455 x1 = x1 - x0
456 y1 = y1 - y0
457 z1 = z1 - z0
458 x2 = x2 - x0
459 y2 = y2 - y0
460 z2 = z2 - z0
461 x3 = x3 - x0
462 y3 = y3 - y0
463 z3 = z3 - z0
464 x4 = x4 - x0
465 y4 = y4 - y0
466 z4 = z4 - z0
467 IF (nir==3) THEN
468 x4 = zero
469 y4 = zero
470 z4 = zero
471 END IF
472
473 xm(1) = x1*h(1) + x2*h(2) + x3*h(3) + x4*h(4)
474 xm(2) = y1*h(1) + y2*h(2) + y3*h(3) + y4*h(4)
475 xm(3) = z1*h(1) + z2*h(2) + z3*h(3) + z4*h(4)
476
477
478
479 rs(1) = xs(1)*e1x + xs(2)*e1y + xs(3)*e1z
480 rs(2) = xs(1)*e2x + xs(2)*e2y + xs(3)*e2z
481 rs(3) = xs(1)*e3x + xs(2)*e3y + xs(3)*e3z
482
483 rm(1) = xm(1)*e1x + xm(2)*e1y + xm(3)*e1z
484 rm(2) = xm(1)*e2x + xm(2)*e2y + xm(3)*e2z
485 rm(3) = xm(1)*e3x + xm(2)*e3y + xm(3)*e3z
486
487 rx(1) = e1x*x1 + e1y*y1 + e1z*z1
488 ry(1) = e2x*x1 + e2y*y1 + e2z*z1
489 rz(1) = e3x*x1 + e3y*y1 + e3z*z1
490 rx(2) = e1x*x2 + e1y*y2 + e1z*z2
491 ry(2) = e2x*x2 + e2y*y2 + e2z*z2
492 rz(2) = e3x*x2 + e3y*y2 + e3z*z2
493 rx(3) = e1x*x3 + e1y*y3 + e1z*z3
494 ry(3) = e2x*x3 + e2y*y3 + e2z*z3
495 rz(3) = e3x*x3 + e3y*y3 + e3z*z3
496 rx(4) = e1x*x4 + e1y*y4 + e1z*z4
497 ry(4) = e2x*x4 + e2y*y4 + e2z*z4
498 rz(4) = e3x*x4 + e3y*y4 + e3z*z4
499
500
501 CALL i2cin_rot27(stbrk,rs,rm,rx(1),ry(1),rz(1),rx(2),ry(2),rz(2),rx(3),ry(3),rz(3),
502 . rx(4),ry(4),rz(4),dpara(1,ii),dwdu,e1x,e1y,e1z,e2x,e2y,e2z,e3x,e3y,e3z,
503 . nir,betax,betay)
504
505
506
507 flx = fs(1)*e1x + fs(2)*e1y + fs(3)*e1z
508 fly = fs(1)*e2x + fs(2)*e2y + fs(3)*e2z
509 flz = fs(1)*e3x + fs(2)*e3y + fs(3)*e3z
510
511 DO j=1,4
512 fmx(j) = h(j)*flx
513 fmy(j) = h(j)*fly
514 fmz(j) = h(j)*flz
515 ENDDO
516
517
518 IF (iroddl==1) THEN
519 mxs = mxi*e1x + myi*e1y + mzi
520 mys = mxi*e2x + myi*e2y + mzi*e2z
521 mzs = mxi*e3x + myi*e3y + mzi*e3z
522
523
525 . fmx ,fmy ,fmz ,h ,stifm ,
526 . mxs ,mys ,mzs ,stifmr ,betax ,
527 . betay)
528 ELSE
529 mxs = zero
530 mys = zero
531 mzs = zero
532
533
535 . fmx ,fmy ,fmz ,h ,stifm ,
536 . mxs ,mys ,mzs ,stifmr ,betax ,
537 . betay)
538 stifmr = zero
539 ENDIF
540
541
542
543 DO j=1,4
544 fx(j) = e1x*fmx(j) + e2x*fmy(j) + e3x*fmz(j)
545 fy(j) = e1y*fmx(j) + e2y*fmy(j) + e3y*fmz(j)
546 fz(j) = e1z*fmx(j) + e2z*fmy(j) + e3z*fmz(j)
547 ENDDO
548
549 ELSE
550
551
552
553
554 stifm=zero
555 stbrk=zero
556 stifmr=zero
557 dwdu=zero
558
559 DO j=1,4
560 fx(j) = fs(1)*h(j)
561 fy(j) = fs(2)*h(j)
562 fz(j) = fs(3)*h(j)
563 ENDDO
564
565 ENDIF
566
567 IF(iroddl/=0)THEN
568 DO jj=1,nir
569 j=irect(jj,l)
570 a(1,j)=a(1,j)+fx(jj)
571 a(2,j)=a(2,j)+fy(jj)
572 a(3,j)=a(3,j)+fz(jj)
573 ms(j)=ms(j)+xmsi*h2(jj)
574 stifn(j)=stifn(j)+weight(i)*(stifn(i)*(one+stbrk)*(abs(h(jj))+stifm)+stifr(i)*stifmr*dwdu)
575 ENDDO
576 ELSE
577 DO jj=1,nir
578 j=irect(jj,l)
579 a(1,j)=a(1,j)+fx(jj)
580 a(2,j)=a(2,j)+fy(jj)
581 a(3,j)=a(3,j)+fz(jj)
582 ms(j)=ms(j)+xmsi*h2(jj)
583 stifn(j)=stifn(j)+weight(i)*(stifn(i)*(one+stbrk)*(abs(h(jj))+stifm))
584 ENDDO
585 END IF
586
587 IF(idtmins==2.OR.idtmins_int/=0) THEN
588
589 t2fac_sms(i) = (one+stbrk)*(one+stifm)
590 ENDIF
591
592
594 . irect(1,l),nir ,fsav ,fncont ,fncontp,
595 . ftcontp ,weight ,h3d_data,i ,h)
596
597 IF(iroddl==0)THEN
598 IF(idel2/=0.AND.ms(i)/=0.) smass(ii)=ms(i)
599 ms(i)=zero
600 stifn(i)=em20
601 a(1,i)=zero
602 a(2,i)=zero
603 a(3,i)=zero
604 ENDIF
605
606 ENDIF
607
608 ENDDO
609 ENDIF
610
611 RETURN
subroutine i2cin_rot27(stbrk, rs, rm, x1, y1, z1, x2, y2, z2, x3, y3, z3, x4, y4, z4, dpara, dwdu, e1x, e1y, e1z, e2x, e2y, e2z, e3x, e3y, e3z, nir, betax, betay)
subroutine i2forces(x, fs, fx, fy, fz, irect, nir, fsav, fncont, fncontp, ftcontp, weight, h3d_data, nsl, h)
subroutine i2loceq_27(nir, rs, rx, ry, rz, fmx, fmy, fmz, h, stifm, mxs, mys, mzs, stifmr, betax, betay)
subroutine i2rep(x1, x2, x3, x4, y1, y2, y3, y4, z1, z2, z3, z4, e1x, e1y, e1z, e2x, e2y, e2z, e3x, e3y, e3z, nir)