OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
i25cor3e.F File Reference
#include "implicit_f.inc"
#include "mvsiz_p.inc"
#include "param_c.inc"
#include "sms_c.inc"
#include "assert.inc"
#include "i25edge_c.inc"

Go to the source code of this file.

Functions/Subroutines

subroutine i25cor3e (jlt, ledge, irect, x, v, cand_s, cand_m, stfe, ms, ex, ey, ez, fx, fy, fz, stif, xxs1, xxs2, xys1, xys2, xzs1, xzs2, xxm1, xxm2, xym1, xym2, xzm1, xzm2, vxs1, vxs2, vys1, vys2, vzs1, vzs2, vxm1, vxm2, vym1, vym2, vzm1, vzm2, ms1, ms2, mm1, mm2, n1, n2, m1, m2, nedge, nin, stfac, nodnx_sms, nsms, gape, gapve, iedge, admsr, lbound, edg_bisector, vtx_bisector, igap0, iam, jam, ibm, jbm, ias, jas, ibs, jbs, itab, edge_id, intfric, ipartfric_e, ipartfricsi, ipartfricmi, igap, gap_e_l, igsti, kmin, kmax, istif_msdt, dtstif, stifmsdt_edg, parameters)

Function/Subroutine Documentation

◆ i25cor3e()

subroutine i25cor3e ( integer jlt,
integer, dimension(nledge,*) ledge,
integer, dimension(4,*) irect,
x,
v,
integer, dimension(*) cand_s,
integer, dimension(*) cand_m,
stfe,
ms,
ex,
ey,
ez,
fx,
fy,
fz,
stif,
xxs1,
xxs2,
xys1,
xys2,
xzs1,
xzs2,
xxm1,
xxm2,
xym1,
xym2,
xzm1,
xzm2,
vxs1,
vxs2,
vys1,
vys2,
vzs1,
vzs2,
vxm1,
vxm2,
vym1,
vym2,
vzm1,
vzm2,
ms1,
ms2,
mm1,
mm2,
integer, dimension(mvsiz) n1,
integer, dimension(mvsiz) n2,
integer, dimension(mvsiz) m1,
integer, dimension(mvsiz) m2,
integer nedge,
integer nin,
stfac,
integer, dimension(*) nodnx_sms,
integer, dimension(mvsiz) nsms,
gape,
gapve,
integer iedge,
integer, dimension(4,*) admsr,
integer, dimension(*) lbound,
real*4, dimension(3,4,*) edg_bisector,
real*4, dimension(3,2,*) vtx_bisector,
integer igap0,
integer, dimension(mvsiz) iam,
integer, dimension(mvsiz) jam,
integer, dimension(mvsiz) ibm,
integer, dimension(mvsiz) jbm,
integer, dimension(mvsiz) ias,
integer, dimension(mvsiz) jas,
integer, dimension(mvsiz) ibs,
integer, dimension(mvsiz) jbs,
integer, dimension(*) itab,
integer, dimension(2,mvsiz) edge_id,
integer intfric,
integer, dimension(*) ipartfric_e,
integer, dimension(mvsiz) ipartfricsi,
integer, dimension(mvsiz) ipartfricmi,
integer igap,
gap_e_l,
integer, intent(in) igsti,
intent(in) kmin,
intent(in) kmax,
integer, intent(in) istif_msdt,
intent(in) dtstif,
dimension(nedge), intent(in) stifmsdt_edg,
type (parameters_), intent(in) parameters )

Definition at line 32 of file i25cor3e.F.

