OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
i25cor3e.F
Go to the documentation of this file.
1Copyright> OpenRadioss
2Copyright> Copyright (C) 1986-2025 Altair Engineering Inc.
3Copyright>
4Copyright> This program is free software: you can redistribute it and/or modify
5Copyright> it under the terms of the GNU Affero General Public License as published by
6Copyright> the Free Software Foundation, either version 3 of the License, or
7Copyright> (at your option) any later version.
8Copyright>
9Copyright> This program is distributed in the hope that it will be useful,
10Copyright> but WITHOUT ANY WARRANTY; without even the implied warranty of
11Copyright> MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12Copyright> GNU Affero General Public License for more details.
13Copyright>
14Copyright> You should have received a copy of the GNU Affero General Public License
15Copyright> along with this program. If not, see <https://www.gnu.org/licenses/>.
16Copyright>
17Copyright>
18Copyright> Commercial Alternative: Altair Radioss Software
19Copyright>
20Copyright> As an alternative to this open-source version, Altair also offers Altair Radioss
21Copyright> software under a commercial license. Contact Altair to discuss further if the
22Copyright> commercial version may interest you: https://www.altair.com/radioss/.
23!||====================================================================
24!|| i25cor3e ../engine/source/interfaces/int25/i25cor3e.F
25!||--- called by ------------------------------------------------------
26!|| i25mainf ../engine/source/interfaces/int25/i25mainf.F
27!||--- uses -----------------------------------------------------
28!|| parameters_mod ../common_source/modules/interfaces/parameters_mod.f
29!|| tri25ebox ../engine/share/modules/tri25ebox.F
30!|| tri7box ../engine/share/modules/tri7box.F
31!||====================================================================
32 SUBROUTINE i25cor3e(
33 1 JLT ,LEDGE ,IRECT ,X ,V ,
34 2 CAND_S ,CAND_M ,STFE ,MS ,EX ,
35 3 EY ,EZ ,FX ,FY ,FZ ,
36 4 STIF ,XXS1 ,XXS2 ,XYS1 ,XYS2 ,
37 5 XZS1 ,XZS2 ,XXM1 ,XXM2 ,XYM1 ,
38 6 XYM2 ,XZM1 ,XZM2 ,VXS1 ,VXS2 ,
39 7 VYS1 ,VYS2 ,VZS1 ,VZS2 ,VXM1 ,
40 8 VXM2 ,VYM1 ,VYM2 ,VZM1 ,VZM2 ,
41 9 MS1 ,MS2 ,MM1 ,MM2 ,N1 ,
42 A N2 ,M1 ,M2 ,NEDGE ,NIN ,
43 C STFAC ,NODNX_SMS,NSMS ,GAPE ,GAPVE ,
44 D IEDGE ,ADMSR ,LBOUND ,EDG_BISECTOR ,
45 E VTX_BISECTOR,IGAP0,
46 F IAM ,JAM ,IBM ,JBM ,IAS ,
47 G JAS ,IBS ,JBS ,ITAB , EDGE_ID,
48 H INTFRIC,IPARTFRIC_E,IPARTFRICSI,IPARTFRICMI,
49 I IGAP ,GAP_E_L,IGSTI ,KMIN ,KMAX ,
50 J ISTIF_MSDT,DTSTIF,STIFMSDT_EDG ,PARAMETERS)
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
618 END
i25cor3e
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)
Definition i25cor3e.F:51
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