56
57
58
59
60
61
62
63
64
65
66
67
68
71 USE intbuf_fric_mod
73 use element_mod , only :nixs,nixc,nixtg,nixt,nixp,nixr
74
75
76
77#include "implicit_f.inc"
78
79
80
81#include "units_c.inc"
82#include "param_c.inc"
83#include "com01_c.inc"
84#include "com04_c.inc"
85#include "scr08_c.inc"
86#include "scr23_c.inc"
87#include "r2r_c.inc"
88
89
90
91 INTEGER,INTENT(IN) :: NRT, NINTR, NTY, NOINT,IGAP,INTTH,INTFRIC,IREM_GAP
92 INTEGER,INTENT(IN) :: IRECT(2,*), IXS(NIXS,NUMELS), IXC(NIXC,NUMELC),
93 . IXTG(NIXTG,NUMELTG),IXT(NIXT,NUMELT),IXP(NIXP,NUMELP),IXR(NIXR,NUMELR),
94 . IPARTC(NUMELC), IPARTTG(NUMELTG),NOD2EL1D(*),KNOD2EL1D(*),ITAB(NUMNOD),
95 . IXS10(6,*),KXX(NIXX,*),IXX(*),IGEO(NPROPGI,NUMGEO),
96 . KNOD2ELS(*), KNOD2ELC(*), KNOD2ELTG(*),
97 . NOD2ELS(*), NOD2ELC(*), NOD2ELTG(*),IWORKSH(3,*),
98 . TAGPRT_FRIC(*),IPARTS(*),
99 . IPARTT(*) ,IPARTP(*) ,IPARTX(*) ,IPARTR(*)
100 INTEGER,INTENT(INOUT) :: IPARTFRIC(*),IELEC(*)
101 INTEGER,INTENT(IN) :: ID,IT19
102 my_real,
INTENT(IN) :: x(3,numnod), pm(npropm,nummat), geo(npropg,numgeo),
103 . ms(*),thk(*),thk_part(*),
104 . lelx(*), fillsol(*),pm_stack(20,*)
105 my_real,
INTENT(IN) :: slsfac,gap0,percent_size
106 my_real,
INTENT(INOUT) :: gap_l(*),stf(*),gap_sm(*),
area(*),drad,gapinf,gapmax,gapmin,bgapsmx
107 CHARACTER(LEN=NCHARTITLE) ,INTENT(IN) :: TITR
108 TYPE(INTBUF_FRIC_STRUCT_),INTENT(IN) :: INTBUF_FRIC_TAB(*)
109
110
111
112 INTEGER NDX, I, INRT, NELS, MT, JJ, JJJ, NELC, J,
113 . MG, NELTG,NELT,NELP,NELR,
114 . IGTYP, IP,N1,N2,K,T,P,R,NELX,IPGMAT,IGMAT,IE,
115 . JJ1,JJ2,IEC,K1,K2,IPL,
116 . IE1(50,2),IE2(50,2),ISUBSTACK,IPG,N3,N4,N5,N6,N7,N8,ICONTR
117 my_real dxm, gapmx, gapmn, areass, vol, dx,gap1,gaps1,gaptmp,xl,
118 . sx1,sy1,sz1,sx2,sy2,sz2,sx3,sy3,sz3,
119 . x1,y1,z1,x2,y2,z2,x3,y3,z3,x4,y4,z4,x5,y5,z5,x6,y6,z6,x7,y7,z7,x8,y8,z8,
120 . xx1(4),xx2(4), xx3(4) ,face(6),
121 . n_1, n_2, n_3,dx1,xl2,bulk
122
124
125
126
127 dxm=zero
128 ndx=0
129 gaps1=zero
130 gapmx=ep30
131 gapmn=ep30
132 ipgmat = 700
133 IF(igap == 3)THEN
134 DO i=1,nrt
135 gap_l(i)=ep30
136 ENDDO
137 ENDIF
138
139 DO i=1,nrt
140 stf(i)=zero
141 gap_sm(i)=zero
142 inrt=i
143 CALL i11gmx3(x,irect,inrt,gapmx,xl2)
144
145
146
147 CALL i11sol(x,irect,ixs,nintr,nels,inrt,areass,noint,knod2els,nod2els,ixs10)
148 IF(nels /= 0) THEN
149 mt=ixs(1,nels)
150 mg=ixs(nixs-1,nels)
151 icontr = igeo(97,mg)
152 IF(mt > 0)THEN
153 DO jj=1,8
154 jjj=ixs(jj+1,nels)
155 xc(jj)=x(1,jjj)
156 yc(jj)=x(2,jjj)
157 zc(jj)=x(3,jjj)
158 ENDDO
160 IF (icontr==1 ) THEN
161 bulk = pm(107,mt)
162 ELSE
163 bulk = pm(32,mt)
164 END IF
165 IF(xl2 > 0.0)THEN
166 stf(i)=slsfac*fillsol(nels)*vol*bulk/xl2
167 ELSE
168 stf(i)=zero
169 ENDIF
170 ELSE
171 IF(nintr >= 0) THEN
172 CALL ancmsg(msgid=95,msgtype=msgwarning,anmode=aninfo_blind_2,
174 . c1=titr,
175 . i2=ixs(nixs,nels),
176 . c2='SOLID',
177 . i3=i)
178 ENDIF
179 IF(nintr < 0) THEN
180 CALL ancmsg(msgid=96,msgtype=msgwarning,anmode=aninfo_blind_2,
182 . c1=titr,
183 . i2=ixs(nixs,nels),
184 . c2='SOLID',
185 . i3=i)
186 ENDIF
187 ENDIF
188
189
190 IF(intfric > 0) THEN
191 ip= iparts(nels)
192 ipg = tagprt_fric(ip)
193 IF(ipg > 0) THEN
195 . ipg,intbuf_fric_tab(intfric)%S_TABPARTS_FRIC,
196 . intbuf_fric_tab(intfric)%TABPARTS_FRIC,ipl )
197 ipartfric(inrt) = ipl
198 ENDIF
199 ENDIF
200
201 ENDIF
202 CALL i11coq(irect ,ixc ,ixtg ,nintr ,nelc ,
203 . neltg ,inrt ,geo ,pm ,thk ,
204 . igeo ,knod2elc,knod2eltg,nod2elc,nod2eltg,
205 . pm_stack,iworksh)
206 IF(neltg /= 0) THEN
207 mt=ixtg(1,neltg)
208 mg=ixtg(5,neltg)
209 igtyp = igeo(11,mg)
210 ip = iparttg(neltg)
211 igmat = igeo(98,mg)
212 IF ( thk_part(ip) /= zero .AND. iintthick == 0) THEN
213 dx=thk_part(ip)
214 ELSEIF (thk(numelc+neltg) /= zero .AND. iintthick ==0)THEN
215 dx=thk(numelc+neltg)
216 ELSEIF(igtyp == 17 .OR. igtyp == 51 .OR. igtyp == 52) THEN
217 dx=thk(numelc+neltg)
218 ELSE
219 dx=geo(1,mg)
220 ENDIF
221 gap_sm(i)=half*dx
222 gaps1=
max(gaps1,gap_sm(i))
223 gapmn =
min(gapmn,dx)
224 dxm=dxm+dx
225 ndx=ndx+1
226 IF(mt > 0)THEN
227 IF(igtyp == 11 .AND. igmat > 0) THEN
228 stf(i)=slsfac*dx*geo(ipgmat + 2 ,mg)
229 ELSEIF(igtyp == 52 .OR. ((igtyp == 17 .OR. igtyp == 51) .AND. igmat > 0)) THEN
230 isubstack = iworksh(3,numelc+neltg)
231 stf(i)=slsfac*dx*pm_stack( 2 ,isubstack)
232 ELSE
233 stf(i)=slsfac*dx*pm(20,mt)
234 ENDIF
235 ELSE
236 IF(nintr >= 0) THEN
237 CALL ancmsg(msgid=95,msgtype=msgwarning,anmode=aninfo_blind_2,
239 . c1=titr,
240 . i2=ixs(nixs,nels),
241 . c2='SOLID',
242 . i3=i)
243 ENDIF
244 IF(nintr < 0) THEN
245 CALL ancmsg(msgid=96,msgtype=msgwarning,anmode=aninfo_blind_2,
247 . c1=titr,
248 . i2=ixs(nixs,nels),
249 . c2='SOLID',
250 . i3=i)
251 ENDIF
252
253 ENDIF
254
255
256 IF(intfric > 0) THEN
257 ip= iparttg(neltg)
258 ipg = tagprt_fric(ip)
259 IF(ipg > 0) THEN
261 . ipg,intbuf_fric_tab(intfric)%S_TABPARTS_FRIC,
262 . intbuf_fric_tab(intfric)%TABPARTS_FRIC,ipl )
263 ipartfric(inrt) = ipl
264 ENDIF
265 ENDIF
266
267
268 ELSEIF(nelc /= 0) THEN
269 mt=ixc(1,nelc)
270 mg=ixc(6,nelc)
271 igtyp = igeo(11,mg)
272 ip = ipartc(nelc)
273 igmat = igeo(98,mg)
274 IF (thk_part(ip) /= zero .AND. iintthick == 0) THEN
275 dx=thk_part(ip)
276 ELSEIF ( thk(nelc) /= zero .AND. iintthick == 0 ) THEN
277 dx=thk(nelc)
278 ELSEIF(igtyp ==17 .OR. igtyp == 51 .OR. igtyp == 52)THEN
279 dx=thk(nelc)
280 ELSE
281 dx=geo(1,mg)
282 ENDIF
283 gap_sm(i)=half*dx
284 gaps1=
max(gaps1,gap_sm(i))
285 gapmn =
min(gapmn,dx)
286 dxm=dxm+dx
287 ndx=ndx+1
288 IF(mt > 0)THEN
289 IF(igtyp == 11 .AND. igmat > 0) THEN
290 stf(i)=slsfac*dx*geo(ipgmat + 2 ,mg)
291 ELSEIF(igtyp == 52 .OR. ((igtyp == 17 .OR. igtyp == 51) .AND. igmat > 0)) THEN
292 isubstack = iworksh(3,nelc)
293 stf(i)=slsfac*dx*pm_stack( 2 ,isubstack)
294 ELSE
295 stf(i)=slsfac*dx*pm(20,mt)
296 ENDIF
297 ELSE
298 IF(nintr >= 0) THEN
299 CALL ancmsg(msgid=95,msgtype=msgwarning,anmode=aninfo_blind_2,
301 . c1=titr,
302 . i2=ixc(nixc,nelc),
303 . c2='SHELL',
304 . i3=i)
305 ENDIF
306 IF(nintr < 0) THEN
307 CALL ancmsg(msgid=96,msgtype=msgwarning, anmode=aninfo_blind_2,
309 . c1=titr,
310 . i2=ixc(nixc,nelc),
311 . c2='SHELL',
312 . i3=i)
313 ENDIF
314 ENDIF
315
316
317 IF(intfric > 0) THEN
318 ip= ipartc(nelc)
319
320 ipg = tagprt_fric(ip)
321 IF(ipg > 0) THEN
323 . ipg,intbuf_fric_tab(intfric)%S_TABPARTS_FRIC,
324 . intbuf_fric_tab(intfric)%TABPARTS_FRIC,ipl )
325 ipartfric(inrt) = ipl
326 ENDIF
327 ENDIF
328
329 ENDIF
330 CALL i11fil(irect,ixt ,ixp,ixr,nintr,nelt ,
331 . nelp,nelr,nelx,inrt,nod2el1d,
332 . knod2el1d,kxx,ixx)
333 IF(nelt /= 0) THEN
334 mt=ixt(1,nelt)
335 mg=ixt(4,nelt)
336 ip = ipartt(nelt)
337 IF (thk_part(ip) > zero ) THEN
338 dx1=thk_part(ip)
339 ELSE
340 dx1=sqrt(geo(1,mg))
341 END IF
342 dx=sqrt(geo(1,mg))
343 gap_sm(i)=
max(gap_sm(i),half*dx1)
344 gaps1=
max(gaps1,gap_sm(i))
345 gapmn =
min(gapmn,dx1)
346 dxm=dxm+dx
347 ndx=ndx+1
348 IF(mt > 0)THEN
349 stf(i)=slsfac*dx*pm(20,mt)
350 IF (nsubdom>0) stf(i)=slsfac*dx*pm_r2r(mt)
351 ELSE
352 IF(nintr >= 0) THEN
353 CALL ancmsg(msgid=95,msgtype=msgwarning,anmode=aninfo_blind_2,
355 . c1=titr,
356 . i2=ixt(nixt,nelt),
357 . c2='TRUSS',
358 . i3=i)
359 ENDIF
360 IF(nintr < 0) THEN
361 CALL ancmsg(msgid=96,msgtype=msgwarning,anmode=aninfo_blind_2,
363 . c1=titr,
364 . i2=ixt(nixt,nelt),
365 . c2='TRUSS',
366 . i3=i)
367 ENDIF
368 ENDIF
369
370 IF(intfric > 0) THEN
371 ip= ipartt(nelt)
372 ipg = tagprt_fric(ip)
373 IF(ipg > 0) THEN
375 . ipg,intbuf_fric_tab(intfric)%S_TABPARTS_FRIC,
376 . intbuf_fric_tab(intfric)%TABPARTS_FRIC,ipl )
377 ipartfric(inrt) = ipl
378 ENDIF
379 ENDIF
380
381 ELSEIF(nelp /= 0) THEN
382 mt=ixp(1,nelp)
383 mg=ixp(5,nelp)
384 ip = ipartp(nelp)
385 IF (thk_part(ip) > zero ) THEN
386 dx1=thk_part(ip)
387 ELSE
388 dx1=sqrt(geo(1,mg))
389 END IF
390 dx=sqrt(geo(1,mg))
391 gap_sm(i)=
max(gap_sm(i),half*dx1)
392 gaps1=
max(gaps1,gap_sm(i))
393 gapmn =
min(gapmn,dx1)
394 dxm=dxm+dx
395 ndx=ndx+1
396 IF(mt > 0)THEN
397 stf(i)=slsfac*dx*pm(20,mt)
398 IF (nsubdom>0) stf(i)=slsfac*dx*pm_r2r(mt)
399 ELSE
400 IF(nintr >= 0) THEN
401 CALL ancmsg(msgid=95,msgtype=msgwarning,anmode=aninfo_blind_2,
403 . c1=titr,
404 . i2=ixp(nixp,nelp),
405 . c2='BEAM',
406 . i3=i)
407 ENDIF
408 IF(nintr<0) THEN
409 CALL ancmsg(msgid=96,msgtype=msgwarning,anmode=aninfo_blind_2,
411 . c1=titr,
412 . i2=ixp(nixp,nelp),
413 . c2='BEAM',
414 . i3=i)
415 ENDIF
416 ENDIF
417
418 IF(intfric > 0) THEN
419 ip= ipartp(nelp)
420 ipg = tagprt_fric(ip)
421 IF(ipg > 0) THEN
423 . ipg,intbuf_fric_tab(intfric)%S_TABPARTS_FRIC,
424 . intbuf_fric_tab(intfric)%TABPARTS_FRIC,ipl )
425 ipartfric(inrt) = ipl
426 ENDIF
427 ENDIF
428
429 ELSEIF(nelr /= 0) THEN
430 mg=ixr(1,nelr)
431 mt = ixr(5,nelr)
432 ip = ipartr(nelr)
433 IF(ip > 0)THEN
434 IF (thk_part(ip) > zero ) THEN
435 dx1=thk_part(ip)
436 gap_sm(i)=
max(gap_sm(i),half*dx1)
437 gaps1=
max(gaps1,gap_sm(i))
438 gapmn =
min(gapmn,dx1)
439 END IF
440 ENDIF
441 IF(mg > 0)THEN
442 igtyp=nint(geo(12,mg))
443 IF(igtyp==4.OR.igtyp==12)THEN
444 stf(i)=slsfac*geo(2,mg)
445 ELSEIF(igtyp==8.OR.igtyp==13)THEN
446 stf(i)=slsfac*
max(geo(3,mg),geo(10,mg),geo(15,mg))
447 ELSEIF(igtyp == 23)THEN
448 stf(i)=slsfac*
max(pm(191,mt),pm(192,mt),pm(193,mt))
449 ELSEIF(igtyp==25)THEN
450 stf(i)=slsfac*geo(10,mg)
451 ELSEIF(igtyp>=29)THEN
452 stf(i)=slsfac*geo(3,mg)
453 ELSE
454 WRITE(6,'(A)') 'INTERNAL ERROR 987'
456 ENDIF
457 ELSE
458 IF(nintr >= 0) THEN
459 CALL ancmsg(msgid=95,msgtype=msgwarning,anmode=aninfo_blind_2,
461 . c1=titr,
462 . i2=ixr(nixr,nelr),
463 . c2='SPRING',
464 . i3=i)
465 ENDIF
466 IF(nintr < 0) THEN
467 CALL ancmsg(msgid=96,msgtype=msgwarning,anmode=aninfo_blind_2,
469 . c1=titr,
470 . i2=ixr(nixr,nelr),
471 . c2='SPRING',
472 . i3=i)
473 ENDIF
474 ENDIF
475
476 IF(intfric > 0) THEN
477 ip= ipartr(nelr)
478 ipg = tagprt_fric(ip)
479 IF(ipg > 0) THEN
481 . ipg,intbuf_fric_tab(intfric)%S_TABPARTS_FRIC,
482 . intbuf_fric_tab(intfric)%TABPARTS_FRIC,ipl )
483 ipartfric(inrt) = ipl
484 ENDIF
485 ENDIF
486
487 ELSEIF(nelx /= 0) THEN
488 mg=kxx(2,nelx)
489 IF(mg>0)THEN
490 stf(i)=slsfac*get_u_geo(4,mg)*(kxx(3,nelx)-1)/lelx(nelx)
491 ELSE
492 IF(nintr >= 0) THEN
493 CALL ancmsg(msgid=95,msgtype=msgwarning,anmode=aninfo_blind_2,
495 . c1=titr,
496 . i2=kxx(nixx,nelx),
497 . c2='XELEM',
498 . i3=i)
499 ENDIF
500 IF(nintr < 0) THEN
501 CALL ancmsg(msgid=96,msgtype=msgwarning,anmode=aninfo_blind_2,
503 . c1=titr,
504 . i2=kxx(nixx,nelx),
505 . c2='XELEM',
506 . i3=i)
507 ENDIF
508 ENDIF
509
510 IF(intfric > 0) THEN
511 ip= ipartx(nelx)
512 ipg = tagprt_fric(ip)
513 IF(ipg > 0) THEN
515 . ipg,intbuf_fric_tab(intfric)%S_TABPARTS_FRIC,
516 . intbuf_fric_tab(intfric)%TABPARTS_FRIC,ipl )
517 ipartfric(inrt) = ipl
518 ENDIF
519 ENDIF
520
521 ENDIF
522
523 IF(nels+nelc+neltg+nelt+nelp+nelr+numelx==0.)THEN
524 IF(nintr > 0) THEN
525 CALL ancmsg(msgid=481,msgtype=msgerror,anmode=aninfo_blind_2,
527 . c1=titr,
528 . i2=i)
529 ENDIF
530 IF(nintr < 0) THEN
531 CALL ancmsg(msgid=482,msgtype=msgerror,anmode=aninfo_blind_2,
533 . c1=titr,
534 . i2=i)
535 ENDIF
536 ENDIF
537 enddo
538
539
540
541
542 IF(igap == 3)THEN
543 DO i=1,nrt
544 xl = ep30
545 n1=irect(1,i)
546 n2=irect(2,i)
547 IF(n1 /= n2 .AND. n1 /= 0)
548 . xl=
min(xl,sqrt((x(1,n1)-x(1,n2))**2+(x(2,n1)-x(2,n2))**2+
549 . (x(3,n1)-x(3,n2))**2))
550
551 DO j=1,2
552 n1=irect(j,i)
553 DO k=knod2el1d(n1)+1,knod2el1d(n1+1)
554 IF (nod2el1d(k) <= numelt .AND. nod2el1d(k) /= zero) THEN
555 t = nod2el1d(k)
556 xl=
min(xl,sqrt((x(1,ixt(2,t))-x(1,ixt(3,t)))**2+
557 . (x(2,ixt(2,t))-x(2,ixt(3,t)))**2+
558 . (x(3,ixt(2,t))-x(3,ixt(3,t)))**2))
559 ELSEIF (nod2el1d(k) <= numelt+numelp .AND. nod2el1d(k) /= zero) THEN
560 p = nod2el1d(k) - numelt
561 xl=
min(xl,sqrt((x(1,ixp(2,p))-x(1,ixp(3,p)))**2+
562 . (x(2,ixp(2,p))-x(2,ixp(3,p)))**2+
563 . (x(3,ixp(2,p))-x(3,ixp(3,p)))**2))
564 ELSEIF (nod2el1d(k) <= numelt+numelp+numelr .AND. nod2el1d(k) /= zero) THEN
565 r = nod2el1d(k) - numelt - numelp
566 xl=
min(xl,sqrt((x(1,ixr(2,r))-x(1,ixr(3,r)))**2+
567 . (x(2,ixr(2,r))-x(2,ixr(3,r)))**2+
568 . (x(3,ixr(2,r))-x(3,ixr(3,r)))**2))
569 ENDIF
570 ENDDO
571 ENDDO
572 DO j=1,2
573 gap_l(i) =
min(gap_l(i),percent_size*xl)
574 ENDDO
575 ENDDO
576 ENDIF
577
578
579
580 gapmx=sqrt(gapmx)
581 IF(igap==0)THEN
582
583
584
585 IF(gap0 > zero)THEN
586 gap1 = gap0
587 ELSE
588 IF(ndx/=0)THEN
589 gap1 =
min(half*gapmx,dxm/ndx)
590 ELSE
591 gap1 = em01* gapmx
592 ENDIF
593 IF ((nintr<0).AND.(it19==0)) WRITE(iout,1300)half*(gapmin+gap1)
594 ENDIF
595
596 IF(nintr < 0) gap1 = half*(gapmin+gap1)
597 gapmin = gap1
598 gapmax = gap1
599
600 IF ((gap1 > half*gapmx) .AND. (irem_gap /= 2)) THEN
601 gaptmp = half*gapmx
602 CALL ancmsg(msgid=94,msgtype=msgwarning,anmode=aninfo_blind_2,
604 . c1=titr,
605 . r1=gap1,
606 . r2=gaptmp)
607 ENDIF
608 ELSE
609
610
611
612 bgapsmx = zero
613 IF(gap0 > zero)THEN
614 gap1 = gap0
615 ELSE
616 IF(ndx /= 0)THEN
617 gap1 =
min(half*gapmx,gapmn)
618 ELSE
619 gap1 = em01 * gapmx
620 ENDIF
621 IF ((nintr<0) .AND. (it19==0)) WRITE(iout,1300)half
622 ENDIF
623
624 IF(nintr>0)THEN
625 gapmin = gap1
626 gapmax = gaps1
627 ELSE
628 gapmin = half*(gapmin+gap1)
629 gapmax =
max(gapmax+gaps1,gapmin)
630 bgapsmx =
max(bgapsmx,gaps1)
631 ENDIF
632
633 IF ((gapmax>half*gapmx) .AND. (igap/=3) .AND. (irem_gap/=2))THEN
634 gaptmp = half*gapmx
635 CALL ancmsg(msgid=94,msgtype=msgwarning,anmode=aninfo_blind_2,
637 . c1=titr,
638 . r1=gapmax,
639 . r2=gaptmp)
640 ENDIF
641 ENDIF
642
643
644
645 IF(slsfac < zero)THEN
646 DO i=1,nrt
647 stf(i)=-slsfac
648 ENDDO
649 ENDIF
650
651
652
653
654 IF(igap == 0) THEN
655 gapinf=gapmax
656 ELSEIF(igap==1 .OR. igap==2) THEN
657 DO i = 1, nrt
658 gapinf =
min(gapinf,gap_sm(i))
659 ENDDO
660 ELSEIF(igap==3) THEN
661 DO i = 1, nrt
662 gapinf =
min(gapinf,
min(gap_sm(i),gap_l(i)))
663 ENDDO
664 ENDIF
665
666 IF(intth /= 0)THEN
667 IF(drad == zero)THEN
668
670 ELSEIF(drad < gap1)THEN
671
672 drad=gap1
673 END IF
674 WRITE(iout,2001)drad
675
676
677 IF(drad > gapmx)THEN
678 CALL ancmsg(msgid=918, msgtype=msgwarning, anmode=aninfo_blind_2,
680 . c1=titr,
681 . r1=drad ,
682 . r2=gapmx,
684 END IF
685 END IF
686
687
688 IF(intth > 0 ) THEN
689 IF(nelc /=0 .OR. neltg /= 0) THEN
690 DO i=1,nrt
692 jj1 = 0
693 DO j=knod2elc(irect(1,i))+1,knod2elc(irect(1,i)+1)
694 jj1 = jj1 +1
695 ie1(jj1,1) = nod2elc(j)
696 ie1(jj1,2) = 4
697 ENDDO
698 DO j= knod2eltg(irect(1,i))+1,knod2eltg(irect(1,i)+1)
699 jj1 = jj1 +1
700 ie1(jj1,1) = nod2eltg(j)
701 ie1(jj1,2) = 3
702 ENDDO
703 jj2 = 0
704 DO j=knod2elc(irect(2,i))+1,knod2elc(irect(2,i)+1)
705 jj2 = jj2 +1
706 ie2(jj2,1) = nod2elc(j)
707 ie2(jj2,2) = 4
708 ENDDO
709 DO j= knod2eltg(irect(2,i))+1,knod2eltg(irect(2,i)+1)
710 jj2 = jj2 +1
711 ie2(jj2,1) = nod2eltg(j)
712 ie2(jj2,2) = 3
713 ENDDO
714 iec = 0
715 DO j=1,jj1
716 DO k=1,jj2
717 IF(ie1(j,1) == ie2(k,1)) THEN
718 ie = ie1(j,1)
719 iec = iec +1
720 IF(ie1(j,2)==4) THEN
721 sx1 = x(1,ixc(4,ie)) - x(1,ixc(2,ie))
722 sy1 = x(2,ixc(4,ie)) - x(2,ixc(2,ie))
723 sz1 = x(3,ixc(4,ie)) - x(3,ixc(2,ie))
724 sx2 = x(1,ixc(5,ie)) - x(1,ixc(3,ie))
725 sy2 = x(2,ixc(5,ie)) - x(2,ixc(3,ie))
726 sz2 = x(3,ixc(5,ie)) - x(3,ixc(3,ie))
727 sx3 = sy1*sz2 - sz1*sy2
728 sy3 = sz1*sx2 - sx1*sz2
729 sz3 = sx1*sy2 - sy1*sx2
730 area(i) =
area(i)+ one_over_8*sqrt(sx3*sx3+sy3*sy3+sz3*sz3)
731 ielec(i) = ixc(1,ie)
732 mg = ixc(6,ie)
733 ENDIF
734 IF(ie1(j,2) == 3) THEN
735 sx1 = x(1,ixtg(3,ie)) - x(1,ixtg(2,ie))
736 sy1 = x(2,ixtg(3,ie)) - x(2,ixtg(2,ie))
737 sz1 = x(3,ixtg(3,ie)) - x(3,ixtg(2,ie))
738 sx2 = x(1,ixtg(4,ie)) - x(1,ixtg(2,ie))
739 sy2 = x(2,ixtg(4,ie)) - x(2,ixtg(2,ie))
740 sz2 = x(3,ixtg(4,ie)) - x(3,ixtg(2,ie))
741 sx3 = sy1*sz2 - sz1*sy2
742 sy3 = sz1*sx2 - sx1*sz2
743 sz3 = sx1*sy2 - sy1*sx2
744 area(i) =
area(i) + one_over_6*sqrt(sx3*sx3+sy3*sy3+sz3*sz3)
745 ielec(i) = ixtg(1,ie)
746 mg = ixtg(5,ie)
747 ENDIF
748 ENDIF
749 ENDDO
750 ENDDO
752 IF(iec == 1) THEN
753 dx=geo(1,mg)
755 ENDIF
756 ENDDO
757
758 ELSEIF(nels /= 0) THEN
759 DO i=1,nrt
761 jj1 = 0
762 DO j=knod2els(irect(1,i))+1,knod2els(irect(1,i)+1)
763 jj1 = jj1 +1
764 ie1(jj1,1) = nod2els(j)
765 ENDDO
766 jj2 = 0
767 DO j=knod2els(irect(2,i))+1,knod2els(irect(2,i)+1)
768 jj2 = jj2 +1
769 ie2(jj2,1) = nod2els(j)
770 ENDDO
771 iec = 0
772 DO j=1,jj1
773 DO k=1,jj2
774 IF(ie1(j,1) == ie2(k,1)) THEN
775 ie = ie1(j,1)
776 iec= iec +1
777 ielec(i) = ixs(1,ie)
778
779 n1=ixs(2,ie)
780 n2=ixs(3,ie)
781 n3=ixs(4,ie)
782 n4=ixs(5,ie)
783 n5=ixs(6,ie)
784 n6=ixs(7,ie)
785 n7=ixs(8,ie)
786 n8=ixs(9,ie)
787
788 x1=x(1,n1)
789 y1=x(2,n1)
790 z1=x(3,n1)
791 x2=x(1,n2)
792 y2=x(2,n2)
793 z2=x(3,n2)
794 x3=x(1,n3)
795 y3=x(2,n3)
796 z3=x(3,n3)
797 x4=x(1,n4)
798 y4=x(2,n4)
799 z4=x(3,n4)
800 x5=x(1,n5)
801 y5=x(2,n5)
802 z5=x(3,n5)
803 x6=x(1,n6)
804 y6=x(2,n6)
805 z6=x(3,n6)
806 x7=x(1,n7)
807 y7=x(2,n7)
808 z7=x(3,n7)
809 x8=x(1,n8)
810 y8=x(2,n8)
811 z8=x(3,n8)
812
813
814 xx1(1)=x1
815 xx2(1)=y1
816 xx3(1)=z1
817 xx1(2)=x2
818 xx2(2)=y2
819 xx3(2)=z2
820 xx1(3)=x3
821 xx2(3)=y3
822 xx3(3)=z3
823 xx1(4)=x4
824 xx2(4)=y4
825 xx3(4)=z4
826 CALL norma1(n_1,n_2,n_3,face(1),xx1,xx2,xx3)
827 IF( n4/=n3
828 . .AND.n3/=n2
829 . .AND.n2/=n1
830 . .AND.n1/=n4)THEN
831 face(1) = fourth*face(1)
832 ELSE
833 face(1) = third*face(1)
834 ENDIF
835
836 xx1(1)=x5
837 xx2(1)=y5
838 xx3(1)=z5
839 xx1(2)=x6
840 xx2(2)=y6
841 xx3(2)=z6
842 xx1(3)=x7
843 xx2(3)=y7
844 xx3(3)=z7
845 xx1(4)=x8
846 xx2(4)=y8
847 xx3(4)=z8
848 CALL norma1(n_1,n_2,n_3,face(2),xx1,xx2,xx3)
849 IF( n8/=n7
850 . .AND.n7/=n6
851 . .AND.n6/=n5
852 . .AND.n5/=n8)THEN
853 face(2) = fourth*face(2)
854 ELSE
855 face(2) = third*face(2)
856 ENDIF
857
858 xx1(1)=x2
859 xx2(1)=y2
860 xx3(1)=z2
861 xx1(2)=x3
862 xx2(2)=y3
863 xx3(2)=z3
864 xx1(3)=x7
865 xx2(3)=y7
866 xx3(3)=z7
867 xx1(4)=x6
868 xx2(4)=y6
869 xx3(4)=z6
870 CALL norma1(n_1,n_2,n_3,face(3),xx1,xx2,xx3)
871 IF( n6/=n7
872 . .AND.n7/=n3
873 . .AND.n3/=n2
874 . .AND.n2/=n6)THEN
875 face(3) = fourth*face(3)
876 ELSE
877 face(3) = third*face(3)
878 ENDIF
879
880 xx1(1)=x1
881 xx2(1)=y1
882 xx3(1)=z1
883 xx1(2)=x4
884 xx2(2)=y4
885 xx3(2)=z4
886 xx1(3)=x8
887 xx2(3)=y8
888 xx3(3)=z8
889 xx1(4)=x5
890 xx2(4)=y5
891 xx3(4)=z5
892 CALL norma1(n_1,n_2,n_3,face(4),xx1,xx2,xx3)
893 IF( n5/=n8
894 . .AND.n8/=n4
895 . .AND.n4/=n1
896 . .AND.n1/=n5)THEN
897 face(4) = fourth*face(4)
898 ELSE
899 face(4) = third*face(4)
900 ENDIF
901
902 xx1(1)=x1
903 xx2(1)=y1
904 xx3(1)=z1
905 xx1(2)=x2
906 xx2(2)=y2
907 xx3(2)=z2
908 xx1(3)=x6
909 xx2(3)=y6
910 xx3(3)=z6
911 xx1(4)=x5
912 xx2(4)=y5
913 xx3(4)=z5
914 CALL norma1(n_1,n_2,n_3,face(5),xx1,xx2,xx3)
915 IF( n5/=n6
916 . .AND.n6/=n2
917 . .AND.n2/=n1
918 . .AND.n1/=n5)THEN
919 face(5) = fourth*face(5)
920 ELSE
921 face(5) = third*face(5)
922 ENDIF
923
924 xx1(1)=x4
925 xx2(1)=y4
926 xx3(1)=z4
927 xx1(2)=x3
928 xx2(2)=y3
929 xx3(2)=z3
930 xx1(3)=x7
931 xx2(3)=y7
932 xx3(3)=z7
933 xx1(4)=x8
934 xx2(4)=y8
935 xx3(4)=z8
936 CALL norma1(n_1,n_2,n_3,face(6),xx1,xx2,xx3)
937 IF( n8/=n7
938 . .AND.n7/=n3
939 . .AND.n3/=n4
940 . .AND.n4/=n8)THEN
941 face(6) = fourth*face(6)
942 ELSE
943 face(6) = third*face(6)
944 ENDIF
945
946 DO k1=1,8
947 n1 = ixs(k1+1,ie)
948 IF (n1 == irect(1,i).AND.n2 == irect(2,i)) THEN
949 DO k2=1,8
950 n2 = ixs(k2+1,ie)
951 IF (k1 == 1 .AND. k2 == 2) THEN
952 area(i) =
area(i) + (face(1)+face(5))
953 ELSEIF (k1 == 1 .AND. k2 == 4) THEN
954 area(i) =
area(i) + (face(1)+face(4))
955 ELSEIF (k1 == 1 .AND. k2 == 5) THEN
956 area(i) =
area(i) + (face(4)+face(5))
957 ELSEIF (k1 == 2 .AND. k2 == 3) THEN
958 area(i) =
area(i) + (face(1)+face(3))
959 ELSEIF (k1 == 2 .AND. k2 == 6) THEN
960 area(i) =
area(i) + (face(5)+face(3))
961 ELSEIF (k1 == 3 .AND. k2 == 4) THEN
962 area(i) =
area(i) + (face(1)+face(6))
963 ELSEIF (k1 == 3 .AND. k2 == 7) THEN
964 area(i) =
area(i) + (face(3)+face(6))
965 ELSEIF (k1 == 4 .AND. k2 == 1) THEN
966 area(i) =
area(i) + (face(1)+face(4))
967 ELSEIF (k1 == 4 .AND. k2 == 8) THEN
968 area(i) =
area(i) + (face(6)+face(4))
969 ELSEIF (k1 == 5 .AND. k2 == 6) THEN
970 area(i) =
area(i) + (face(2)+face(5))
971 ELSEIF (k1 == 6 .AND. k2 == 7) THEN
972 area(i) =
area(i) + (face(2)+face(5))
973 ELSEIF (k1 == 7 .AND. k2 == 8) THEN
974 area(i) =
area(i) + (face(2)+face(6))
975 ELSEIF (k1 == 8.AND. k2 == 5) THEN
976 area(i) =
area(i) + (face(2)+face(4))
977 ENDIF
978 ENDDO
979 ENDIF
980 ENDDO
981 ENDIF
982 ENDDO
983 ENDDO
984 ENDDO
985 ELSEIF(nelp /= 0) THEN
986 DO i=1,nrt
988 jj1 = 0
989 DO j=knod2el1d(irect(1,i))+1,knod2el1d(irect(1,i)+1)
990 ie = nod2el1d(j)
991 jj1 = jj1+1
992 mg = ixp(5,ie)
993 dx = sqrt(geo(1,mg))
995 ielec(i) = ixp(1,ie)
996 ENDDO
997 DO j=knod2el1d(irect(2,i))+1,knod2el1d(irect(2,i)+1)
998 ie = nod2el1d(j)
999 jj1 = jj1+1
1000 mg = ixp(5,ie)
1001 dx = sqrt(geo(1,mg))
1003 ENDDO
1005 ENDDO
1006 ELSEIF(nelt /= 0) THEN
1007 DO i=1,nrt
1009 jj1 = 0
1010 DO j=knod2el1d(irect(1,i))+1,knod2el1d(irect(1,i)+1)
1011 ie = nod2el1d(j)
1012 jj1 = jj1+1
1013 mg=ixt(4,ie)
1014 dx=sqrt(geo(1,mg))
1016 ielec(i) = ixt(1,ie)
1017 ENDDO
1018 DO j=knod2el1d(irect(2,i))+1,knod2el1d(irect(2,i)+1)
1019 ie = nod2el1d(j)
1020 jj1 = jj1+1
1021 mg = ixt(4,ie)
1022 dx = sqrt(geo(1,mg))
1024 ENDDO
1026 ENDDO
1027 ENDIF
1028
1029 ENDIF
1030
1031
1032 RETURN
1033
1034
1035 1300 FORMAT(2x,'COMPUTED GAP = ',1pg20.13)
1036 2001 FORMAT(2x,'Maximum distance for radiation computation = ',1pg20.13)
1037
subroutine area(d1, x, x2, y, y2, eint, stif0)
subroutine i11coq(irect, ixc, ixtg, nint, nel, neltg, is, geo, pm, thk, igeo, knod2elc, knod2eltg, nod2elc, nod2eltg, pm_stack, iworksh)
subroutine i11fil(irect, ixt, ixp, ixr, nint, nelt, nelp, nelr, nelx, is, nod2el1d, knod2el1d, kxx, ixx)
subroutine i11gmx3(x, irect, i, gapmax, xl2)
subroutine i11sol(x, irect, ixs, nint, nel, i, area, noint, knod2els, nod2els, ixs10)
subroutine friction_parts_search(ip, npartsfric, partsfric, ipl)
integer, parameter nchartitle
subroutine norma1(n1, n2, n3, area, xx1, xx2, xx3)
subroutine ancmsg(msgid, msgtype, anmode, i1, i2, i3, i4, i5, i6, i7, i8, i9, i10, i11, i12, i13, i14, i15, i16, i17, i18, i19, i20, r1, r2, r3, r4, r5, r6, r7, r8, r9, c1, c2, c3, c4, c5, c6, c7, c8, c9, prmode)