60
61
62
63 USE my_alloc_mod
65 USE intbuf_fric_mod
68 use element_mod , only :nixs,nixc,nixtg,nixt,nixp,nixr
69
70
71
72#include "implicit_f.inc"
73
74
75
76#include "com01_c.inc"
77#include "com04_c.inc"
78#include "param_c.inc"
79#include "remesh_c.inc"
80#include "scr08_c.inc"
81#include "scr17_c.inc"
82#include "units_c.inc"
83
84
85
86 INTEGER NRT, NINT, NTY, NOINT,NSN,IGAP,
87 . INACTI,INTFRIC, NRT_IGE
88 INTEGER IRECT(4,*), IXS(NIXS,*), IXC(NIXC,*),
89 . NSV(*), IXTG(NIXTG,*), IXT(NIXT,*), IXP(NIXP,*),
90 . KNOD2ELS(*), KNOD2ELC(*), KNOD2ELTG(*), NOD2ELS(*), NOD2ELC(*),
91 . NOD2ELTG(*), IELES(*), INTTH, IELEC(*),
92 . SH3TREE(KSH3TREE,*), SH4TREE(KSH4TREE,*),IXR(NIXR,*) ,
93 . IPART(LIPART1,*),IPARTC(*),IPARTTG(*),NOD2EL1D(*),KNOD2EL1D(*),
94 . ITAB(*), IXS10(6,*), IXS16(*), IXS20(*),IDDLEVEL,IGEO(NPROPGI,*),
95 . IWORKSH(3,*),IT19,KXIG3D(NIXIG3D,*),IXIG3D(*),TAGPRT_FRIC(*),
96 . IPARTFRICS(*),IPARTFRICM(*),IPARTS(*),IREM_GAP
98 . stfac, gap,gapmin,gapinf, gapmax,percent_size, bgapsmx,
99 . gapinf_s, gapinf_m, drad, gapm_mx, gaps_mx, gaps_l_mx, gapm_l_mx
101 . x(3,*), stf(*), pm(npropm,*), geo(npropg,*), stfn(*),
102 . ms(*),wa(*),gap_s(*),gap_m(*),
103 . areas(*),thk(*),thk_part(*),
104 . gap_s_l(*),gap_m_l(*), fillsol(*),pm_stack(20,*)
105 INTEGER, DIMENSION(NUMELT), INTENT(IN) :: IPARTT
106 INTEGER, DIMENSION(NUMELP), INTENT(IN) :: IPARTP
107 INTEGER, DIMENSION(NUMELR), INTENT(IN) :: IPARTR
108 INTEGER, INTENT(IN) :: ID
109 CHARACTER(LEN=NCHARTITLE)::TITR
110 TYPE(INTBUF_FRIC_STRUCT_) INTBUF_FRIC_TAB(*)
111 TYPE (SURF_) :: IGRSURF
112 INTEGER, DIMENSION(NUMELS), INTENT(INOUT):: ELEM_LINKED_TO_SEGMENT
113 INTEGER, INTENT(IN) :: FLAG_ELEM_INTER25(NINTER25,NUMELS)
114
115
116
117 INTEGER NDX, I, J, INRT, NELS, MT, JJ, JJJ, NELC,
118 . MG, L, NELTG,N1,N2,IE,
119 . IP, NLEV, MYLEV, K, P, R, T,IGTYP,IPGMAT,IGMAT,
120 . ISUBSTACK,NELIG3D, COIN_IGE(8), IPID, PX, PY, PZ ,IPFMAX,,
121 . IPFLMAX,IPG,NINV,ICONTR,NIN25
122 INTEGER JPERM(4)
123 LOGICAL TYPE18
125 . dxm, gapmx, gapmn,
area, vol, dx, gapm,
126 . gaptmp, gapscale,sx1,sy1,sz1,sx2,sy2,sz2,sx3,sy3,sz3,
127 . slsfac,xl,bulk
128 my_real,
dimension(:),
allocatable :: gap_s_l_tmp
129 INTEGER, DIMENSION(:),ALLOCATABLE :: TAGELEMS,INDEXE
130 LOGICAL :: PRINT_ERROR
131 INTEGER, DIMENSION(4) :: NODE_ID
132 INTEGER :: IELEM(2)
133 DATA jperm/2,3,4,1/
134
135
136
137
138
139
140
141 ALLOCATE( gap_s_l_tmp(numnod) )
142 type18 = .false.
143 IF(inacti == 7)type18=.true.
144 ipgmat = 700
145 IF(nty == 20)THEN
146 slsfac = one
147 ELSE
148 slsfac = stfac
149 ENDIF
150 IF(igap == 3)THEN
151 DO i=1,numnod
152 gap_s_l_tmp(i)=zero
153 ENDDO
154 DO i=1,nrt+nrt_ige
155 gap_m_l(i)=ep30
156 ENDDO
157 DO i=1,nsn
158 gap_s_l(i)=ep30
159 ENDDO
160 ENDIF
161 dxm=0.
162 ndx=0
163 gapmx=ep30
164 gapmn=ep30
165 gaps_mx=zero
166 gaps_l_mx=zero
167 gapm_mx=zero
168 gapm_l_mx=zero
169
170 IF(igap == 2)THEN
171 IF(iddlevel == 1) igap = 1
172 gapscale = gapmin
173 gapmin = zero
174 ELSEIF(igap == 3)THEN
175 gapscale=gapmin
176 gapmin=zero
177 ELSE
178 gapscale = one
179 ENDIF
180
181 IF(igap == 3)THEN
182 DO i=1,nrt+nrt_ige
183 xl = ep30
184 DO j=1,4
185 n1=irect(j,i)
186 n2=irect(jperm(j),i)
187 IF(n1 /= n2 .AND. n1 /= 0)xl=
min(xl,sqrt((x(1,n1)-x(1,n2))**2+(x(2,n1)-x(2,n2))**2+(x(3,n1)-x(3,n2))**2))
188 ENDDO
189 DO j=1,4
190 IF(gap_s_l_tmp(irect(j,i)) == zero) THEN
191 gap_s_l_tmp(irect(j,i))= percent_size*xl
192 ELSE
193 gap_s_l_tmp(irect(j,i))=
min(gap_s_l_tmp(irect(j,i)),percent_size*xl)
194 ENDIF
195 ENDDO
196
197 DO j=1,4
198 n1=irect(j,i)
199 DO k=knod2el1d(n1)+1,knod2el1d(n1+1)
200 IF (nod2el1d(k) <= numelt .AND. nod2el1d(k) /= zero) THEN
201 t=nod2el1d(k)
202 xl=
min(xl,sqrt((x(1,ixt(2,t))-x(1,ixt(3,t)))**2 + (x(2,ixt(2,t))-x(2,ixt(3,t)))**2 + (x(3,ixt(2,t))-x(3,ixt(3,t)))**2))
203 ELSEIF (nod2el1d(k) <= numelt+numelp .AND. nod2el1d(k) /= zero) THEN
204 p=nod2el1d(k) - numelt
205 xl=
min(xl,sqrt((x(1,ixp(2,p))-x(1,ixp(3,p)))**2 + (x(2,ixp(2,p))-x(2,ixp(3,p)))**2 + (x(3,ixp(2,p))-x(3,ixp(3,p)))**2))
206 ELSEIF (nod2el1d(k) <= numelt+numelp+numelr .AND. nod2el1d(k) /= zero) THEN
207 r=nod2el1d(k) - numelt - numelp
208 xl=
min(xl,sqrt((x(1,ixr(2,r))-x(1,ixr(3,r)))**2 + (x(2,ixr(2,r))-x(2,ixr(3,r)))**2 + (x(3,ixr(2,r))-x(3,ixr(3,r)))**2))
209 ENDIF
210 ENDDO
211 ENDDO
212 gap_m_l(i)=percent_size*xl
213 gapm_l_mx =
max(gapm_l_mx,gap_m_l(i))
214 DO j=1,4
215 gap_s_l_tmp(irect(j,i))=
min(gap_s_l_tmp(irect(j,i)),percent_size*xl)
216 ENDDO
217 ENDDO
218 ENDIF
219
220
221
222 IF(igap >= 1)THEN
223 DO i=1,numnod
224 wa(i)=zero
225 ENDDO
226 DO i=1,numelc
227 mg=ixc(6,i)
228 igtyp = igeo(11,mg)
229 ip = ipartc(i)
230 IF ( thk_part(ip) /= zero .AND. iintthick == 0) THEN
231 dx=half*thk_part(ip)
232 ELSEIF ( thk(i) /= zero .AND. iintthick == 0) THEN
233 dx=half*thk(i)
234 ELSEIF(igtyp == 17 .OR. igtyp ==51 .OR. igtyp == 52) THEN
235 dx=half*thk(i)
236 ELSE
237 dx=half*geo(1,mg)
238 ENDIF
239 wa(ixc(2,i))=
max(wa(ixc(2,i)),dx)
240 wa(ixc(3,i))=
max(wa(ixc(3,i)),dx)
241 wa(ixc(4,i))=
max(wa(ixc(4,i)),dx)
242 wa(ixc(5,i))=
max(wa(ixc(5,i)),dx)
243 ENDDO
244 DO i=1,numeltg
245 mg=ixtg(5,i)
246 igtyp = igeo(11,mg)
247 ip = iparttg(i)
248 IF ( thk_part(ip) /= zero .AND. iintthick == 0) THEN
249 dx=half*thk_part(ip)
250 ELSEIF ( thk(numelc+i) /= zero .AND. iintthick == 0) THEN
251 dx=half*thk(numelc+i)
252 ELSEIF(igtyp == 17 .OR. igtyp ==51 .OR. igtyp == 52) THEN
253 dx=half*thk(numelc+i)
254 ELSE
255 dx=half*geo(1,mg)
256 ENDIF
257 wa(ixtg(2,i))=
max(wa(ixtg(2,i)),dx)
258 wa(ixtg(3,i))=
max(wa(ixtg(3,i)),dx)
259 wa(ixtg(4,i))=
max(wa(ixtg(4,i)),dx)
260 ENDDO
261 DO i=1,numelt
262 mg=ixt(4,i)
263 ip = ipartt(i)
264 IF ( thk_part(ip) > zero ) THEN
265 dx=half*thk_part(ip)
266 ELSE
267 dx=half*sqrt(geo(1,mg))
268 END IF
269 wa(ixt(2,i))=
max(wa(ixt(2,i)),dx)
270 wa(ixt(3,i))=
max(wa(ixt(3,i)),dx)
271 ENDDO
272 DO i=1,numelp
273 mg=ixp(5,i)
274 ip = ipartp(i)
275 IF ( thk_part(ip) > zero ) THEN
276 dx=half*thk_part(ip)
277 ELSE
278 dx=half*sqrt(geo(1,mg))
279 END IF
280 wa(ixp(2,i))=
max(wa(ixp(2,i)),dx)
281 wa(ixp(3,i))=
max(wa(ixp(3,i)),dx)
282 ENDDO
283 DO i=1,numelr
284 ip = ipartr(i)
285 IF ( thk_part(ip) > zero ) THEN
286 mg=ixr(1,i)
287 igtyp = igeo(11,mg)
288 dx=half*thk_part(ip)
289 wa(ixr(2,i))=
max(wa(ixr(2,i)),dx)
290 wa(ixr(3,i))=
max(wa(ixr(3,i)),dx)
291 IF (igtyp==12) wa(ixr(4,i))=
max(wa(ixr(4,i)),dx)
292 END IF
293 ENDDO
294 IF(type18)THEN
295 gaps_mx = zero
296 gap_s(1:nsn) = zero
297 ELSE
298 DO i=1,nsn
299 gap_s(i)=gapscale * wa(nsv(i))
300 IF(igap == 3)THEN
301 IF(gap_s_l_tmp(nsv(i)) /= zero)gap_s_l(i)=
min(gap_s_l(i),gap_s_l_tmp(nsv(i)))
302 gaps_mx =
max(gaps_mx,gap_s(i))
303 gaps_l_mx =
max(gaps_l_mx,gap_s_l(i))
304 ELSE
305 gaps_mx=
max(gaps_mx,gap_s(i))
306 END IF
307 ENDDO
308 ENDIF
309 ENDIF
310
311
312 IF(intth > 0 ) THEN
313 IF(nadmesh == 0)THEN
314 DO i = 1,nsn
315 areas(i) = zero
316 DO j= knod2elc(nsv(i))+1,knod2elc(nsv(i)+1)
317 ie = nod2elc(j)
318 sx1 = x(1,ixc(4,ie)) - x(1,ixc(2,ie))
319 sy1 = x(2,ixc(4,ie)) - x(2,ixc(2,ie))
320 sz1 = x(3,ixc(4,ie)) - x(3,ixc(2,ie))
321 sx2 = x(1,ixc(5,ie)) - x(1,ixc(3,ie))
322 sy2 = x(2,ixc(5,ie)) - x(2,ixc(3,ie))
323 sz2 = x(3,ixc(5,ie)) - x(3,ixc(3,ie))
324 sx3 = sy1*sz2 - sz1*sy2
325 sy3 = sz1*sx2 - sx1*sz2
326 sz3 = sx1*sy2 - sy1*sx2
327 areas(i) = areas(i) + one_over_8*sqrt(sx3*sx3+sy3*sy3+sz3*sz3)
328
329 ielec(i) = ixc(1,ie)
330 END DO
331
332 DO j= knod2eltg(nsv(i))+1,knod2eltg(nsv(i)+1)
333 ie = nod2eltg(j)
334 sx1 = x(1,ixtg(3,ie)) - x(1,ixtg(2,ie))
335 sy1 = x(2,ixtg(3,ie)) - x(2,ixtg(2,ie))
336 sz1 = x(3,ixtg(3,ie)) - x(3,ixtg(2,ie))
337 sx2 = x(1,ixtg(4,ie)) - x(1,ixtg(2,ie))
338 sy2 = x(2,ixtg(4,ie)) - x(2,ixtg(2,ie))
339 sz2 = x(3,ixtg(4,ie)) - x(3,ixtg(2,ie))
340 sx3 = sy1*sz2 - sz1*sy2
341 sy3 = sz1*sx2 - sx1*sz2
342 sz3 = sx1*sy2 - sy1*sx2
343 areas(i) = areas(i) + one_over_6*sqrt(sx3*sx3+sy3*sy3+sz3*sz3)
344
345 ielec(i) = ixtg(1,ie)
346 END DO
347 END DO
348 ELSE
349 DO i = 1,nsn
350 areas(i) = zero
351 DO j=knod2elc(nsv(i))+1,knod2elc(nsv(i)+1)
352 ie = nod2elc(j)
353
354 ip = ipartc(ie)
355 nlev =ipart(10,ip)
356 mylev=sh4tree(3,ie)
357 IF(mylev < 0) mylev=-(mylev+1)
358
359 IF(mylev == nlev)THEN
360 sx1 = x(1,ixc(4,ie)) - x(1,ixc(2,ie))
361 sy1 = x(2,ixc(4,ie)) - x(2,ixc(2,ie))
362 sz1 = x(3,ixc(4,ie)) - x(3,ixc(2,ie))
363 sx2 = x(1,ixc(5,ie)) - x(1,ixc(3,ie))
364 sy2 = x(2,ixc(5,ie)) - x(2,ixc(3,ie))
365 sz2 = x(3,ixc(5,ie)) - x(3,ixc(3,ie))
366 sx3 = sy1*sz2 - sz1*sy2
367 sy3 = sz1*sx2 - sx1*sz2
368 sz3 = sx1*sy2 - sy1*sx2
369 areas(i) = areas(i) + one_over_8*sqrt(sx3*sx3+sy3*sy3+sz3*sz3)
370
371 ielec(i) = ixc(1,ie)
372 END IF
373
374 END DO
375
376 DO j= knod2eltg(nsv(i))+1,knod2eltg(nsv(i)+1)
377 ie = nod2eltg(j)
378 ip = iparttg(ie)
379 nlev =ipart(10,ip)
380 mylev=sh3tree(3,ie)
381 IF(mylev < 0) mylev=-(mylev+1)
382 IF(mylev == nlev)THEN
383 sx1 = x(1,ixtg(3,ie)) - x(1,ixtg(2,ie))
384 sy1 = x(2,ixtg(3,ie)) - x(2,ixtg(2,ie))
385 sz1 = x(3,ixtg(3,ie)) - x(3,ixtg(2,ie))
386 sx2 = x(1,ixtg(4,ie)) - x(1,ixtg(2,ie))
387 sy2 = x(2,ixtg(4,ie)) - x(2,ixtg(2,ie))
388 sz2 = x(3,ixtg(4,ie)) - x(3,ixtg(2,ie))
389 sx3 = sy1*sz2 - sz1*sy2
390 sy3 = sz1*sx2 - sx1*sz2
391 sz3 = sx1*sy2 - sy1*sx2
392 areas(i) = areas(i) + one_over_6*sqrt(sx3*sx3+sy3*sy3+sz3*sz3)
393
394 ielec(i) = ixtg(1,ie)
395 END IF
396
397 END DO
398 END DO
399 END IF
400 END IF
401
402
403
404
405 IF(numels > 0) THEN
406 CALL my_alloc(tagelems,numels)
407 tagelems = 0
408 CALL my_alloc(indexe,numels)
409 indexe = 0
410 ENDIF
411 ninv = 0
412 DO i=1,nrt
413 stf(i)=zero
414 IF(intth > 0 ) ieles(i) = 0
415 IF(slsfac < zero)stf(i)=slsfac
416 gapm=zero
417 inrt=i
418 CALL i4gmx3(x,irect,inrt,gapmx)
419
420 CALL inelts(x ,irect,ixs ,nint,nels ,
421 . inrt ,
area ,noint,0 ,igrsurf%ELTYP,
422 . igrsurf%ELEM)
423
424
425
426
427 IF(nels /= 0)THEN
428 mt=ixs(1,nels)
429 mg=ixs(nixs-1,nels)
430 icontr = igeo(97,mg)
431 IF(mt > 0)THEN
432 DO jj=1,8
433 jjj=ixs(jj+1,nels)
434 xc(jj)=x(1,jjj)
435 yc(jj)=x(2,jjj)
436 zc(jj)=x(3,jjj)
437 END DO
439 IF (icontr==1 ) THEN
440 bulk = pm(107,mt)
441 ELSE
442 bulk = pm(32,mt)
443 END IF
444 stf(i)=slsfac*fillsol(nels)*
area*
area*bulk/vol
445 ELSE
446 IF(nint >= 0) THEN
448 . msgtype=msgwarning,
449 . anmode=aninfo_blind_2,
451 . c1=titr,
452 . i2=ixs(nixs,nels),
453 . c2='SOLID',
454 . i3=i)
455 ENDIF
456 IF(nint < 0) THEN
458 . msgtype=msgwarning,
459 . anmode=aninfo_blind_2,
461 . c1=titr,
462 . i2=ixs(nixs,nels),
463 . c2='SOLID',
464 . i3=i)
465 ENDIF
466 ENDIF
467 IF(igap /= 0 .OR. (nty /=7 .AND. nty /= 20)) gap_m(i)=gapm
468
469 IF(intfric > 0) THEN
470 ip= iparts(nels)
471 ipg = tagprt_fric(ip)
472 IF(ipg > 0) THEN
474 . ipg,intbuf_fric_tab(intfric)%S_TABPARTS_FRIC,
475 . intbuf_fric_tab(intfric)%TABPARTS_FRIC,ipl )
476 ipartfricm(inrt) = ipl
477 ENDIF
478 ENDIF
479
480 cycle
481 ENDIF
482
483
484 CALL ineltc(nelc ,neltg ,inrt ,igrsurf%ELTYP, igrsurf%ELEM)
485
486
487
488 IF(neltg /= 0) THEN
489 mt=ixtg(1,neltg)
490 mg=ixtg(5,neltg)
491 igtyp = igeo(11,mg)
492 igmat = igeo(98,mg)
493 ip = iparttg(neltg)
494 IF (thk_part(ip) /= zero .AND. iintthick == 0) THEN
495 dx=thk_part(ip)*gapscale
496 ELSEIF ( thk(numelc+neltg) /= zero .AND. iintthick == 0)THEN
497 dx=thk(numelc+neltg)*gapscale
498 ELSEIF(igtyp == 17 .OR. igtyp == 51 .OR. igtyp == 52) THEN
499 dx=thk(numelc+neltg)*gapscale
500 ELSE
501 dx=geo(1,mg)*gapscale
502 ENDIF
503 gapm=half*dx
504 gapm_mx=
max(gapm_mx,gapm)
506 dxm=dxm+dx
507 ndx=ndx+1
508 IF(mt > 0)THEN
509 IF(igtyp == 11 .AND. igmat > 0) THEN
510 IF(slsfac < zero)THEN
511 stf(i)=-slsfac
512 ELSEIF ( thk(numelc+neltg) /= zero .AND. iintthick == 0)THEN
513 stf(i)=slsfac*thk(numelc+neltg)*geo(ipgmat + 2 ,mg)
514 ELSE
515 stf(i)=slsfac*geo(1,mg)*geo(ipgmat + 2 ,mg)
516 ENDIF
517 ELSEIF(igtyp == 52 .OR. ((igtyp == 17 .OR. igtyp == 51) .AND. igmat > 0))THEN
518 isubstack = iworksh(3,numelc + neltg)
519 IF(slsfac < zero)THEN
520 stf(i)=-slsfac
521 ELSE
522 stf(i)=slsfac*thk(numelc+neltg)*pm_stack(2 ,isubstack)
523 ENDIF
524 ELSE
525 IF(slsfac < zero)THEN
526 stf(i)=-slsfac
527 ELSEIF ( thk(numelc+neltg) /= zero .AND. iintthick == 0)THEN
528 stf(i)=slsfac*thk(numelc+neltg)*pm(20,mt)
529 ELSE
530 stf(i)=slsfac*geo(1,mg)*pm(20,mt)
531 ENDIF
532 ENDIF
533 ELSE
534 IF(nint >= 0) THEN
536 . msgtype=msgwarning,
537 . anmode=aninfo_blind_2,
539 . c1=titr,
540 . i2=ixtg(nixtg,neltg),
541 . c2='SHELL',
542 . i3=i)
543 END IF
544 IF(nint < 0) THEN
546 . msgtype=msgwarning,
547 . anmode=aninfo_blind_2,
549 . c1=titr,
550 . i2=ixtg(nixtg,neltg),
551 . c2='SHELL',
552 . i3=i)
553 END IF
554 END IF
555 IF(igap /= 0 .OR. (nty /= 7 .AND. nty /= 20)) gap_m(i)=gapm
556
557 IF(intfric > 0) THEN
558 ip= iparttg(neltg)
559 ipg = tagprt_fric(ip)
560 IF(ipg > 0) THEN
562 . ipg,intbuf_fric_tab(intfric)%S_TABPARTS_FRIC,
563 . intbuf_fric_tab(intfric)%TABPARTS_FRIC,ipl )
564 ipartfricm(inrt) = ipl
565 ENDIF
566 ENDIF
567
568 cycle
569 ENDIF
570
571
572
573
574 IF(nelc /= 0) THEN
575 mt=ixc(1,nelc)
576 mg=ixc(6,nelc)
577 ip = ipartc(nelc)
578 igtyp = igeo(11,mg)
579 igmat = igeo(98,mg)
580 IF (thk_part(ip) /= zero .AND. iintthick == 0) THEN
581 dx=thk_part(ip)*gapscale
582 ELSEIF (thk(nelc) /= zero .AND. iintthick == 0) THEN
583 dx=thk(nelc)*gapscale
584 ELSEIF(igtyp == 17 .OR. igtyp == 51 .OR. igtyp == 52) THEN
585 dx=thk(nelc)*gapscale
586 ELSE
587 dx=geo(1,mg)*gapscale
588 ENDIF
589 gapm=half*dx
590 gapm_mx=
max(gapm_mx,gapm)
591 gapmn =
min(gapmn,dx)
592 dxm=dxm+dx
593 ndx=ndx+1
594
595 IF(mt > 0)THEN
596 IF(igtyp == 11 .AND. igmat > 0) THEN
597 IF(slsfac < zero)THEN
598 stf(i)=-slsfac
599 ELSEIF ( thk(nelc) /= zero .AND. iintthick == 0) THEN
600 stf(i)=slsfac*thk(nelc)*geo(ipgmat + 2 ,mg)
601 ELSE
602 stf(i)=slsfac*geo(1,mg)*geo(ipgmat + 2 ,mg)
603 ENDIF
604 ELSEIF(igtyp == 52 .OR. ((igtyp == 17 .OR. igtyp == 51) .AND. igmat > 0))THEN
605 isubstack = iworksh(3,nelc)
606 IF(slsfac < zero)THEN
607 stf(i)=-slsfac
608 ELSE
609 stf(i)=slsfac*thk(nelc)*pm_stack(2 ,isubstack )
610 ENDIF
611 ELSE
612 IF(slsfac < zero)THEN
613 stf(i)=-slsfac
614 ELSEIF (thk(nelc) /= zero .AND. iintthick == 0) THEN
615 stf(i)=slsfac*thk(nelc)*pm(20,mt)
616 ELSE
617 stf(i)=slsfac*geo(1,mg)*pm(20,mt)
618 ENDIF
619 ENDIF
620 ELSE
621 IF(nint >= 0) THEN
623 . msgtype=msgwarning,
624 . anmode=aninfo_blind_2,
626 . c1=titr,
627 . i2=ixc(nixc,nelc),
628 . c2='SHELL',
629 . i3=i)
630 END IF
631 IF(nint < 0) THEN
633 . msgtype=msgwarning,
634 . anmode=aninfo_blind_2,
636 . c1=titr,
637 . i2=ixc(nixc,nelc),
638 . c2='SHELL',
639 . i3=i)
640 END IF
641 END IF
642 IF(igap /=0 .OR. (nty /=7 .AND. nty /= 20)) gap_m(i)=gapm
643
644 IF(intfric > 0) THEN
645 ip= ipartc(nelc)
646 ipg = tagprt_fric(ip)
647 IF(ipg > 0) THEN
649 . ipg,intbuf_fric_tab(intfric)%S_TABPARTS_FRIC,
650 . intbuf_fric_tab(intfric)%TABPARTS_FRIC,ipl )
651 ipartfricm(inrt) = ipl
652 ENDIF
653 ENDIF
654
655 cycle
656 END IF
657
658
659
660
661 print_error = .false.
662 nin25 = 0
663 CALL insol3d(x ,irect ,ixs ,nint ,nels,inrt,
664 .
area ,noint ,knod2els ,nod2els ,0 ,
665 . ixs10 ,ixs16 ,ixs20 ,tagelems,indexe,ninv,ielem,
666 . elem_linked_to_segment ,print_error,
667 . nin25,nty, flag_elem_inter25 )
668 IF(print_error) THEN
669 node_id(1:4) = itab(irect(1:4,inrt))
670
672 . msgtype=msgwarning,
673 . anmode=aninfo_blind_1,
675 . i2=node_id(1),
676 . i3=node_id(2),
677 . i4=node_id(3),
678 . i5=node_id(4),
679 . c1=titr ,
680 . prmod=msg_print)
681 ENDIF
682
683
684
685
686 IF(nels /= 0) THEN
687 mt=ixs(1,nels)
688 IF(intth > 0 ) ieles(i) = nels
689 IF(mt > 0)THEN
690 DO jj=1,8
691 jjj=ixs(jj+1,nels)
692 xc(jj)=x(1,jjj)
693 yc(jj)=x(2,jjj)
694 zc(jj)=x(3,jjj)
695 ENDDO
697 stf(i)=slsfac*fillsol(nels)*
area*
area*pm(32,mt)/vol
698 ELSE
699 IF(nint >= 0) THEN
701 . msgtype=msgwarning,
702 . anmode=aninfo_blind_2,
704 . c1=titr,
705 . i2=ixs(nixs,nels),
706 . c2='SOLID',
707 . i3=i)
708 ENDIF
709 IF(nint < 0) THEN
711 . msgtype=msgwarning,
712 . anmode=aninfo_blind_2,
714 . c1=titr,
715 . i2=ixs(nixs,nels),
716 . c2='SOLID',
717 . i3=i)
718 ENDIF
719 ENDIF
720
721 IF(intfric > 0) THEN
722 ip= iparts(nels)
723 ipg = tagprt_fric(ip)
724 IF(ipg > 0) THEN
726 . ipg , intbuf_fric_tab(intfric)%S_TABPARTS_FRIC,
727 . intbuf_fric_tab(intfric)%TABPARTS_FRIC, ipl )
728 ipartfricm(inrt) = ipl
729 ENDIF
730 ENDIF
731
732 ENDIF
733
734
735
736
737 CALL incoq3(irect ,ixc ,ixtg ,nint ,nelc ,
738 . neltg ,inrt ,geo ,pm ,knod2elc ,
739 . knod2eltg ,nod2elc ,nod2eltg ,thk ,nty ,
740 . igeo ,pm_stack ,iworksh )
741
742
743
744 IF(neltg /= 0) THEN
745 mt=ixtg(1,neltg)
746 mg=ixtg(5,neltg)
747 igtyp = igeo(11,mg)
748 igmat = igeo(98,mg)
749 ip = iparttg(neltg)
750 IF ( thk_part(ip) /= zero .AND. iintthick == 0) THEN
751 dx=thk_part(ip)*gapscale
752 ELSEIF (thk(numelc+neltg) /= zero .AND. iintthick == 0)THEN
753 dx=thk(numelc+neltg)*gapscale
754 ELSEIF(igtyp == 17 .OR. igtyp == 51 .OR. igtyp == 52) THEN
755 dx=thk(numelc+neltg)*gapscale
756 ELSE
757 dx=geo(1,mg)*gapscale
758 ENDIF
759 gapm=half*dx
760 gapm_mx=
max(gapm_mx,gapm)
761 gapmn =
min(gapmn,dx)
762 dxm=dxm+dx
763 ndx=ndx+1
764 IF(mt > 0)THEN
765 IF(igtyp == 11 .AND. igmat > 0) THEN
766 IF(slsfac < zero)THEN
767 stf(i)=-slsfac
768 ELSEIF ( thk(numelc+neltg) /= zero .AND. iintthick == 0) THEN
769 stf(i)=slsfac*thk(numelc+neltg)*geo(ipgmat + 2 ,mg)
770 ELSE
771 stf(i)=slsfac*geo(1,mg)*geo(ipgmat + 2 ,mg)
772 ENDIF
773 ELSEIF(igtyp == 52 .OR. ((igtyp == 17 .OR. igtyp == 51) .AND. igmat > 0))THEN
774 isubstack = iworksh(3,numelc+neltg)
775 IF(slsfac < zero)THEN
776 stf(i)=-slsfac
777 ELSE
778 stf(i)=slsfac*thk(numelc+neltg)*pm_stack(2 ,isubstack)
779 ENDIF
780 ELSE
781 IF(slsfac < zero)THEN
782 stf(i)=-slsfac
783 ELSEIF ( thk(numelc+neltg) /= zero .AND. iintthick == 0) THEN
784 stf(i)=slsfac*thk(numelc+neltg)*pm(20,mt)
785 ELSE
786 stf(i)=slsfac*geo(1,mg)*pm(20,mt)
787 ENDIF
788 ENDIF
789 ELSE
790 IF(nint >= 0) THEN
792 . msgtype=msgwarning,
793 . anmode=aninfo_blind_2,
795 . c1=titr,
796 . i2=ixtg(nixtg,neltg),
797 . c2='SHELL',
798 . i3=i)
799 ENDIF
800 IF(nint < 0) THEN
802 . msgtype=msgwarning,
803 . anmode=aninfo_blind_2,
805 . c1=titr,
806 . i2=ixtg(nixtg,neltg),
807 . c2='SHELL',
808 . i3=i)
809 ENDIF
810 ENDIF
811
812 IF(intfric > 0) THEN
813 ip= iparttg(neltg)
814 ipg = tagprt_fric(ip)
815 IF(ip > 0) THEN
817 . ipg,intbuf_fric_tab(intfric)%S_TABPARTS_FRIC,
818 . intbuf_fric_tab(intfric)%TABPARTS_FRIC,ipl )
819 ipartfricm(inrt) = ipl
820 ENDIF
821 ENDIF
822
823
824
825
826
827 ELSEIF(nelc /= 0) THEN
828 mt=ixc(1,nelc)
829 mg=ixc(6,nelc)
830 igtyp = igeo(11,mg)
831 igmat = igeo(98,mg)
832 ip = ipartc(nelc)
833 IF ( thk_part(ip) /= zero .AND. iintthick == 0) THEN
834 dx=thk_part(ip)*gapscale
835 ELSEIF ( thk(nelc) /= zero .AND. iintthick == 0) THEN
836 dx=thk(nelc)*gapscale
837 ELSEIF(igtyp == 17 .OR. igtyp == 51 .OR. igtyp == 52) THEN
838 dx=thk(nelc)*gapscale
839 ELSE
840 dx=geo(1,mg)*gapscale
841 ENDIF
842 gapm=half*dx
843 gapm_mx=
max(gapm_mx,gapm)
844 gapmn =
min(gapmn,dx)
845 dxm=dxm+dx
846 ndx=ndx+1
847 IF(mt > 0)THEN
848 IF(igtyp == 11 .AND. igmat > 0) THEN
849 IF(slsfac < zero)THEN
850 stf(i)=-slsfac
851 ELSEIF (thk(nelc) /= zero .AND. iintthick == 0) THEN
852 stf(i)=slsfac*thk(nelc)*geo(ipgmat + 2 ,mg)
853 ELSE
854 stf(i)=slsfac*geo(1,mg)*geo(ipgmat + 2 ,mg)
855 ENDIF
856 ELSEIF(igtyp == 52 .OR. ((igtyp == 17 .OR. igtyp == 51) .AND. igmat > 0))THEN
857 isubstack = iworksh(3,nelc)
858 IF(slsfac < zero)THEN
859 stf(i)=-slsfac
860 ELSE
861 stf(i)=slsfac*thk(nelc)*pm_stack(2 ,isubstack)
862 ENDIF
863 ELSE
864 IF(slsfac < zero)THEN
865 stf(i)=-slsfac
866 ELSEIF ( thk(nelc) /= zero .AND. iintthick == 0) THEN
867 stf(i)=slsfac*thk(nelc)*pm(20,mt)
868 ELSE
869 stf(i)=slsfac*geo(1,mg)*pm(20,mt)
870 ENDIF
871 ENDIF
872 ELSE
873 IF(nint >= 0) THEN
875 . msgtype=msgwarning,
876 . anmode=aninfo_blind_2,
878 . c1=titr,
879 . i2=ixc(nixc,nelc),
880 . c2='SHELL',
881 . i3=i)
882 ENDIF
883 IF(nint < 0) THEN
885 . msgtype=msgwarning,
886 . anmode=aninfo_blind_2,
888 . c1=titr,
889 . i2=ixc(nixc,nelc),
890 . c2='SHELL',
891 . i3=i)
892 ENDIF
893 ENDIF
894
895 IF(intfric > 0) THEN
896 ip= ipartc(nelc)
897 ipg = tagprt_fric(ip)
898 IF(ipg > 0) THEN
900 . ipg,intbuf_fric_tab(intfric)%S_TABPARTS_FRIC,
901 . intbuf_fric_tab(intfric)%TABPARTS_FRIC,ipl )
902 ipartfricm(inrt) = ipl
903 ENDIF
904 ENDIF
905
906 ENDIF
907
908 IF(igap /= 0 .OR. (nty /= 7 .AND. nty /= 20)) gap_m(i)=gapm
909
910
911
912
913
914 IF(nels+nelc+neltg == 0)THEN
915
916
917 IF(nint > 0) THEN
919 . msgtype=msgerror,
920 . anmode=aninfo_blind_2,
922 . c1=titr,
923 . i2=i)
924 ENDIF
925 IF(nint < 0) THEN
927 . msgtype=msgerror,
928 . anmode=aninfo_blind_2,
930 . c1=titr,
931 . i2=i)
932 ENDIF
933
934 ENDIF
935
936 enddo
937
938 IF(numels > 0) DEALLOCATE(tagelems,indexe)
939
941 . msgtype=msgwarning,
942 . anmode=aninfo_blind_1,
944 . c1=titr,
945 . prmod=msg_print)
947 . msgtype=msgwarning,
948 . anmode=aninfo_blind_1,
950 . c1=titr,
951 . prmod=msg_print)
952 IF(ninv > 0 .AND.nint>0)
954 . msgtype=msgwarning,
955 . anmode=aninfo_blind_1,
957 . c1=titr,
958 . i2=ninv)
959
960 IF(ninv > 0 .AND.nint< 0)
962 . msgtype=msgwarning,
963 . anmode=aninfo_blind_1,
965 . c1=titr,
966 . i2=ninv)
967
968
969
970 DO i=nrt+1,nrt+nrt_ige
971 stf(i)=zero
972 IF(intth > 0) ieles(i) = 0
973 IF(slsfac < zero)stf(i)=slsfac
974 gapm =zero
975 inrt=i
976 CALL i4gmx3(x,irect,inrt,gapmx)
977
978
979
980 CALL ineltigeo(x ,irect ,ixs ,nint ,nelig3d ,
981 . inrt ,
area ,noint ,0 ,igrsurf%ELTYP_IGE,
982 . ixig3d ,kxig3d ,igeo ,igrsurf%ELEM_IGE)
983 IF(nelig3d /= 0)THEN
984 mt=kxig3d(1,nelig3d)
985 IF(mt > 0)THEN
986 ipid = kxig3d(2,nelig3d)
987 px = igeo(41,ipid)-1
988 py = igeo(42,ipid)-1
989 pz = igeo(43,ipid)-1
990 coin_ige(1) = (px+1)*py+1
991 coin_ige(2) = (px+1)*(py+1)
992 coin_ige(3) = px+1
993 coin_ige(4) = 1
994 coin_ige(5) = (px+1)*(py+1)*pz+(px+1)*py+1
995 coin_ige(6) = (px+1)*(py+1)*(pz+1)
996 coin_ige(7) = (px+1)*(py+1)*pz+px+1
997 coin_ige(8) = (px+1)*(py+1)*pz+1
998 DO jj=1,8
999 xc(jj)=x(1,ixig3d(kxig3d(4,nelig3d)+coin_ige(jj)-1))
1000 yc(jj)=x(2,ixig3d(kxig3d(4,nelig3d)+coin_ige(jj)-1))
1001 zc(jj)=x(3,ixig3d(kxig3d(4,nelig3d)+coin_ige(jj)-1))
1002 END DO
1004 stf(i)=slsfac*
area*
area*pm(32,mt)/vol
1005 stf(i)=stf(i)*((px+1)*(py+1)+(py+1)*(pz+1)+(pz+1)*(px+1))/3
1006 ELSE
1007 IF(nint >= 0) THEN
1009 . msgtype=msgwarning,
1010 . anmode=aninfo_blind_2,
1012 . c1=titr,
1013 . i2=kxig3d(5,nelig3d),
1014 . c2='ISOGEOMETRIC SOLID',
1015 . i3=i)
1016 ENDIF
1017 IF(nint < 0) THEN
1019 . msgtype=msgwarning,
1020 . anmode=aninfo_blind_2,
1022 . c1=titr,
1023 . i2=kxig3d(5,nelig3d),
1024 . c2='ISOGEOMETRIC SOLID',
1025 . i3=i)
1026 ENDIF
1027 ENDIF
1028 ELSEIF(nelig3d == 0)THEN
1029
1030
1031 IF(nint > 0) THEN
1033 . msgtype=msgerror,
1034 . anmode=aninfo_blind_2,
1036 . c1=titr,
1037 . i2=i)
1038 ENDIF
1039 IF(nint < 0) THEN
1041 . msgtype=msgerror,
1042 . anmode=aninfo_blind_2,
1044 . c1=titr,
1045 . i2=i)
1046 ENDIF
1047
1048 ENDIF
1049 enddo
1050
1051
1052
1053
1054 gapmx=sqrt(gapmx)
1055 IF(igap == 0)THEN
1056
1057 IF(gap <= zero)THEN
1058 IF(ndx /= 0)THEN
1059 gap = dxm/ndx
1060 gap =
min(half*gapmx,gap)
1061 ELSE
1062 gap = em01 * gapmx
1063 ENDIF
1064 IF (it19 <= 0 .AND. .NOT.type18) WRITE(iout,1300)gap
1065 ENDIF
1066 gapmin = gap
1067
1068 IF (gapmin <= 0) THEN
1070 . msgtype=msgerror,
1071 . anmode=aninfo,
1073 . c1=titr)
1074 ENDIF
1075 IF ((inacti /= 7).AND.(gap > 0.5*gapmx) .AND. (irem_gap /= 2)) THEN
1076 gaptmp = half*gapmx
1078 . msgtype=msgwarning,
1079 . anmode=aninfo_blind_2,
1081 . c1=titr,
1082 . r1=gap,
1083 . r2=gaptmp)
1084 ENDIF
1085 ELSE
1086
1087
1088
1089 IF(gap <= zero)THEN
1090 IF(ndx /= 0)THEN
1091 gapmin = gapmn
1092 gapmin =
min(half*gapmx,gapmin)
1093 ELSE
1094 gapmin = em01 * gapmx
1095 ENDIF
1096 IF (gapmin <= 0) THEN
1098 . msgtype=msgerror,
1099 . anmode=aninfo,
1101 . c1=titr)
1102 ENDIF
1103 IF (it19 <= 0 .AND. .NOT.type18) WRITE(iout,1300)gapmin
1104 ELSE
1105 gapmin = gap
1106 ENDIF
1107
1108 IF(igap == 3) THEN
1109 gap =
max(
min(gaps_mx+gapm_mx,gaps_l_mx+gapm_l_mx) ,gapmin)
1110 ELSE
1111 gap =
max(gaps_mx+gapm_mx,gapmin)
1112 ENDIF
1114 IF ((igap /= 3).AND.(irem_gap /= 2)) THEN
1115 IF(inacti /= 7.AND.gap > half*gapmx .AND. iddlevel == 1)THEN
1116 gaptmp = 0.5*gapmx
1118 . msgtype=msgwarning,
1119 . anmode=aninfo_blind_2,
1121 . c1=titr,
1122 . r1=gap)
1123 ENDIF
1124 ENDIF
1125 ENDIF
1126
1127 IF(intth /= 0)THEN
1128 IF(drad == zero)THEN
1129
1131 ELSEIF(drad < gap)THEN
1132
1133 drad=gap
1134 END IF
1135 WRITE(iout,2001)drad
1136
1137
1138 IF(drad > gapmx)THEN
1140 . msgtype=msgwarning,
1141 . anmode=aninfo_blind_2,
1143 . c1=titr,
1144 . r1=drad ,
1145 . r2=gapmx,
1147 END IF
1148 END IF
1149
1150
1151 IF(intfric > 0) THEN
1152 IF(numels /= 0)THEN
1153 DO i = 1,nsn
1154 ipfmax = 0
1155 ipflmax = 0
1156 DO j= knod2els(nsv(i))+1,knod2els(nsv(i)+1)
1157 ie = nod2els(j)
1158 ip = iparts(ie)
1159 ipg = tagprt_fric(ip)
1160 IF(ipg > 0 .AND. ip > ipfmax) THEN
1162 . ipg,intbuf_fric_tab(intfric)%S_TABPARTS_FRIC,
1163 . intbuf_fric_tab(intfric)%TABPARTS_FRIC,ipl )
1164 IF(ipl /= 0) THEN
1165 ipfmax = ip
1166 ipflmax = ipl
1167 ENDIF
1168 ENDIF
1169 ENDDO
1170 IF(ipfmax /= 0) THEN
1171 ipartfrics(i) = ipflmax
1172 ENDIF
1173
1174 ENDDO
1175 ENDIF
1176
1177 IF(numelc /= 0 .OR. numeltg /= 0) THEN
1178 DO i = 1,nsn
1179 ipfmax = 0
1180 ipflmax = 0
1181 DO j= knod2elc(nsv(i))+1,knod2elc(nsv(i)+1)
1182 ie = nod2elc(j)
1183 ip = ipartc(ie)
1184 ipg = tagprt_fric(ip)
1185 IF(ipg > 0 .AND. ip > ipfmax) THEN
1187 . ipg,intbuf_fric_tab(intfric)%S_TABPARTS_FRIC,
1188 . intbuf_fric_tab(intfric)%TABPARTS_FRIC,ipl )
1189 IF(ipl /= 0) THEN
1190 ipfmax = ip
1191 ipflmax = ipl
1192 ENDIF
1193 ENDIF
1194 ENDDO
1195
1196 DO j= knod2eltg(nsv(i))+1,knod2eltg(nsv(i)+1)
1197 ie = nod2eltg(j)
1198 ip = iparttg(ie)
1199 ipg = tagprt_fric(ip)
1200 IF(ipg > 0.AND.ip > ipfmax) THEN
1202 . ipg,intbuf_fric_tab(intfric)%S_TABPARTS_FRIC,
1203 . intbuf_fric_tab(intfric)%TABPARTS_FRIC,ipl
1204
1205 IF(ipl /= 0) THEN
1206 ipfmax = ip
1207 ipflmax = ipl
1208 ENDIF
1209 ENDIF
1210 ENDDO
1211 IF(ipfmax /= 0) THEN
1212 ipartfrics(i) = ipflmax
1213 ENDIF
1214
1215 ENDDO
1216 ENDIF
1217 ENDIF
1218
1219
1220
1221
1222
1223 DO l=1,nsn
1224 stfn(l) = 1.
1225 ENDDO
1226
1227
1228
1229 bgapsmx = zero
1230 IF (igap == 0) THEN
1231 gapinf=gap
1232 ELSE
1233 gapinf = ep30
1234 gapinf_s = ep30
1235 gapinf_m = ep30
1236 DO i = 1, nsn
1237 gapinf_s =
min(gapinf_s,gap_s(i))
1238 bgapsmx =
max(bgapsmx,gap_s(i))
1239 ENDDO
1240 DO i = 1, nrt+nrt_ige
1241 gapinf_m =
min(gapinf_m,gap_m(i))
1242 ENDDO
1243 gapinf=gapinf_s+gapinf_m
1244 gapinf=
max(gapinf,gapmin)
1245 ENDIF
1246 DEALLOCATE( gap_s_l_tmp )
1247 RETURN
1248 1300 FORMAT(2x,'GAP MIN = ',1pg20.13)
1249 2001 FORMAT(2x,'Maximum distance for radiation computation = ',
1250 . 1pg20.13)
subroutine area(d1, x, x2, y, y2, eint, stif0)
subroutine i4gmx3(x, irect, i, gapmax)
subroutine friction_parts_search(ip, npartsfric, partsfric, ipl)
subroutine incoq3(irect, ixc, ixtg, nint, nel, neltg, is, geo, pm, knod2elc, knod2eltg, nod2elc, nod2eltg, thk, nty, igeo, pm_stack, iworksh)
subroutine inelts(x, irect, ixs, nint, nel, i, area, noint, ir, surf_eltyp, surf_elem)
subroutine ineltigeo(xe, irect, ixs, nint, nelig3d, i, area, noint, ir, surf_eltyp_ige, ixig3d, kxig3d, igeo, surf_elem_ige)
subroutine ineltc(nelc, neltg, is, surf_eltyp, surf_elem)
subroutine insol3d(x, irect, ixs, nint, nel, i, area, noint, knod2els, nod2els, ir, ixs10, ixs16, ixs20, tagelems, indexe, ninv, ielem_m, elem_linked_to_segment, print_error, nin25, nty, flag_elem_inter25)
integer, parameter nchartitle
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)