49
50
51
55 USE intbuf_fric_mod
56 use element_mod , only :nixs,nixc,nixtg
57
58
59
60#include "implicit_f.inc"
61
62
63
64#include "com04_c.inc"
65#include "param_c.inc"
66#include "scr17_c.inc"
67#include "remesh_c.inc"
68
69
70
71 INTEGER NRTS, NRTM, NINT, NTY, NOINT, NSN, NMN
72 INTEGER , INTENT(IN) :: INTFRIC
73 INTEGER IRECTS(4,*), IRECTM(4,*), IXS(NIXS,*), IXC(NIXC,*),
74 . (*), IXTG(NIXTG,*),
75 . KNOD2ELS(*), KNOD2ELC(*), KNOD2ELTG(*), NOD2ELS(*), NOD2ELC(*),
76 . NOD2ELTG(*),
77 . INTTH, IELES(*), MSR(*), IELEM21(*), IXS10(*),
78 . IXS16(*), IXS20(*),IGEO(*),SH3TREE(KSH3TREE,*), SH4TREE(KSH4TREE,*),
79 . IPART(LIPART1,*),IPARTC(*),IPARTTG(*),IWORKSH(*)
80
82 . x(3,*), pm(npropm,*), geo(npropg,*), areas(*),thk(*),
83 . as(*), bs(*),pm_stack(*)
84 INTEGER ID
85 CHARACTER(LEN=NCHARTITLE) :: TITR
86 TYPE (SURF_) :: IGRSURFS
87 TYPE (SURF_) :: IGRSURFM
88 TYPE(INTBUF_FRIC_STRUCT_),INTENT(INOUT) :: INTBUF_FRIC_TAB(NINTERFRIC)
89 INTEGER,INTENT(INOUT) :: IPARTFRICS(),IPARTFRICM(NRTM)
90 INTEGER,INTENT(IN) :: TAGPRT_FRIC(NPART)
91 INTEGER, DIMENSION(NUMELS), INTENT(IN) :: IPARTS
92
93
94
95 INTEGER I, J, INRT, NELS, NELC, NELTG, IE, II, MAT,N1,N2,N3,N4
96 INTEGER ITMP(NUMNOD),NLEV, MYLEV,IP,NELEM,STAT,IPG,IPL,IPFMAX,IPFLMAX
97 INTEGER, DIMENSION(:),ALLOCATABLE ::INRTIE
98
100 .
area,sx1,sy1,sz1,sx2,sy2,sz2,sx3,sy3,sz3
101
102 INTEGER :: NB_CONTRIB
103 INTEGER, DIMENSION(:), ALLOCATABLE :: CONTRIB_KEY, CONTRIB_VALUE
104
105
106 IF(intth > 0)THEN
107 nelem = numelc+numeltg+numels+numelr
108 + + numelp+numelt+numelq+numelx+numelig3d
109 ALLOCATE(inrtie(nelem),stat=stat)
110 IF (stat /= 0)
CALL ancmsg(msgid=268,anmode=aninfo,
111 . msgtype=msgerror,
112 . c1='INRTIE')
113 inrtie=0
114 ALLOCATE(contrib_key(nelem),contrib_value(nelem))
115 END IF
116
117
118
119
120 DO 250 i=1,nrts
121 inrt=i
122
123 CALL inelts(x ,irects,ixs ,nint,nels ,
124 . inrt ,
area ,noint,0 ,igrsurfs%ELTYP,
125 . igrsurfs%ELEM)
126 IF(nels/=0)THEN
127 ieles(i)=nels
128 IF(intth > 0) inrtie(nels) = inrt
129 GO TO 250
130 ELSE
131 CALL ineltc(nelc ,neltg ,inrt ,igrsurfs%ELTYP, igrsurfs%ELEM)
132 IF(neltg/=0) THEN
133 ieles(i)=neltg+numels+numelc
134 GO TO 250
135 ELSEIF(nelc/=0) THEN
136 ieles(i)=nelc+numels
137 GO TO 250
138 END IF
139 END IF
140
141
142
143 CALL insol3(x,irects,ixs,nint,nels,inrt,
144 .
area,noint,knod2els ,nod2els ,0 ,ixs10,
145 . ixs16,ixs20)
146 IF(nels/=0) ieles(i)=nels
147
148
149
150 CALL incoq3(irects,ixc ,ixtg ,nint ,nelc ,
151 . neltg,inrt,geo ,pm ,knod2elc ,
152 . knod2eltg ,nod2elc ,nod2eltg,thk,nty,igeo ,
153 . pm_stack , iworksh )
154 IF(neltg/=0) THEN
155 ieles(i)=neltg
156 ELSEIF(nelc/=0) THEN
157 ieles(i)=nelc
158 ENDIF
159
160 IF(nels+nelc+neltg==0)THEN
161 IF(nint>0) THEN
163 . msgtype=msgerror,
164 . anmode=aninfo_blind_2,
166 . c1=titr,
167 . i2=i)
168 ENDIF
169 IF(nint<0) THEN
171 . msgtype=msgerror,
172 . anmode=aninfo_blind_2,
174 . c1=titr,
175 . i2=i)
176 ENDIF
177 ENDIF
178 250 CONTINUE
179 DO 500 i=1,nrtm
180 inrt=i
181
182 CALL inelts(x ,irectm,ixs ,nint,nels ,
183 . inrt ,
area ,noint,0 ,igrsurfm%ELTYP,
184 . igrsurfm%ELEM)
185 IF(nels/=0)THEN
186 ielem21(nels)=1
187
188 IF(intfric > 0) THEN
189 ip= iparts(nels)
190 ipg = tagprt_fric(ip)
191 IF(ipg > 0) THEN
193 . ipg,intbuf_fric_tab(intfric)%S_TABPARTS_FRIC,
194 . intbuf_fric_tab(intfric)%TABPARTS_FRIC,ipl )
195 ipartfricm(i) = ipl
196 ENDIF
197 ENDIF
198 GO TO 500
199 ELSE
200 CALL ineltc(nelc ,neltg ,inrt ,igrsurfm%ELTYP, igrsurfm%ELEM)
201 IF(neltg/=0) THEN
202 ielem21(numels+numelq+numelc+numelt
203 . +numelp+numelr+neltg)=1
204
205 IF(intfric > 0) THEN
206 ip= iparttg(neltg)
207 ipg = tagprt_fric(ip)
208 IF(ipg > 0) THEN
210 . ipg,intbuf_fric_tab(intfric)%S_TABPARTS_FRIC,
211 . intbuf_fric_tab(intfric)%TABPARTS_FRIC,ipl )
212 ipartfricm(i) = ipl
213 ENDIF
214 ENDIF
215 GO TO 500
216 ELSEIF(nelc/=0) THEN
217 ielem21(numels+numelq+nelc)=1
218
219 IF(intfric > 0) THEN
220 ip= ipartc(nelc)
221 ipg = tagprt_fric(ip)
222 IF(ipg > 0) THEN
224 . ipg,intbuf_fric_tab(intfric)%S_TABPARTS_FRIC,
225 . intbuf_fric_tab(intfric)%TABPARTS_FRIC,ipl )
226 ipartfricm(i) = ipl
227 ENDIF
228 ENDIF
229 GO TO 500
230 END IF
231 END IF
232
233
234
235 CALL insol3(x,irectm,ixs,nint,nels,inrt,
236 .
area,noint,knod2els ,nod2els ,0 ,ixs10,
237 . ixs16,ixs20)
238
239
240
241 CALL incoq3(irectm,ixc ,ixtg ,nint ,nelc ,
242 . neltg,inrt,geo ,pm ,knod2elc ,
243 . knod2eltg ,nod2elc ,nod2eltg,thk,nty,igeo,
244 . pm_stack , iworksh )
245
246 IF(nels+nelc+neltg==0)THEN
247
248 IF(nint>0) THEN
250 . msgtype=msgerror,
251 . anmode=aninfo_blind_2,
253 . c1=titr,
254 . i2=i)
255 ENDIF
256 IF(nint<0) THEN
258 . msgtype=msgerror,
259 . anmode=aninfo_blind_2,
261 . c1=titr,
262 . i2=i)
263 ENDIF
264
265 ELSE
266 IF(nels/=0) THEN
267 ielem21(nels)=1
268
269 IF(intfric > 0) THEN
270 ip= iparts(nels)
271 ipg = tagprt_fric(ip)
272 IF(ipg > 0) THEN
274 . ipg,intbuf_fric_tab(intfric)%S_TABPARTS_FRIC,
275 . intbuf_fric_tab(intfric)%TABPARTS_FRIC,ipl )
276 ipartfricm(i) = ipl
277 ENDIF
278 ENDIF
279 ENDIF
280 IF(neltg/=0) THEN
281 ielem21(numels+numelq+numelc+numelt
282 . +numelp+numelr+neltg)=1
283
284 IF(intfric > 0) THEN
285 ip= iparttg(neltg)
286 ipg = tagprt_fric(ip)
287 IF(ipg > 0) THEN
289 . ipg,intbuf_fric_tab(intfric)%S_TABPARTS_FRIC,
290 . intbuf_fric_tab(intfric)%TABPARTS_FRIC,ipl )
291 ipartfricm(i) = ipl
292 ENDIF
293 ENDIF
294 ENDIF
295 IF(nelc/=0) THEN
296 ielem21(numels+numelq+nelc)=1
297
298 IF(intfric > 0) THEN
299 ip= ipartc(nelc)
300 ipg = tagprt_fric(ip)
301 IF(ipg > 0) THEN
303 . ipg,intbuf_fric_tab(intfric)%S_TABPARTS_FRIC,
304 . intbuf_fric_tab(intfric)%TABPARTS_FRIC,ipl )
305 ipartfricm(i) = ipl
306 ENDIF
307 ENDIF
308 ENDIF
309 END IF
310 500 CONTINUE
311
312
313
314 itmp=0
315 DO i=1,nmn
316 ii=msr(i)
317 itmp(ii)=i
318 END DO
319
320 DO i=1,nrtm
321 DO j=1,4
322 irectm(j,i)=itmp(irectm(j,i))
323 END DO
324 END DO
325
326
327
328 IF(intth > 0 ) THEN
329
330 IF(numelc/=0) THEN
331
332 IF(nadmesh==0)THEN
333 DO i = 1,nsn
334 areas(i) = zero
335
336
337
338 nb_contrib = 0
339 DO j= knod2elc(nsv(i))+1,knod2elc(nsv(i)+1)
340 nb_contrib = nb_contrib + 1
341 ie = nod2elc(j)
342 contrib_key(nb_contrib) = ixc(nixc,ie)
343 contrib_value(nb_contrib) = ie
344 ENDDO
346
347 DO j = 1,nb_contrib
348 ie = contrib_value(j)
349 sx1 = x(1,ixc(4,ie)) - x(1,ixc(2,ie))
350 sy1 = x(2,ixc(4,ie)) - x(2,ixc(2,ie))
351 sz1 = x(3,ixc(4,ie)) - x(3,ixc(2,ie))
352 sx2 = x(1,ixc(5,ie)) - x(1,ixc(3,ie))
353 sy2 = x(2,ixc(5,ie)) - x(2,ixc(3,ie))
354 sz2 = x(3,ixc(5,ie)) - x(3,ixc(3,ie))
355 sx3 = sy1*sz2 - sz1*sy2
356 sy3 = sz1*sx2 - sx1*sz2
357 sz3 = sx1*sy2 - sy1*sx2
358 area = one_over_8*sqrt(sx3*sx3+sy3*sy3+sz3*sz3)
359 areas(i) = areas(i) +
area
360
361 mat = ixc(1,ie)
362 as(i)= as(i)+pm(75,mat)*
area
363 bs(i)= bs(i)+pm(76,mat)*
area
364 END DO
365
366
367
368 nb_contrib = 0
369 DO j= knod2eltg(nsv(i))+1,knod2eltg(nsv(i)+1)
370 nb_contrib = nb_contrib + 1
371 ie = nod2eltg(j)
372 contrib_key(nb_contrib) = ixtg(nixtg,ie)
373 contrib_value(nb_contrib) = ie
374 ENDDO
376
377
378 DO j = 1,nb_contrib
379 ie = contrib_value(j)
380 sx1 = x(1,ixtg(3,ie)) - x(1,ixtg(2,ie))
381 sy1 = x(2,ixtg(3,ie)) - x(2,ixtg(2,ie))
382 sz1 = x(3,ixtg(3,ie)) - x(3,ixtg(2,ie))
383 sx2 = x(1,ixtg(4,ie)) - x(1,ixtg(2,ie))
384 sy2 = x(2,ixtg(4,ie)) - x(2,ixtg(2,ie))
385 sz2 = x(3,ixtg(4,ie)) - x(3,ixtg(2,ie))
386 sx3 = sy1*sz2 - sz1*sy2
387 sy3 = sz1*sx2 - sx1*sz2
388 sz3 = sx1*sy2 - sy1*sx2
389 area = one_over_6*sqrt(sx3*sx3+sy3*sy3+sz3*sz3)
390 areas(i) = areas(i)+
area
391
392 mat = ixtg(1,ie)
393 as(i)= as(i)+pm(75,mat)*
area
394 bs(i)= bs(i)+pm(76,mat)*
area
395 END DO
396 as(i)=as(i)/
max(em20,areas(i))
397 bs(i)=bs(i)/
max(em20,areas(i))
398 END DO
399 ELSE
400 DO i = 1,nsn
401 areas(i) = zero
402
403
404
405 nb_contrib = 0
406 DO j= knod2elc(nsv(i))+1,knod2elc(nsv(i)+1)
407 nb_contrib = nb_contrib + 1
408 ie = nod2elc(j)
409 contrib_key(nb_contrib) = ixc(nixc,ie)
410 contrib_value(nb_contrib) = ie
411 ENDDO
413
414 DO j = 1,nb_contrib
415 ie = contrib_value(j)
416
417 ip = ipartc(ie)
418 nlev =ipart(10,ip)
419 mylev=sh4tree(3,ie)
420 IF(mylev < 0) mylev=-(mylev+1)
421
422 IF(mylev==nlev)THEN
423 sx1 = x(1,ixc(4,ie)) - x(1,ixc(2,ie))
424 sy1 = x(2,ixc(4,ie)) - x(2,ixc(2,ie))
425 sz1 = x(3,ixc(4,ie)) - x(3,ixc(2,ie))
426 sx2 = x(1,ixc(5,ie)) - x(1,ixc(3,ie))
427 sy2 = x(2,ixc(5,ie)) - x(2,ixc(3,ie))
428 sz2 = x(3,ixc(5,ie)) - x(3,ixc(3,ie))
429 sx3 = sy1*sz2 - sz1*sy2
430 sy3 = sz1*sx2 - sx1*sz2
431 sz3 = sx1*sy2 - sy1*sx2
432 area = one_over_8*sqrt(sx3*sx3+sy3*sy3+sz3*sz3)
433 areas(i) = areas(i) +
area
434
435 mat = ixc(1,ie)
436 as(i)= as(i)+pm(75,mat)*
area
437 bs(i)= bs(i)+pm(76,mat)*
area
438 ENDIF
439 END DO
440
441
442
443 nb_contrib = 0
444 DO j= knod2eltg(nsv(i))+1,knod2eltg(nsv(i)+1)
445 nb_contrib = nb_contrib + 1
446 ie = nod2eltg(j)
447 contrib_key(nb_contrib) = ixtg(nixtg,ie)
448 contrib_value(nb_contrib) = ie
449 ENDDO
451 DO j = 1,nb_contrib
452 ie = contrib_value(j)
453 ip = iparttg(ie)
454 nlev =ipart(10,ip)
455 mylev=sh3tree(3,ie)
456 IF(mylev < 0) mylev=-(mylev+1)
457
458 IF(mylev==nlev)THEN
459 sx1 = x(1,ixtg(3,ie)) - x(1,ixtg(2,ie))
460 sy1 = x(2,ixtg(3,ie)) - x(2,ixtg(2,ie))
461 sz1 = x(3,ixtg(3,ie)) - x(3,ixtg(2,ie))
462 sx2 = x(1,ixtg(4,ie)) - x(1,ixtg(2,ie))
463 sy2 = x(2,ixtg(4,ie)) - x(2,ixtg(2,ie))
464 sz2 = x(3,ixtg(4,ie)) - x(3,ixtg(2,ie))
465 sx3 = sy1*sz2 - sz1*sy2
466 sy3 = sz1*sx2 - sx1*sz2
467 sz3 = sx1*sy2 - sy1*sx2
468 area = one_over_6*sqrt(sx3*sx3+sy3*sy3+sz3*sz3)
469 areas(i) = areas(i)+
area
470
471 mat = ixtg(1,ie)
472 as(i)= as(i)+pm(75,mat)*
area
473 bs(i)= bs(i)+pm(76,mat)*
area
474 END IF
475
476 END DO
477 as(i)=as(i)/
max(em20,areas(i))
478 bs(i)=bs(i)/
max(em20,areas(i))
479 END DO
480 END IF
481 ENDIF
482
483 IF(numels/=0)THEN
484 DO i = 1,nsn
485 areas(i) = zero
486
487
488 nb_contrib = 0
489 DO j= knod2els(nsv(i))+1,knod2els(nsv(i)+1)
490 nb_contrib = nb_contrib + 1
491 ie = nod2els(j)
492 contrib_key(nb_contrib) = ixs(nixs,ie)
493 contrib_value(nb_contrib) = ie
494 ENDDO
496 DO j = 1,nb_contrib
497 ie = contrib_value(j)
498 inrt = inrtie(ie)
499 IF(inrt/=0)THEN
500 n1=irects(1,inrt)
501 n2=irects(2,inrt)
502 n3=irects(3,inrt)
503 n4=irects(4,inrt)
504 sx1 = x(1,n3) - x(1,n1)
505 sy1 = x(2,n3) - x(2,n1)
506 sz1 = x(3,n3) - x(3,n1)
507 sx2 = x(1,n4) - x(1,n2)
508 sy2 = x(2,n4) - x(2,n2)
509 sz2 = x(3,n4) - x(3,n2)
510 sx3 = sy1*sz2 - sz1*sy2
511 sy3 = sz1*sx2 - sx1*sz2
512 sz3 = sx1*sy2 - sy1*sx2
513 area = one_over_8*sqrt(sx3*sx3+sy3*sy3+sz3*sz3)
514 areas(i) = areas(i) +
area
515
516 mat = ixs(1,ie)
517 as(i)= as(i)+pm(75,mat)*
area
518 bs(i)= bs(i)+pm(76,mat)*
area
519 ENDIF
520 END DO
521 as(i)=as(i)/
max(em20,areas(i))
522 bs(i)=bs(i)/
max(em20,areas(i))
523 ENDDO
524 ENDIF
525 ENDIF
526
527 IF(intth > 0) THEN
528 DEALLOCATE(contrib_key,contrib_value)
529 ENDIF
530
531
532 IF(intfric > 0) THEN
533 IF(numels /= 0)THEN
534 DO i = 1,nsn
535 ipfmax = 0
536 ipflmax = 0
537 DO j= knod2els(nsv(i))+1,knod2els(nsv(i)+1)
538 ie = nod2els(j)
539 ip = iparts(ie)
540 ipg = tagprt_fric(ip)
541 IF(ipg > 0 .AND. ip > ipfmax) THEN
543 . ipg,intbuf_fric_tab(intfric)%S_TABPARTS_FRIC,
544 . intbuf_fric_tab(intfric)%TABPARTS_FRIC,ipl )
545 IF(ipl /= 0) THEN
546 ipfmax = ip
547 ipflmax = ipl
548 ENDIF
549 ENDIF
550 ENDDO
551 IF(ipfmax /= 0) THEN
552 ipartfrics(i) = ipflmax
553 ENDIF
554
555 ENDDO
556 ENDIF
557
558 IF(numelc /= 0 .OR. numeltg /= 0) THEN
559 DO i = 1,nsn
560 ipfmax = 0
561 ipflmax = 0
562 DO j= knod2elc(nsv(i))+1,knod2elc(nsv(i)+1)
563 ie = nod2elc(j)
564 ip = ipartc(ie)
565 ipg = tagprt_fric(ip)
566 IF(ipg > 0 .AND. ip > ipfmax) THEN
568 . ipg,intbuf_fric_tab(intfric)%S_TABPARTS_FRIC,
569 . intbuf_fric_tab(intfric)%TABPARTS_FRIC,ipl )
570 IF(ipl /= 0) THEN
571 ipfmax = ip
572 ipflmax = ipl
573 ENDIF
574 ENDIF
575 ENDDO
576
577 DO j= knod2eltg(nsv(i))+1,knod2eltg(nsv(i)+1)
578 ie = nod2eltg(j)
579 ip = iparttg(ie)
580 ipg = tagprt_fric(ip)
581 IF(ipg > 0.AND.ip > ipfmax) THEN
583 . ipg,intbuf_fric_tab(intfric)%S_TABPARTS_FRIC,
584 . intbuf_fric_tab(intfric)%TABPARTS_FRIC,ipl )
585
586 IF(ipl /= 0) THEN
587 ipfmax = ip
588 ipflmax = ipl
589 ENDIF
590 ENDIF
591 ENDDO
592 IF(ipfmax /= 0) THEN
593 ipartfrics(i) = ipflmax
594 ENDIF
595
596 ENDDO
597 ENDIF
598 ENDIF
599
600
601 RETURN
void stlsort_int_int(int *len, int *keys, int *values)
subroutine area(d1, x, x2, y, y2, eint, stif0)
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 ineltc(nelc, neltg, is, surf_eltyp, surf_elem)
subroutine insol3(x, irect, ixs, nint, nel, i, area, noint, knod2els, nod2els, ir, ixs10, ixs16, ixs20)
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)