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