51C-----------------------------------------------
52C M o d u l e s
53C-----------------------------------------------
54 USE tri25ebox
55 USE tri7box
56 USE parameters_mod
57C-----------------------------------------------
58C I m p l i c i t T y p e s
59C-----------------------------------------------
60#include "implicit_f.inc"
61C-----------------------------------------------
62C G l o b a l P a r a m e t e r s
63C-----------------------------------------------
64#include "mvsiz_p.inc"
65C-----------------------------------------------
66C C o m m o n B l o c k s
67C-----------------------------------------------
68#include "param_c.inc"
69#include "sms_c.inc"
70#include "assert.inc"
71#include "i25edge_c.inc"
72C-----------------------------------------------
73C D u m m y A r g u m e n t s
74C-----------------------------------------------
75 INTEGER LEDGE(NLEDGE,*), IRECT(4,*), CAND_M(*), CAND_S(*), ADMSR(4,*),
76 . LBOUND(*), JLT, NRTS, NIN, IEDGE, IGAP0,INTFRIC,IGAP,
77 . N1(MVSIZ), N2(MVSIZ),
78 . M1(MVSIZ), M2(MVSIZ),
79 . NODNX_SMS(*), NSMS(MVSIZ),
80 . ITAB(*),
81 . IAM(MVSIZ),JAM(MVSIZ),IBM(MVSIZ),JBM(MVSIZ),
82 . IAS(MVSIZ),JAS(MVSIZ),IBS(MVSIZ),JBS(MVSIZ),
83 . IPARTFRIC_E(*),IPARTFRICSI(MVSIZ),IPARTFRICMI(MVSIZ)
84 INTEGER :: EDGE_ID(2,MVSIZ)
85 INTEGER NEDGE
86 INTEGER , INTENT(IN) :: IGSTI
87 INTEGER , INTENT(IN) :: ISTIF_MSDT
88C REAL
89 my_real
90 . x(3,*), stfe(*), ms(*), v(3,*), gape(*),
91 . xxs1(*), xxs2(*), xys1(*), xys2(*),
92 . xzs1(*), xzs2(*), xxm1(*), xxm2(*),
93 . xym1(*), xym2(*), xzm1(*), xzm2(*),
94 . vxs1(*), vxs2(*), vys1(*), vys2(*),
95 . vzs1(*), vzs2(*), vxm1(*), vxm2(*),
96 . vym1(*), vym2(*), vzm1(*), vzm2(*),
97 . ms1(*), ms2(*), mm1(*), mm2(*),
98 . stif(*), stfac, sts, stm, gapve(*),
99 . ex(*), ey(*), ez(*), fx(*), fy(*), fz(*),
100 . gap_e_l(nedge)
101 real*4 edg_bisector(3,4,*), vtx_bisector(3,2,*)
102 my_real , INTENT(IN) :: kmin, kmax
103 my_real , INTENT(IN) :: dtstif
104 my_real , INTENT(IN) :: stifmsdt_edg(nedge)
105 TYPE (PARAMETERS_) ,INTENT(IN):: PARAMETERS
106C-----------------------------------------------
107C L o c a l V a r i a b l e s
108C-----------------------------------------------
109 INTEGER I ,NN, J, JRM, K, KRM, I1, J1, I2, J2,
110 . EM, IE, ES, JE, SOL_EDGE, SH_EDGE,
111 . TYPEDGS(MVSIZ), IM(MVSIZ), IS(MVSIZ)
112 INTEGER :: NOD1S(MVSIZ),NOD2S(MVSIZ)
113 INTEGER :: NOD1M(MVSIZ),NOD2M(MVSIZ)
114
115 INTEGER :: ISR ! orientation SECONDARY edge
116 my_real
117 . aaa, dx, dy, dz, dd, nni, ni2, invcos,dts
118 my_real
119 . gape_m(mvsiz), gape_s(mvsiz), stif_msdt(mvsiz)
120C-----------------------------------------------
121 DO i=1,jlt
122
123 em =cand_m(i)
124C debug ! global id of MAIN edge
125 edge_id(1,i) = ledge(8,em)
126
127 iam(i) = abs(ledge(1,em))
128 jam(i)=ledge(2,em)
129 ibm(i)=ledge(3,em)
130 jbm(i)=ledge(4,em)
131
132 m1(i)=ledge(5,em)
133 m2(i)=ledge(6,em)
134 im(i) = ledge(10,em)
135 nod1m(i) = ledge(11,em)
136 nod2m(i) = ledge(12,em)
137
138CC IF(IAM(I) > 0) THEN
139CC IF(IRECT(JAM(I),IAM(I))==M1(I).AND.IRECT(MOD(JAM(I),4)+1,IAM(I))==M2(I))THEN
140CC IM(I)= 1
141CC ELSEIF(IRECT(JAM(I),IAM(I))==M2(I).AND.IRECT(MOD(JAM(I),4)+1,IAM(I))==M1(I))THEN
142CC IM(I)=-1
143CC ELSE
144CC print *,'i25cor3e - internal problem',EM,M1(I),M2(I),
145CC . IRECT(JAM(I),IAM(I)),IRECT(MOD(JAM(I),4)+1,IAM(I))
146CC END IF
147CC ELSE
148CCC Faire le cas frontiere
149CC IM(I) = 1
150CC ENDIF
151CC ASSERT(IM(I) == LEDGE(10,EM))
152c IF ( IM(I) /= LEDGE(10,EM) ) THEN
153c WRITE(6,'(4I10)') I,LEDGE(10,EM),LEDGE(11,EM),IM(I)
154c ENDIF
155 stm=stfe(em)
156
157 xxm1(i) = x(1,m1(i))
158 xym1(i) = x(2,m1(i))
159 xzm1(i) = x(3,m1(i))
160 xxm2(i) = x(1,m2(i))
161 xym2(i) = x(2,m2(i))
162 xzm2(i) = x(3,m2(i))
163 vxm1(i) = v(1,m1(i))
164 vym1(i) = v(2,m1(i))
165 vzm1(i) = v(3,m1(i))
166 vxm2(i) = v(1,m2(i))
167 vym2(i) = v(2,m2(i))
168 vzm2(i) = v(3,m2(i))
169 mm1(i) = ms(m1(i))
170 mm2(i) = ms(m2(i))
171C
172 IF(cand_s(i)<=nedge) THEN
173
174 es =cand_s(i)
175 edge_id(2,i) = ledge(8,es)
176 ias(i)=abs( ledge(1,es) )
177 jas(i)=ledge(2,es)
178 ibs(i)=ledge(3,es)
179 jbs(i)=ledge(4,es)
180 n1(i)=ledge(5,es)
181 n2(i)=ledge(6,es)
182 nod1s(i) = ledge(11,es)
183 nod2s(i) = ledge(12,es)
184
185 is(i) = ledge(10,es)
186CC IF(IAS(I) > 0) THEN
187CC IF(IRECT(JAS(I),IAS(I))==N1(I).AND.IRECT(MOD(JAS(I),4)+1,IAS(I))==N2(I))THEN
188CC IS(I)= 1
189CC ELSEIF(IRECT(JAS(I),IAS(I))==N2(I).AND.IRECT(MOD(JAS(I),4)+1,IAS(I))==N1(I))THEN
190CC IS(I)=-1
191CC ELSE
192CC print *,'i25cor3e - internal problem',ES,N1(I),N2(I),
193CC . IRECT(JAS(I),IAS(I)),IRECT(MOD(JAS(I),4)+1,IAS(I))
194CC END IF
195CC ELSE
196CCC faire le SPMD
197CC IS(I) = 1
198CC ENDIF
199CC ASSERT(IS(I) == LEDGE(10,ES))
200
201 sts=stfe(es)
202 stif(i)=sts*stm / max(em20,sts+stm)
203c STIF(I)=MAX(KMIN,MIN(STIF(I),KMAX))
204 xxs1(i) = x(1,n1(i))
205 xys1(i) = x(2,n1(i))
206 xzs1(i) = x(3,n1(i))
207 xxs2(i) = x(1,n2(i))
208 xys2(i) = x(2,n2(i))
209 xzs2(i) = x(3,n2(i))
210 vxs1(i) = v(1,n1(i))
211 vys1(i) = v(2,n1(i))
212 vzs1(i) = v(3,n1(i))
213 vxs2(i) = v(1,n2(i))
214 vys2(i) = v(2,n2(i))
215 vzs2(i) = v(3,n2(i))
216 ms1(i) = ms(n1(i))
217 ms2(i) = ms(n2(i))
218
219C
220 typedgs(i)=ledge(7,es)
221C
222 ELSE
223 ! IBUF_EDGE(P)%p(E_LEFT_SEG + L) = LEDGE(1,I)
224 ! IBUF_EDGE(P)%p(E_LEFT_ID + L) = LEDGE(2,I)
225 ! IBUF_EDGE(P)%p(E_RIGHT_SEG + L) = LEDGE(3,I)
226 ! IBUF_EDGE(P)%p(E_RIGHT_ID + L) = LEDGE(4,I)
227 ! IBUF_EDGE(P)%p(E_NODE1_ID + L) = LEDGE(5,I)
228 ! IBUF_EDGE(P)%p(E_NODE2_ID + L) = LEDGE(6,I)
229 ! IBUF_EDGE(P)%p(E_TYPE + L) = LEDGE(7,I)
230
231 nn = cand_s(i) - nedge
232 is(i) = ledge_fie(nin)%P(e_im,nn)
233C IF(IS(i) == 1) THEN
234 n1(i)=2*(nn-1)+1
235 n2(i)=2*nn
236c ELSE
237c N2(I)=2*(NN-1)+1
238c N1(I)=2*NN
239c ENDIF
240
241 edge_id(2,i) = ledge_fie(nin)%P(e_global_id,nn)
242
243 ias(i)=abs( ledge_fie(nin)%P(e_left_seg ,nn) )
244 jas(i)=ledge_fie(nin)%P(e_left_id ,nn)
245 ibs(i)=ledge_fie(nin)%P(e_right_seg ,nn)
246 jbs(i)=ledge_fie(nin)%P(e_right_id ,nn)
247 typedgs(i)=ledge_fie(nin)%P(e_type,nn)
248 sts=stifie(nin)%P(nn)
249 stif(i)=sts*stm / max(em20,sts+stm)
250c STIF(I)=MAX(KMIN,MIN(STIF(I),KMAX))
251
252
253c STS=STFE(ES)
254c STIF(I)=STS*STM / MAX(EM20,STS+STM)
255c STIF(I)=ABS(STIFIE(NIN)%P(NN))*STFM
256c / MAX(EM20,ABS(STIFIE(NIN)%P(NN))+STM)
257c
258c TYPEDGS(I)=LEDGE(7,CAND_S(I))
259c
260 xxs1(i) = xfie(nin)%P(1,n1(i))
261 xys1(i) = xfie(nin)%P(2,n1(i))
262 xzs1(i) = xfie(nin)%P(3,n1(i))
263 xxs2(i) = xfie(nin)%P(1,n2(i))
264 xys2(i) = xfie(nin)%P(2,n2(i))
265 xzs2(i) = xfie(nin)%P(3,n2(i))
266 vxs1(i) = vfie(nin)%P(1,n1(i))
267 vys1(i) = vfie(nin)%P(2,n1(i))
268 vzs1(i) = vfie(nin)%P(3,n1(i))
269 vxs2(i) = vfie(nin)%P(1,n2(i))
270 vys2(i) = vfie(nin)%P(2,n2(i))
271 vzs2(i) = vfie(nin)%P(3,n2(i))
272 ms1(i) = msfie(nin)%P(n1(i))
273 ms2(i) = msfie(nin)%P(n2(i))
274 END IF
275c DEBUG_E2E(EDGE_ID(1,I) == D_EM.AND. EDGE_ID(2,I) == D_ES,XXS1(i))
276c DEBUG_E2E(EDGE_ID(1,I) == D_EM.AND. EDGE_ID(2,I) == D_ES,XXS2(i))
277 END DO
278
279C------------------------------------------
280C Stiffness based on mass and time step
281C------------------------------------------
282
283 IF(istif_msdt > 0) THEN
284 IF(dtstif > zero) THEN
285 dts = dtstif
286 ELSE
287 dts = parameters%DT_STIFINT
288 ENDIF
289 DO i=1,jlt
290 em =cand_m(i)
291 IF(cand_s(i)<=nedge) THEN
292 es =cand_s(i)
293 stif_msdt(i) = stifmsdt_edg(es)
294 ELSE
295 nn = cand_s(i) - nedge
296 stif_msdt(i) = abs(stife_msdt_fi(nin)%P(nn))
297 ENDIF
298 stif_msdt(i) = stifmsdt_edg(em)*stif_msdt(i)/(stifmsdt_edg(em)+stif_msdt(i))
299
300 stif_msdt(i) = stif_msdt(i)/(dts*dts)
301 stif(i)=max(stif(i),stif_msdt(i))
302 ENDDO
303 ENDIF
304C
305 DO i=1,jlt
306 stif(i)=max(kmin,min(stif(i),kmax))
307 ENDDO
308
309 sol_edge=iedge/10 ! solids
310 sh_edge =iedge-10*sol_edge ! shells
311
312 ex(1:jlt)=zero
313 ey(1:jlt)=zero
314 ez(1:jlt)=zero
315 fx(1:jlt)=zero
316 fy(1:jlt)=zero
317 fz(1:jlt)=zero
318 IF(sh_edge/=0)THEN
319 DO i=1,jlt
320
321CC BOUNDARY_EDGE :
322C IF(IAM(I) < 0) CYCLE
323C
324 ex(i) = edg_bisector(1,jam(i),iam(i))
325 ey(i) = edg_bisector(2,jam(i),iam(i))
326 ez(i) = edg_bisector(3,jam(i),iam(i))
327C
328 IF(iabs(typedgs(i))/=1)THEN
329 IF(cand_s(i)<=nedge) THEN
330 fx(i) = edg_bisector(1,jas(i),ias(i))
331 fy(i) = edg_bisector(2,jas(i),ias(i))
332 fz(i) = edg_bisector(3,jas(i),ias(i))
333 ELSE
334 fx(i) = edg_bisector_fie(nin)%P(1,1,cand_s(i)-nedge)
335 fy(i) = edg_bisector_fie(nin)%P(2,1,cand_s(i)-nedge)
336 fz(i) = edg_bisector_fie(nin)%P(3,1,cand_s(i)-nedge)
337 END IF
338 END IF
339 END DO
340 END IF
341
342 DO i=1,jlt
343 gape_m(i)=gape(cand_m(i))
344 IF(cand_s(i)<=nedge) THEN
345 gape_s(i)=gape(cand_s(i))
346 ELSE
347 gape_s(i)= gapfie(nin)%P(cand_s(i) - nedge)
348 END IF
349 debug_e2e(edge_id(1,i) == d_em .AND. edge_id(2,i) == d_es,gape_s(i))
350 gapve(i)=gape_m(i)+gape_s(i)
351 END DO
352 IF(igap == 3) THEN
353 DO i=1,jlt
354 gape_m(i)=min(gape_m(i),gap_e_l(cand_m(i)))
355 IF(cand_s(i)<=nedge) THEN
356 gape_s(i)=min(gape_s(i),gap_e_l(cand_s(i)))
357 gapve(i)=min(gap_e_l(cand_m(i))+gap_e_l(cand_s(i)),gapve(i))
358 ELSE
359 gape_s(i)=min(gape_s(i),gape_l_fie(nin)%P(cand_s(i) - nedge))
360 gapve(i)=min(gap_e_l(cand_m(i))+gape_l_fie(nin)%P(cand_s(i) - nedge),gapve(i))
361 ENDIF
362 ENDDO
363 ENDIF
364C
365C Always shift MAIN edge toward inner side
366 DO i=1,jlt
367
368 IF(ibm(i)/=0)cycle
369
370 ie=iam(i)
371 je=jam(i)
372
373 i1 = nod1m(i)
374 i2 = nod2m(i)
375 ex(i) = edg_bisector(1,je,ie)
376 ey(i) = edg_bisector(2,je,ie)
377 ez(i) = edg_bisector(3,je,ie)
378C
379C DX, DY, DZ Direction for shifting
380 dx = vtx_bisector(1,1,i1)+vtx_bisector(1,2,i1)
381 dy = vtx_bisector(2,1,i1)+vtx_bisector(2,2,i1)
382 dz = vtx_bisector(3,1,i1)+vtx_bisector(3,2,i1)
383C
384 nni = ex(i)*dx + ey(i)*dy + ez(i)*dz
385 ni2 = dx*dx + dy*dy + dz*dz
386
387 debug_e2e( edge_id(1,i) == d_em.AND. edge_id(2,i) == d_es ,ex(i))
388 debug_e2e( edge_id(1,i) == d_em.AND. edge_id(2,i) == d_es ,ey(i))
389 debug_e2e( edge_id(1,i) == d_em.AND. edge_id(2,i) == d_es ,ez(i))
390 debug_e2e( edge_id(1,i) == d_em.AND. edge_id(2,i) == d_es ,dx)
391 debug_e2e( edge_id(1,i) == d_em.AND. edge_id(2,i) == d_es ,dy)
392 debug_e2e( edge_id(1,i) == d_em.AND. edge_id(2,i) == d_es ,dz)
393
394 IF(nni < zero)THEN
395 dx=dx-two*nni*ex(i)
396 dy=dy-two*nni*ey(i)
397 dz=dz-two*nni*ez(i)
398 nni=-nni
399 END IF
400
401 IF(two*nni*nni < ni2)THEN
402c scharp angle bound nodal normal to 45 from edge normal
403 aaa = sqrt(max(zero,ni2-nni*nni)) - nni
404 dx = dx + aaa*ex(i)
405 dy = dy + aaa*ey(i)
406 dz = dz + aaa*ez(i)
407 ENDIF
408 dd=one/max(em20,sqrt(dx*dx+dy*dy+dz*dz))
409 dx = dx*dd
410 dy = dy*dd
411 dz = dz*dd
412 invcos = one / max(em20,ex(i)*dx + ey(i)*dy + ez(i)*dz)
413 dx = dx*invcos
414 dy = dy*invcos
415 dz = dz*invcos
416C
417 xxm1(i) = xxm1(i)-gape_m(i)*dx
418 xym1(i) = xym1(i)-gape_m(i)*dy
419 xzm1(i) = xzm1(i)-gape_m(i)*dz
420C
421C DX, DY, DZ Direction for shifting
422 dx = vtx_bisector(1,1,i2)+vtx_bisector(1,2,i2)
423 dy = vtx_bisector(2,1,i2)+vtx_bisector(2,2,i2)
424 dz = vtx_bisector(3,1,i2)+vtx_bisector(3,2,i2)
425C
426 nni = ex(i)*dx + ey(i)*dy + ez(i)*dz
427 ni2 = dx*dx + dy*dy + dz*dz
428
429 IF(nni < zero)THEN
430 dx=dx-two*nni*ex(i)
431 dy=dy-two*nni*ey(i)
432 dz=dz-two*nni*ez(i)
433 nni=-nni
434 END IF
435
436 IF(two*nni*nni < ni2)THEN
437c scharp angle bound nodal normal to 45 from edge normal
438 aaa = sqrt(max(zero,ni2-nni*nni)) - nni
439 dx = dx + aaa*ex(i)
440 dy = dy + aaa*ey(i)
441 dz = dz + aaa*ez(i)
442 ENDIF
443 dd=one/max(em20,sqrt(dx*dx+dy*dy+dz*dz))
444 dx = dx*dd
445 dy = dy*dd
446 dz = dz*dd
447 invcos = one / max(em20,ex(i)*dx + ey(i)*dy + ez(i)*dz)
448 dx = dx*invcos
449 dy = dy*invcos
450 dz = dz*invcos
451C
452 xxm2(i) = xxm2(i)-gape_m(i)*dx
453 xym2(i) = xym2(i)-gape_m(i)*dy
454 xzm2(i) = xzm2(i)-gape_m(i)*dz
455C
456 END DO
457
458 IF(igap0/=0)THEN
459C
460C Shift SECONDARY edge toward inner side
461 DO i=1,jlt
462 debug_e2e(edge_id(1,i) == d_em .AND. edge_id(2,i) == d_es,ibs(i))
463
464 IF(ibs(i)/=0)cycle
465
466 IF(cand_s(i)<=nedge) THEN
467 ie=ias(i)
468 je=jas(i)
469
470 i1 = nod1s(i)
471 i2 = nod2s(i)
472
473 fx(i) = edg_bisector(1,je,ie)
474 fy(i) = edg_bisector(2,je,ie)
475 fz(i) = edg_bisector(3,je,ie)
476C
477C DX, DY, DZ Direction for shifting
478 dx = vtx_bisector(1,1,i1)+vtx_bisector(1,2,i1)
479 dy = vtx_bisector(2,1,i1)+vtx_bisector(2,2,i1)
480 dz = vtx_bisector(3,1,i1)+vtx_bisector(3,2,i1)
481 assert(fx(i) == fx(i))
482 assert(fy(i) == fy(i))
483 assert(fz(i) == fz(i))
484 ELSE ! edge SECONDARY remote
485 fx(i) = edg_bisector_fie(nin)%P(1,1,cand_s(i)-nedge)
486 fy(i) = edg_bisector_fie(nin)%P(2,1,cand_s(i)-nedge)
487 fz(i) = edg_bisector_fie(nin)%P(3,1,cand_s(i)-nedge)
488
489 assert(fx(i) == fx(i))
490 assert(fy(i) == fy(i))
491 assert(fz(i) == fz(i))
492C DX, DY, DZ Direction for shifting
493 dx = vtx_bisector_fie(nin)%P(1,1,cand_s(i)-nedge)+vtx_bisector_fie(nin)%P(1,2,cand_s(i)-nedge)
494 dy = vtx_bisector_fie(nin)%P(2,1,cand_s(i)-nedge)+vtx_bisector_fie(nin)%P(2,2,cand_s(i)-nedge)
495 dz = vtx_bisector_fie(nin)%P(3,1,cand_s(i)-nedge)+vtx_bisector_fie(nin)%P(3,2,cand_s(i)-nedge)
496 assert(dx == dx)
497 assert(dy == dy)
498 assert(dz == dz)
499 ENDIF
500
501 nni = fx(i)*dx + fy(i)*dy + fz(i)*dz
502 ni2 = dx*dx + dy*dy + dz*dz
503
504 IF(nni < zero)THEN
505 dx=dx-two*nni*fx(i)
506 dy=dy-two*nni*fy(i)
507 dz=dz-two*nni*fz(i)
508 nni=-nni
509 END IF
510
511 IF(two*nni*nni < ni2)THEN
512c scharp angle bound nodal normal to 45 from edge normal
513 aaa = sqrt(max(zero,ni2-nni*nni)) - nni
514 dx = dx + aaa*fx(i)
515 dy = dy + aaa*fy(i)
516 dz = dz + aaa*fz(i)
517 ENDIF
518 dd=one/max(em20,sqrt(dx*dx+dy*dy+dz*dz))
519 dx = dx*dd
520 dy = dy*dd
521 dz = dz*dd
522 invcos = one / max(em20,fx(i)*dx + fy(i)*dy + fz(i)*dz)
523 dx = dx*invcos
524 dy = dy*invcos
525 dz = dz*invcos
526C
527 xxs1(i) = xxs1(i)-gape_s(i)*dx
528 xys1(i) = xys1(i)-gape_s(i)*dy
529 xzs1(i) = xzs1(i)-gape_s(i)*dz
530C
531 IF(cand_s(i)<=nedge) THEN
532C DX, DY, DZ Direction for shifting
533 dx = vtx_bisector(1,1,i2)+vtx_bisector(1,2,i2)
534 dy = vtx_bisector(2,1,i2)+vtx_bisector(2,2,i2)
535 dz = vtx_bisector(3,1,i2)+vtx_bisector(3,2,i2)
536 assert(dx == dx)
537 assert(dy == dy)
538 assert(dz == dz)
539 ELSE
540 dx = vtx_bisector_fie(nin)%P(1,3,cand_s(i)-nedge)+vtx_bisector_fie(nin)%P(1,4,cand_s(i)-nedge)
541 dy = vtx_bisector_fie(nin)%P(2,3,cand_s(i)-nedge)+vtx_bisector_fie(nin)%P(2,4,cand_s(i)-nedge)
542 dz = vtx_bisector_fie(nin)%P(3,3,cand_s(i)-nedge)+vtx_bisector_fie(nin)%P(3,4,cand_s(i)-nedge)
543 assert(dx == dx)
544 assert(dy == dy)
545 assert(dz == dz)
546 ENDIF
547C
548 nni = fx(i)*dx + fy(i)*dy + fz(i)*dz
549 ni2 = dx*dx + dy*dy + dz*dz
550
551 IF(nni < zero)THEN
552 dx=dx-two*nni*fx(i)
553 dy=dy-two*nni*fy(i)
554 dz=dz-two*nni*fz(i)
555 nni=-nni
556 END IF
557
558 IF(two*nni*nni < ni2)THEN
559c scharp angle bound nodal normal to 45 from edge normal
560 aaa = sqrt(max(zero,ni2-nni*nni)) - nni
561 dx = dx + aaa*fx(i)
562 dy = dy + aaa*fy(i)
563 dz = dz + aaa*fz(i)
564 ENDIF
565 dd=one/max(em20,sqrt(dx*dx+dy*dy+dz*dz))
566 dx = dx*dd
567 dy = dy*dd
568 dz = dz*dd
569 invcos = one / max(em20,fx(i)*dx + fy(i)*dy + fz(i)*dz)
570 dx = dx*invcos
571 dy = dy*invcos
572 dz = dz*invcos
573C
574 xxs2(i) = xxs2(i)-gape_s(i)*dx
575 xys2(i) = xys2(i)-gape_s(i)*dy
576 xzs2(i) = xzs2(i)-gape_s(i)*dz
577C
578 END DO
579 END IF
580C
581 IF(idtmins==2)THEN
582 DO i=1,jlt
583 IF(cand_s(i)<=nedge)THEN
584 nsms(i)=nodnx_sms(n1(i))+nodnx_sms(n2(i))+
585 . nodnx_sms(m1(i))+nodnx_sms(m2(i))
586 ELSE
587 nsms(i)=nodnxfie(nin)%P(n1(i))+nodnxfie(nin)%P(n2(i))+
588 . nodnx_sms(m1(i))+nodnx_sms(m2(i))
589 END IF
590 ENDDO
591 IF(idtmins_int/=0)THEN
592 DO i=1,jlt
593 IF(nsms(i)==0)nsms(i)=-1
594 ENDDO
595 END IF
596 ELSEIF(idtmins_int/=0)THEN
597 DO i=1,jlt
598 nsms(i)=-1
599 ENDDO
600 ENDIF
601C
602C----Friction model : secnd part IDs---------
603 IF(intfric > 0) THEN
604 DO i=1,jlt
605
606 IF(cand_s(i)<=nedge)THEN
607 ipartfricsi(i) = ipartfric_e(cand_s(i))
608 ELSE
609 nn = cand_s(i) - nedge
610 ipartfricsi(i)= ipartfric_fie(nin)%P(nn)
611 ENDIF
612
613C
614 ipartfricmi(i) = ipartfric_e(cand_m(i))
615 ENDDO
616 ENDIF
617 RETURN
my_real
#define my_real
Definition cppsort.cpp:32
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
parameters_mod
Definition parameters_mod.F:49
tri25ebox::gape_l_fie
type(real_pointer), dimension(:), allocatable gape_l_fie
Definition tri25ebox.F:86
tri25ebox::edg_bisector_fie
type(real4_pointer3), dimension(:), allocatable edg_bisector_fie
Definition tri25ebox.F:83
tri25ebox::vtx_bisector_fie
type(real4_pointer3), dimension(:), allocatable vtx_bisector_fie
Definition tri25ebox.F:84
tri25ebox::ledge_fie
type(int_pointer2), dimension(:), allocatable ledge_fie
Definition tri25ebox.F:88
tri7box::gapfie
type(real_pointer), dimension(:), allocatable gapfie
Definition tri7box.F:449
tri7box::vfie
type(real_pointer2), dimension(:), allocatable vfie
Definition tri7box.F:459
tri7box::ipartfric_fie
type(int_pointer), dimension(:), allocatable ipartfric_fie
Definition tri7box.F:440
tri7box::xfie
type(real_pointer2), dimension(:), allocatable xfie
Definition tri7box.F:459
tri7box::stifie
type(real_pointer), dimension(:), allocatable stifie
Definition tri7box.F:449
tri7box::nodnxfie
type(int_pointer), dimension(:), allocatable nodnxfie
Definition tri7box.F:440
tri7box::stife_msdt_fi
type(real_pointer), dimension(:), allocatable stife_msdt_fi
Definition tri7box.F:553
tri7box::msfie
type(real_pointer), dimension(:), allocatable msfie
Definition tri7box.F:449