OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
i7keg3.F File Reference
#include "implicit_f.inc"
#include "mvsiz_p.inc"
#include "scr05_c.inc"
#include "com01_c.inc"
#include "impl1_c.inc"

Go to the source code of this file.

Functions/Subroutines

subroutine i7keg3 (jlt, a, v, ms, fric, nx1, nx2, nx3, nx4, ny1, ny2, ny3, ny4, nz1, nz2, nz3, nz4, lb1, lb2, lb3, lb4, lc1, lc2, lc3, lc4, p1, p2, p3, p4, nin, ix1, ix2, ix3, ix4, nsvg, gapv, inacti, cand_p, index, stiglo, stif, vxi, vyi, vzi, msi, ki11, ki12, kj11, kj12, kk11, kk12, kl11, kl12, off, scalk, lrem, icurv, x, x1, x2, x3, x4, y1, y2, y3, y4, z1, z2, z3, z4, xi, yi, zi)
subroutine i7frf3 (jlt, a, v, ms, fric, n1, n2, n3, h1, h2, h3, h4, ix1, ix2, ix3, ix4, index, vxi, vyi, vzi, msi, dxi, dyi, dzi, stif, nin, d, scalk)
subroutine i7kfor3 (jlt, a, v, ms, fric, nx1, nx2, nx3, nx4, ny1, ny2, ny3, ny4, nz1, nz2, nz3, nz4, lb1, lb2, lb3, lb4, lc1, lc2, lc3, lc4, p1, p2, p3, p4, nin, ix1, ix2, ix3, ix4, nsvg, gapv, inacti, cand_p, index, stiglo, stif, vxi, vyi, vzi, msi, x, irect, ce_loc, mfrot, ifq, frot_p, cand_fx, cand_fy, cand_fz, alpha0, dxi, dyi, dzi, d, scalk, x1, x2, x3, x4, y1, y2, y3, y4, z1, z2, z3, z4, xi, yi, zi, icurv)

Function/Subroutine Documentation

◆ i7frf3()

subroutine i7frf3 ( integer jlt,
a,
v,
ms,
fric,
n1,
n2,
n3,
h1,
h2,
h3,
h4,
integer, dimension(mvsiz) ix1,
integer, dimension(mvsiz) ix2,
integer, dimension(mvsiz) ix3,
integer, dimension(mvsiz) ix4,
integer, dimension(mvsiz) index,
vxi,
vyi,
vzi,
msi,
dxi,
dyi,
dzi,
stif,
integer nin,
d,
scalk )

Definition at line 551 of file i7keg3.F.

557C-----------------------------------------------
558C M o d u l e s
559C-----------------------------------------------
560 USE imp_intm
561C-----------------------------------------------
562C I m p l i c i t T y p e s
563C-----------------------------------------------
564#include "implicit_f.inc"
565C-----------------------------------------------
566C G l o b a l P a r a m e t e r s
567C-----------------------------------------------
568#include "mvsiz_p.inc"
569C-----------------------------------------------
570C D u m m y A r g u m e n t s
571C-----------------------------------------------
572 INTEGER JLT, INACTI,NIN
573 INTEGER IX1(MVSIZ), IX2(MVSIZ), IX3(MVSIZ), IX4(MVSIZ),
574 . INDEX(MVSIZ)
575 my_real
576 . a(3,*), ms(*), v(3,*),d(3,*),
577 . fric,scalk,dxi(mvsiz),dyi(mvsiz),dzi(mvsiz),
578 . h1(mvsiz), h2(mvsiz), h3(mvsiz), h4(mvsiz),
579 . vxi(mvsiz),vyi(mvsiz),vzi(mvsiz),msi(mvsiz)
580 my_real
581 . n1(mvsiz), n2(mvsiz), n3(mvsiz),stif(mvsiz)
582C-----------------------------------------------
583C L o c a l V a r i a b l e s
584C-----------------------------------------------
585 INTEGER I, J1, J, K,IG,ISF,NN,NS,NI
586 my_real
587 . vx(mvsiz), vy(mvsiz), vz(mvsiz), vn(mvsiz),
588 . dx(mvsiz), dy(mvsiz), dz(mvsiz), dn(mvsiz),
589 . dni(mvsiz),dt(mvsiz), dti(mvsiz),
590 . s2,facn(mvsiz),facf, fact(mvsiz)
591 my_real
592 . fx,fy,fz,fn,ft,fni,fti,vtx,vty,vtz,vt,
593 . t1(mvsiz), t2(mvsiz), t3(mvsiz),
594 . kt1,kt2,kt3,kt4,q1,q2,h0
595
596C-----------------------------------------------
597C
598 DO i=1,jlt
599 vx(i) = vxi(i) - h1(i)*v(1,ix1(i)) - h2(i)*v(1,ix2(i))
600 . - h3(i)*v(1,ix3(i)) - h4(i)*v(1,ix4(i))
601 vy(i) = vyi(i) - h1(i)*v(2,ix1(i)) - h2(i)*v(2,ix2(i))
602 . - h3(i)*v(2,ix3(i)) - h4(i)*v(2,ix4(i))
603 vz(i) = vzi(i) - h1(i)*v(3,ix1(i)) - h2(i)*v(3,ix2(i))
604 . - h3(i)*v(3,ix3(i)) - h4(i)*v(3,ix4(i))
605 vn(i) = n1(i)*vx(i) + n2(i)*vy(i) + n3(i)*vz(i)
606 dx(i) = dxi(i) - h1(i)*d(1,ix1(i)) - h2(i)*d(1,ix2(i))
607 . - h3(i)*d(1,ix3(i)) - h4(i)*d(1,ix4(i))
608 dy(i) = dyi(i) - h1(i)*d(2,ix1(i)) - h2(i)*d(2,ix2(i))
609 . - h3(i)*d(2,ix3(i)) - h4(i)*d(2,ix4(i))
610 dz(i) = dzi(i) - h1(i)*d(3,ix1(i)) - h2(i)*d(3,ix2(i))
611 . - h3(i)*d(3,ix3(i)) - h4(i)*d(3,ix4(i))
612 dn(i) = n1(i)*dx(i) + n2(i)*dy(i) + n3(i)*dz(i)
613 dni(i) = n1(i)*dxi(i) + n2(i)*dyi(i) + n3(i)*dzi(i)
614 ENDDO
615C-------------------------------------------
616 DO i=1,jlt
617 vtx = vx(i) -vn(i)*n1(i)
618 vty = vy(i) -vn(i)*n2(i)
619 vtz = vz(i) -vn(i)*n3(i)
620 vt = vtx*vtx+vty*vty+vtz*vtz
621 IF (vt>em20) THEN
622 s2=one/sqrt(vt)
623 t1(i)=vtx*s2
624 t2(i)=vty*s2
625 t3(i)=vtz*s2
626 fact(i)=fric
627 ELSE
628 fact(i)=zero
629 t1(i)=one
630 t2(i)=zero
631 t3(i)=zero
632 ENDIF
633 ENDDO
634 DO i=1,jlt
635 dt(i) = t1(i)*dx(i) + t2(i)*dy(i) + t3(i)*dz(i)
636 dti(i) = t1(i)*dxi(i) + t2(i)*dyi(i) + t3(i)*dzi(i)
637 ENDDO
638 IF (scalk<0) THEN
639 isf=1
640 ELSE
641 isf=0
642 ENDIF
643 facf=abs(scalk)
644 IF (isf==1) THEN
645 DO i=1,jlt
646 IF (vn(i)>zero) THEN
647 facn(i)=stif(i)*facf
648 ELSEIF (vn(i)<zero) THEN
649 facn(i)=stif(i)/facf
650 ELSE
651 facn(i)=stif(i)
652 ENDIF
653 fact(i)=facn(i)*fact(i)
654 ENDDO
655 ELSE
656 DO i=1,jlt
657 facn(i)=stif(i)*facf
658 fact(i)=facn(i)*fact(i)
659 ENDDO
660 ENDIF
661C--------partie NML-------
662 DO i=1,jlt
663 fn = -facn(i)*dni(i)
664 fx=fn*n1(i)
665 fy=fn*n2(i)
666 fz=fn*n3(i)
667 IF (fact(i)/=zero) THEN
668 ft = -fact(i)*dti(i)
669 fx = fx + ft*t1(i)
670 fy = fy + ft*t2(i)
671 fz = fz + ft*t3(i)
672 ENDIF
673 a(1,ix1(i))=a(1,ix1(i))+fx*h1(i)
674 a(1,ix2(i))=a(1,ix2(i))+fx*h2(i)
675 a(1,ix3(i))=a(1,ix3(i))+fx*h3(i)
676 a(1,ix4(i))=a(1,ix4(i))+fx*h4(i)
677 a(2,ix1(i))=a(2,ix1(i))+fy*h1(i)
678 a(2,ix2(i))=a(2,ix2(i))+fy*h2(i)
679 a(2,ix3(i))=a(2,ix3(i))+fy*h3(i)
680 a(2,ix4(i))=a(2,ix4(i))+fy*h4(i)
681 a(3,ix1(i))=a(3,ix1(i))+fz*h1(i)
682 a(3,ix2(i))=a(3,ix2(i))+fz*h2(i)
683 a(3,ix3(i))=a(3,ix3(i))+fz*h3(i)
684 a(3,ix4(i))=a(3,ix4(i))+fz*h4(i)
685 ENDDO
686C--------partie NSL-------
687 DO i=1,jlt
688 fni = facn(i)*dn(i)
689 ni = index(i)
690 ffi(1,ni)=ffi(1,ni)+fni*n1(i)
691 ffi(2,ni)=ffi(2,ni)+fni*n2(i)
692 ffi(3,ni)=ffi(3,ni)+fni*n3(i)
693 IF (fact(i)/=zero) THEN
694 fti = fact(i)*dt(i)
695 ffi(1,ni)=ffi(1,ni)+fti*t1(i)
696 ffi(2,ni)=ffi(2,ni)+fti*t2(i)
697 ffi(3,ni)=ffi(3,ni)+fti*t3(i)
698 ENDIF
699 ENDDO
700C
701 RETURN
my_real
#define my_real
Definition cppsort.cpp:32
imp_intm
Definition imp_intm.F:118

◆ i7keg3()

subroutine i7keg3 ( integer jlt,
a,
v,
ms,
fric,
nx1,
nx2,
nx3,
nx4,
ny1,
ny2,
ny3,
ny4,
nz1,
nz2,
nz3,
nz4,
lb1,
lb2,
lb3,
lb4,
lc1,
lc2,
lc3,
lc4,
p1,
p2,
p3,
p4,
integer nin,
integer, dimension(mvsiz) ix1,
integer, dimension(mvsiz) ix2,
integer, dimension(mvsiz) ix3,
integer, dimension(mvsiz) ix4,
integer, dimension(mvsiz) nsvg,
gapv,
integer inacti,
cand_p,
integer, dimension(mvsiz) index,
stiglo,
stif,
vxi,
vyi,
vzi,
msi,
ki11,
ki12,
kj11,
kj12,
kk11,
kk12,
kl11,
kl12,
off,
scalk,
integer lrem,
integer, dimension(3) icurv,
x,
x1,
x2,
x3,
x4,
y1,
y2,
y3,
y4,
z1,
z2,
z3,
z4,
xi,
yi,
zi )

Definition at line 30 of file i7keg3.F.

45C-----------------------------------------------
46C M o d u l e s
47C-----------------------------------------------
48 USE imp_intm
49C-----------------------------------------------
50C I m p l i c i t T y p e s
51C-----------------------------------------------
52#include "implicit_f.inc"
53C-----------------------------------------------
54C G l o b a l P a r a m e t e r s
55C-----------------------------------------------
56#include "mvsiz_p.inc"
57C-----------------------------------------------
58C C o m m o n B l o c k s
59C-----------------------------------------------
60#include "scr05_c.inc"
61Ctmp+1
62#include "com01_c.inc"
63#include "impl1_c.inc"
64C-----------------------------------------------
65C D u m m y A r g u m e n t s
66C-----------------------------------------------
67 INTEGER JLT, INACTI,NIN,LREM
68 INTEGER IX1(MVSIZ), IX2(MVSIZ), IX3(MVSIZ), IX4(MVSIZ),
69 . NSVG(MVSIZ), INDEX(MVSIZ),ICURV(3)
70 my_real
71 . a(3,*), ms(*), v(3,*),
72 . stiglo,cand_p(*),fric,off(*),scalk,
73 . vxi(mvsiz),vyi(mvsiz),vzi(mvsiz),msi(mvsiz)
74 my_real
75 . nx1(mvsiz), nx2(mvsiz), nx3(mvsiz), nx4(mvsiz),
76 . ny1(mvsiz), ny2(mvsiz), ny3(mvsiz), ny4(mvsiz),
77 . nz1(mvsiz), nz2(mvsiz), nz3(mvsiz), nz4(mvsiz),
78 . lb1(mvsiz), lb2(mvsiz), lb3(mvsiz), lb4(mvsiz),
79 . lc1(mvsiz), lc2(mvsiz), lc3(mvsiz), lc4(mvsiz),
80 . p1(mvsiz), p2(mvsiz), p3(mvsiz), p4(mvsiz), stif(mvsiz),
81 . gapv(mvsiz),ki11(3,3,mvsiz),kj11(3,3,mvsiz),
82 . kk11(3,3,mvsiz),kl11(3,3,mvsiz),ki12(3,3,mvsiz),
83 . kj12(3,3,mvsiz),kk12(3,3,mvsiz),kl12(3,3,mvsiz)
84 my_real
85 . x1(mvsiz),y1(mvsiz),z1(mvsiz),
86 . x2(mvsiz),y2(mvsiz),z2(mvsiz),
87 . x3(mvsiz),y3(mvsiz),z3(mvsiz),
88 . x4(mvsiz),y4(mvsiz),z4(mvsiz),
89 . xi(mvsiz),yi(mvsiz),zi(mvsiz),x(3,*)
90C-----------------------------------------------
91C L o c a l V a r i a b l e s
92C-----------------------------------------------
93 INTEGER I, J1, J, K,IG,ISF,NN,NS,JLTF,NE
94 my_real
95 . n1(mvsiz), n2(mvsiz), n3(mvsiz), pene(mvsiz),
96 . h1(mvsiz), h2(mvsiz), h3(mvsiz), h4(mvsiz),
97 . vx(mvsiz), vy(mvsiz), vz(mvsiz), vn(mvsiz),
98 . s2,fac,facf, h0, la1, la2, la3, la4,fact(mvsiz),
99 . d1,d2,d3,d4,a1,a2,a3,a4,kn(4,mvsiz),q(3,3,mvsiz)
100 my_real
101 . prec,q11,q12,q13,q22,q23,q33,h00,vtx,vty,vtz,vt,
102 . kt1,kt2,kt3,kt4,q1,q2
103 INTEGER NA1,NA2
104 my_real
105 . a0x,a0y,a0z,rx,ry,rz,
106 . anx,any,anz,aan,aax,aay,aaz ,rr,rs,aaa ,tm,ts
107C-----------------------------------------------
108C dev021
109 IF (iresp==1) THEN
110 prec = fiveem4
111 ELSE
112 prec = em10
113 ENDIF
114C--------------------------------------------------------
115C CAS DES PAQUETS MIXTES
116C--------------------------------------------------------
117 IF (imp_int7==3) THEN
118 DO i=1,jlt
119 d1 = sqrt(p1(i))
120 p1(i) = fourth*gapv(i)
121 d2 = sqrt(p2(i))
122 p2(i) = fourth*gapv(i)
123 d3 = sqrt(p3(i))
124 p3(i) = fourth*gapv(i)
125 d4 = sqrt(p4(i))
126 p4(i) = fourth*gapv(i)
127 a1 = p1(i)/max(em20,d1)
128 a2 = p2(i)/max(em20,d2)
129 a3 = p3(i)/max(em20,d3)
130 a4 = p4(i)/max(em20,d4)
131 n1(i) = a1*nx1(i) + a2*nx2(i) + a3*nx3(i) + a4*nx4(i)
132 n2(i) = a1*ny1(i) + a2*ny2(i) + a3*ny3(i) + a4*ny4(i)
133 n3(i) = a1*nz1(i) + a2*nz2(i) + a3*nz3(i) + a4*nz4(i)
134 ENDDO
135 ELSE
136 DO i=1,jlt
137C
138 d1 = sqrt(p1(i))
139 p1(i) = max(zero, gapv(i) - d1)
140C
141 d2 = sqrt(p2(i))
142 p2(i) = max(zero, gapv(i) - d2)
143C
144 d3 = sqrt(p3(i))
145 p3(i) = max(zero, gapv(i) - d3)
146C
147 d4 = sqrt(p4(i))
148 p4(i) = max(zero, gapv(i) - d4)
149C
150 a1 = p1(i)/max(em20,d1)
151 a2 = p2(i)/max(em20,d2)
152 a3 = p3(i)/max(em20,d3)
153 a4 = p4(i)/max(em20,d4)
154 n1(i) = a1*nx1(i) + a2*nx2(i) + a3*nx3(i) + a4*nx4(i)
155 n2(i) = a1*ny1(i) + a2*ny2(i) + a3*ny3(i) + a4*ny4(i)
156 n3(i) = a1*nz1(i) + a2*nz2(i) + a3*nz3(i) + a4*nz4(i)
157 ENDDO
158 ENDIF !(IMP_INT7==3)
159C
160 DO i=1,jlt
161 IF(ix3(i)/=ix4(i))THEN
162 pene(i) = max(p1(i),p2(i),p3(i),p4(i))
163C
164 la1 = one - lb1(i) - lc1(i)
165 la2 = one - lb2(i) - lc2(i)
166 la3 = one - lb3(i) - lc3(i)
167 la4 = one - lb4(i) - lc4(i)
168C
169 h0 = fourth *
170 . (p1(i)*la1 + p2(i)*la2 + p3(i)*la3 + p4(i)*la4)
171 h1(i) = h0 + p1(i) * lb1(i) + p4(i) * lc4(i)
172 h2(i) = h0 + p2(i) * lb2(i) + p1(i) * lc1(i)
173 h3(i) = h0 + p3(i) * lb3(i) + p2(i) * lc2(i)
174 h4(i) = h0 + p4(i) * lb4(i) + p3(i) * lc3(i)
175 h00 = one/max(em20,h1(i) + h2(i) + h3(i) + h4(i))
176 h1(i) = h1(i) * h00
177 h2(i) = h2(i) * h00
178 h3(i) = h3(i) * h00
179 h4(i) = h4(i) * h00
180C
181 ELSE
182 pene(i) = p1(i)
183 n1(i) = nx1(i)
184 n2(i) = ny1(i)
185 n3(i) = nz1(i)
186 h1(i) = lb1(i)
187 h2(i) = lc1(i)
188 h3(i) = one - lb1(i) - lc1(i)
189 h4(i) = zero
190 ENDIF
191 ENDDO
192C---------------------
193C COURBURE FIXE
194C---------------------
195 IF(icurv(1)==1)THEN
196C spherique (que concave pour le moment)
197 na1 = icurv(2)
198 DO i=1,jlt
199 rr = 1.e30
200 a0x = x(1,na1)
201 a0y = x(2,na1)
202 a0z = x(3,na1)
203C
204 rx = x1(i)-a0x
205 ry = y1(i)-a0y
206 rz = z1(i)-a0z
207 rr = min(rr , rx*rx + ry*ry + rz*rz)
208 rx = x2(i)-a0x
209 ry = y2(i)-a0y
210 rz = z2(i)-a0z
211 rr = min(rr , rx*rx + ry*ry + rz*rz)
212 rx = x3(i)-a0x
213 ry = y3(i)-a0y
214 rz = z3(i)-a0z
215 rr = min(rr , rx*rx + ry*ry + rz*rz)
216 IF(ix3(i)/=ix4(i))THEN
217 rx = x4(i)-a0x
218 ry = y4(i)-a0y
219 rz = z4(i)-a0z
220 rr = min(rr , rx*rx + ry*ry + rz*rz)
221 ENDIF
222 rx = xi(i)-a0x
223 ry = yi(i)-a0y
224 rz = zi(i)-a0z
225 rs = sqrt(rx*rx + ry*ry + rz*rz)
226 rr = sqrt(rr)
227 IF(rs-rr+gapv(i)<0.)THEN
228 stif(i) = 0.
229 pene(i) = 0.
230 ELSEIF(rs-rr+gapv(i)<pene(i))THEN
231 pene(i) = rs-rr+gapv(i)
232 ENDIF
233 n1(i) = -rx
234 n2(i) = -ry
235 n3(i) = -rz
236 ENDDO
237 ELSEIF(icurv(1)==2)THEN
238C cylindrique (que concave pour le moment)
239 na1 = icurv(2)
240 na2 = icurv(3)
241 DO i=1,jlt
242 rr = 1.e30
243 a0x = x(1,na1)
244 a0y = x(2,na1)
245 a0z = x(3,na1)
246 anx = x(1,na2)-a0x
247 any = x(2,na2)-a0y
248 anz = x(3,na2)-a0z
249 aan = 1. / (anx*anx + any*any + anz*anz)
250C
251 aax = x1(i)-a0x
252 aay = y1(i)-a0y
253 aaz = z1(i)-a0z
254 aaa = (aax*anx + aay*any + aaz*anz) * aan
255 rx = aax - aaa * anx
256 ry = aay - aaa * any
257 rz = aaz - aaa * anz
258 rr = min(rr , rx*rx + ry*ry + rz*rz)
259 aax = x2(i)-a0x
260 aay = y2(i)-a0y
261 aaz = z2(i)-a0z
262 aaa = (aax*anx + aay*any + aaz*anz) * aan
263 rx = aax - aaa * anx
264 ry = aay - aaa * any
265 rz = aaz - aaa * anz
266 rr = min(rr , rx*rx + ry*ry + rz*rz)
267 aax = x3(i)-a0x
268 aay = y3(i)-a0y
269 aaz = z3(i)-a0z
270 aaa = (aax*anx + aay*any + aaz*anz) * aan
271 rx = aax - aaa * anx
272 ry = aay - aaa * any
273 rz = aaz - aaa * anz
274 rr = min(rr , rx*rx + ry*ry + rz*rz)
275 IF(ix3(i)/=ix4(i))THEN
276 aax = x4(i)-a0x
277 aay = y4(i)-a0y
278 aaz = z4(i)-a0z
279 aaa = (aax*anx + aay*any + aaz*anz) * aan
280 rx = aax - aaa * anx
281 ry = aay - aaa * any
282 rz = aaz - aaa * anz
283 rr = min(rr , rx*rx + ry*ry + rz*rz)
284 ENDIF
285 aax = xi(i)-a0x
286 aay = yi(i)-a0y
287 aaz = zi(i)-a0z
288 aaa = (aax*anx + aay*any + aaz*anz) * aan
289 rx = aax - aaa * anx
290 ry = aay - aaa * any
291 rz = aaz - aaa * anz
292 rs = sqrt(rx*rx + ry*ry + rz*rz)
293 rr = sqrt(rr)
294 IF(rs-rr+gapv(i)<0.)THEN
295 stif(i) = 0.
296 pene(i) = 0.
297 ELSEIF(rs-rr+gapv(i)<pene(i))THEN
298 pene(i) = rs-rr+gapv(i)
299 n1(i) = -rx
300 n2(i) = -ry
301 n3(i) = -rz
302 ENDIF
303 ENDDO
304 ELSEIF(icurv(1) == 3)THEN
305c print *,'ICURV=3 NOT AVAILABLE WITH IMPLICIT OPTION'
306 ENDIF
307C
308 DO i=1,jlt
309 s2 = one/max(em30,sqrt(n1(i)**2 + n2(i)**2 + n3(i)**2))
310 n1(i) = n1(i)*s2
311 n2(i) = n2(i)*s2
312 n3(i) = n3(i)*s2
313 ENDDO
314C
315 DO 500 i=1,jlt
316 vx(i) = vxi(i) - h1(i)*v(1,ix1(i)) - h2(i)*v(1,ix2(i))
317 . - h3(i)*v(1,ix3(i)) - h4(i)*v(1,ix4(i))
318 vy(i) = vyi(i) - h1(i)*v(2,ix1(i)) - h2(i)*v(2,ix2(i))
319 . - h3(i)*v(2,ix3(i)) - h4(i)*v(2,ix4(i))
320 vz(i) = vzi(i) - h1(i)*v(3,ix1(i)) - h2(i)*v(3,ix2(i))
321 . - h3(i)*v(3,ix3(i)) - h4(i)*v(3,ix4(i))
322 vn(i) = n1(i)*vx(i) + n2(i)*vy(i) + n3(i)*vz(i)
323C correction hourglass
324 h0 = -.25*(h1(i) - h2(i) + h3(i) - h4(i))
325 h0 = min(h0,h2(i),h4(i))
326 h0 = max(h0,-h1(i),-h3(i))
327 IF(ix3(i)==ix4(i))h0 = zero
328 h1(i) = h1(i) + h0
329 h2(i) = h2(i) - h0
330 h3(i) = h3(i) + h0
331 h4(i) = h4(i) - h0
332 500 CONTINUE
333C---------------------
334C PENE INITIALE
335C---------------------
336 IF(inacti==5)THEN
337 DO i=1,jlt
338 pene(i)=max(zero,pene(i)-cand_p(index(i)))
339 IF( pene(i)==zero ) stif(i) = zero
340 gapv(i)=gapv(i)-cand_p(index(i))
341 ENDDO
342 ELSE IF(inacti==6)THEN
343 DO i=1,jlt
344 cand_p(index(i))=min(cand_p(index(i)),
345 . ( (one-fiveem2)*cand_p(index(i))
346 . +fiveem2*(pene(i)+fiveem2*(gapv(i)-pene(i)))) )
347 pene(i)=max(zero,pene(i)-cand_p(index(i)))
348 IF( pene(i)==zero ) stif(i) = zero
349 gapv(i)=gapv(i)-cand_p(index(i))
350 ENDDO
351 ENDIF
352C-------------------------------------------
353c print *,'IMP_INT7=',IMP_INT7,jlt
354 IF(imp_int7>=2)THEN
355 DO i=1,jlt
356 IF(stiglo<=zero)THEN
357 stif(i) = half*stif(i)
358 ELSEIF(stif(i)/=zero)THEN
359 stif(i) = stiglo
360 ENDIF
361 ENDDO
362 ELSEIF(imp_int7==1)THEN
363 DO i=1,jlt
364 fac = gapv(i)/max( em10,( gapv(i)-pene(i) ) )
365 IF( (gapv(i)-pene(i))/gapv(i) <prec .AND.
366 . stif(i)>0. ) stif(i) = 0.
367 IF(stiglo<=zero)THEN
368 stif(i) = half*stif(i) * fac
369 ELSEIF(stif(i)/=zero)THEN
370 stif(i) = stiglo * fac
371 ENDIF
372 ENDDO
373 ELSE
374C
375 DO i=1,jlt
376 fac = gapv(i)/max( em10,( gapv(i)-pene(i) ) )
377 IF( (gapv(i)-pene(i))/gapv(i) <prec .AND.
378 . stif(i)>0. ) stif(i) = 0.
379 IF(stiglo<=zero)THEN
380 stif(i) = half*stif(i) * fac
381 ELSEIF(stif(i)/=zero)THEN
382 stif(i) = stiglo * fac
383 ENDIF
384 ENDDO
385 DO i=1,jlt
386 stif(i) = stif(i) * gapv(i) /
387 . max((gapv(i) - pene(i)),em10)
388 ENDDO
389 ENDIF
390C---------------------------------
391C ----sans frottement d'abord---
392 DO i=1,jlt
393 vtx = vx(i) -vn(i)*n1(i)
394 vty = vy(i) -vn(i)*n2(i)
395 vtz = vz(i) -vn(i)*n3(i)
396 vt = vtx*vtx+vty*vty+vtz*vtz
397 IF (vt>em20) THEN
398 s2=one/sqrt(vt)
399 q(1,1,i)=vtx*s2
400 q(1,2,i)=vty*s2
401 q(1,3,i)=vtz*s2
402 q(3,1,i)=n1(i)
403 q(3,2,i)=n2(i)
404 q(3,3,i)=n3(i)
405 q(2,1,i)=q(3,2,i)*q(1,3,i)-q(3,3,i)*q(1,2,i)
406 q(2,2,i)=q(3,3,i)*q(1,1,i)-q(3,1,i)*q(1,3,i)
407 q(2,3,i)=q(3,1,i)*q(1,2,i)-q(3,2,i)*q(1,1,i)
408 fact(i)=fric
409 ELSE
410 fact(i)=zero
411 ENDIF
412 ENDDO
413 IF (scalk<0) THEN
414 isf=1
415 ELSE
416 isf=0
417 ENDIF
418 facf=abs(scalk)
419 IF (isf==1) THEN
420 DO i=1,jlt
421 IF (vn(i)>zero) THEN
422 fac=stif(i)/facf
423 ELSEIF (vn(i)<zero) THEN
424 fac=stif(i)*facf
425 ELSE
426 fac=stif(i)
427 ENDIF
428 kn(1,i)=fac*h1(i)
429 kn(2,i)=fac*h2(i)
430 kn(3,i)=fac*h3(i)
431 kn(4,i)=fac*h4(i)
432 fact(i)=fac*fact(i)
433 ENDDO
434 ELSE
435 DO i=1,jlt
436 fac=stif(i)*facf
437 kn(1,i)=fac*h1(i)
438 kn(2,i)=fac*h2(i)
439 kn(3,i)=fac*h3(i)
440 kn(4,i)=fac*h4(i)
441 fact(i)=fac*fact(i)
442 ENDDO
443 ENDIF
444 DO i=1,jlt
445 q11=n1(i)*n1(i)
446 q12=n1(i)*n2(i)
447 q13=n1(i)*n3(i)
448 q22=n2(i)*n2(i)
449 q23=n2(i)*n3(i)
450 q33=n3(i)*n3(i)
451 ki11(1,1,i)=kn(1,i)*q11
452 ki11(1,2,i)=kn(1,i)*q12
453 ki11(1,3,i)=kn(1,i)*q13
454 ki11(2,2,i)=kn(1,i)*q22
455 ki11(2,3,i)=kn(1,i)*q23
456 ki11(3,3,i)=kn(1,i)*q33
457 kj11(1,1,i)=kn(2,i)*q11
458 kj11(1,2,i)=kn(2,i)*q12
459 kj11(1,3,i)=kn(2,i)*q13
460 kj11(2,2,i)=kn(2,i)*q22
461 kj11(2,3,i)=kn(2,i)*q23
462 kj11(3,3,i)=kn(2,i)*q33
463 kk11(1,1,i)=kn(3,i)*q11
464 kk11(1,2,i)=kn(3,i)*q12
465 kk11(1,3,i)=kn(3,i)*q13
466 kk11(2,2,i)=kn(3,i)*q22
467 kk11(2,3,i)=kn(3,i)*q23
468 kk11(3,3,i)=kn(3,i)*q33
469 kl11(1,1,i)=kn(4,i)*q11
470 kl11(1,2,i)=kn(4,i)*q12
471 kl11(1,3,i)=kn(4,i)*q13
472 kl11(2,2,i)=kn(4,i)*q22
473 kl11(2,3,i)=kn(4,i)*q23
474 kl11(3,3,i)=kn(4,i)*q33
475 ENDDO
476C ----avec frottement ---
477 DO j=1,3
478 DO k=j,3
479 DO i=1,jlt
480 IF (fact(i)>zero) THEN
481 q1 =q(1,j,i)*q(1,k,i)
482 q2 =q(2,j,i)*q(2,k,i)
483 fac=fact(i)*(q1+q2)
484 kt1=fac*h1(i)
485 ki11(j,k,i)=ki11(j,k,i)+kt1
486 kt2=fac*h2(i)
487 kj11(j,k,i)=kj11(j,k,i)+kt2
488 kt3=fac*h3(i)
489 kk11(j,k,i)=kk11(j,k,i)+kt3
490 kt4=fac*h4(i)
491 kl11(j,k,i)=kl11(j,k,i)+kt4
492 ENDIF
493 ENDDO
494 ENDDO
495 ENDDO
496C
497 DO j=1,3
498 DO k=j,3
499 DO i=1,jlt
500 ki12(j,k,i)=-ki11(j,k,i)
501 kj12(j,k,i)=-kj11(j,k,i)
502 kk12(j,k,i)=-kk11(j,k,i)
503 kl12(j,k,i)=-kl11(j,k,i)
504 ENDDO
505 ENDDO
506 ENDDO
507 DO j=1,3
508 DO k=j+1,3
509 DO i=1,jlt
510 ki12(k,j,i)=-ki11(j,k,i)
511 kj12(k,j,i)=-kj11(j,k,i)
512 kk12(k,j,i)=-kk11(j,k,i)
513 kl12(k,j,i)=-kl11(j,k,i)
514 ENDDO
515 ENDDO
516 ENDDO
517C
518 DO i=1,jlt
519 off(i)=one
520 ENDDO
521C
522 IF (intp_d <= 0 .AND. nspmd > 1) THEN
523 jltf = 0
524 DO i=1,jlt
525 IF(nsvg(i)<0) THEN
526 nn=-nsvg(i)
527 jltf = jltf + 1
528 ne=shf_int(nin) + jltf +lrem
529 ns=ind_int(nin)%P(nn)
530 stifs(ne)=stif(i)
531 h_e(1,ne)=h1(i)
532 h_e(2,ne)=h2(i)
533 h_e(3,ne)=h3(i)
534 h_e(4,ne)=h4(i)
535 n_e(1,ne)=n1(i)
536 n_e(2,ne)=n2(i)
537 n_e(3,ne)=n3(i)
538 ENDIF
539 ENDDO
540 ENDIF
541C
542 RETURN
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
imp_intm::shf_int
integer, dimension(:), allocatable shf_int
Definition imp_intm.F:136
imp_intm::intp_d
integer intp_d
Definition imp_intm.F:173
imp_intm::ind_int
type(int_pointer2), dimension(:), allocatable ind_int
Definition imp_intm.F:133

◆ i7kfor3()

subroutine i7kfor3 ( integer jlt,
a,
v,
ms,
fric,
nx1,
nx2,
nx3,
nx4,
ny1,
ny2,
ny3,
ny4,
nz1,
nz2,
nz3,
nz4,
lb1,
lb2,
lb3,
lb4,
lc1,
lc2,
lc3,
lc4,
p1,
p2,
p3,
p4,
integer nin,
integer, dimension(mvsiz) ix1,
integer, dimension(mvsiz) ix2,
integer, dimension(mvsiz) ix3,
integer, dimension(mvsiz) ix4,
integer, dimension(mvsiz) nsvg,
gapv,
integer inacti,
cand_p,
integer, dimension(mvsiz) index,
stiglo,
stif,
vxi,
vyi,
vzi,
msi,
x,
integer, dimension(4,*) irect,
integer, dimension(mvsiz) ce_loc,
integer mfrot,
integer ifq,
frot_p,
cand_fx,
cand_fy,
cand_fz,
alpha0,
dxi,
dyi,
dzi,
d,
scalk,
x1,
x2,
x3,
x4,
y1,
y2,
y3,
y4,
z1,
z2,
z3,
z4,
xi,
yi,
zi,
icurv )

Definition at line 710 of file i7keg3.F.

726C-----------------------------------------------
727C M o d u l e s
728C-----------------------------------------------
729 USE imp_intm
730C-----------------------------------------------
731C I m p l i c i t T y p e s
732C-----------------------------------------------
733#include "implicit_f.inc"
734#include "comlock.inc"
735C-----------------------------------------------
736C G l o b a l P a r a m e t e r s
737C-----------------------------------------------
738#include "mvsiz_p.inc"
739C-----------------------------------------------
740C C o m m o n B l o c k s
741C-----------------------------------------------
742#include "com08_c.inc"
743#include "scr05_c.inc"
744#include "impl1_c.inc"
745C-----------------------------------------------
746C D u m m y A r g u m e n t s
747C-----------------------------------------------
748 INTEGER JLT,INACTI,NIN,MFROT, IFQ, NOINT,IRECT(4,*)
749 INTEGER IX1(MVSIZ), IX2(MVSIZ), IX3(MVSIZ), IX4(MVSIZ),
750 . CE_LOC(MVSIZ),INDEX(MVSIZ),NSVG(MVSIZ)
751 my_real
752 . stiglo,cand_p(*),frot_p(*), x(3,*),a(3,*), ms(*), v(3,*),
753 . cand_fx(*),cand_fy(*),cand_fz(*),alpha0,fric,scalk,icurv(3)
754 my_real
755 . nx1(mvsiz), nx2(mvsiz), nx3(mvsiz), nx4(mvsiz),
756 . ny1(mvsiz), ny2(mvsiz), ny3(mvsiz), ny4(mvsiz),
757 . nz1(mvsiz), nz2(mvsiz), nz3(mvsiz), nz4(mvsiz),
758 . lb1(mvsiz), lb2(mvsiz), lb3(mvsiz), lb4(mvsiz),
759 . lc1(mvsiz), lc2(mvsiz), lc3(mvsiz), lc4(mvsiz),
760 . p1(mvsiz), p2(mvsiz), p3(mvsiz), p4(mvsiz), stif(mvsiz),
761 . gapv(mvsiz),vxi(mvsiz),vyi(mvsiz),vzi(mvsiz),msi(mvsiz),
762 . dxi(mvsiz),dyi(mvsiz),dzi(mvsiz),d(3,*)
763 my_real
764 . x1(mvsiz),y1(mvsiz),z1(mvsiz),
765 . x2(mvsiz),y2(mvsiz),z2(mvsiz),
766 . x3(mvsiz),y3(mvsiz),z3(mvsiz),
767 . x4(mvsiz),y4(mvsiz),z4(mvsiz),
768 . xi(mvsiz),yi(mvsiz),zi(mvsiz)
769C-----------------------------------------------
770C L o c a l V a r i a b l e s
771C-----------------------------------------------
772 INTEGER I, J1, IG, J, JG , K0,NBINTER,K1S,K,IL,IE, NN, NI,NS
773 my_real
774 . fxi(mvsiz), fyi(mvsiz), fzi(mvsiz), fni(mvsiz),
775 . fxt(mvsiz),fyt(mvsiz),fzt(mvsiz),stif0(mvsiz),
776 . n1(mvsiz), n2(mvsiz), n3(mvsiz), pene(mvsiz),
777 . h1(mvsiz), h2(mvsiz), h3(mvsiz), h4(mvsiz),
778 . vx(mvsiz), vy(mvsiz), vz(mvsiz), vn(mvsiz),xmu(mvsiz),
779 . dx(mvsiz), dy(mvsiz), dz(mvsiz), dn0(mvsiz),
780 .
781 . vnx, vny, vnz, aa, crit,s2,dist,rdist,
782 . v2, fm2, fac,ff,alphi,alpha,beta,
783 . fx, fy, fz, f2, mas2, m2sk, dti,ft,fn,fmax,ftn,
784 . h0, la1, la2, la3, la4,d1,d2,d3,d4,a1,a2,a3,a4, ff0,
785 . vv,ax1,ax2,ay1,ay2,az1,az2,ax,ay,az,t1,t2,t3,dxt,
786 . area,p,vv1,vv2,v21,dmu, dti2,h00 ,a0x,a0y,a0z,rx,ry,rz,
787 . anx,any,anz,aan,aax,aay,aaz ,rr,rs,aaa ,gap2 ,pene2,prec
788 INTEGER NA1,NA2
789C-----------------------------------------------
790 IF (iresp==1) THEN
791 prec = fiveem4
792 ELSE
793 prec = em10
794 ENDIF
795C---------------------
796C PENE INITIALE
797C---------------------
798 IF(inacti==5.OR.inacti==6)THEN
799 DO i=1,jlt
800 IF(stif(i)==zero)THEN
801 cand_p(index(i))=0
802 ENDIF
803 ENDDO
804 ENDIF
805C--------------------------------------------------------
806C actualise stif
807C--------------------------------------------------------
808 DO i=1,jlt
809 gap2=gapv(i)*gapv(i)
810C
811 d1 = max(zero, gap2 - p1(i))
812 d2 = max(zero, gap2 - p2(i))
813 d3 = max(zero, gap2 - p3(i))
814 d4 = max(zero, gap2 - p4(i))
815 pene2 = max(d1,d2,d3,d4)
816 IF (pene2<=zero) stif(i) = zero
817 ENDDO
818C--------------------------------------------------------
819C CAS DES PAQUETS MIXTES
820C--------------------------------------------------------
821 DO i=1,jlt
822 d1 = sqrt(p1(i))
823 p1(i) = max(zero, gapv(i) - d1)
824C
825 d2 = sqrt(p2(i))
826 p2(i) = max(zero, gapv(i) - d2)
827C
828 d3 = sqrt(p3(i))
829 p3(i) = max(zero, gapv(i) - d3)
830C
831 d4 = sqrt(p4(i))
832 p4(i) = max(zero, gapv(i) - d4)
833C
834 a1 = p1(i)/max(em20,d1)
835 a2 = p2(i)/max(em20,d2)
836 a3 = p3(i)/max(em20,d3)
837 a4 = p4(i)/max(em20,d4)
838 n1(i) = a1*nx1(i) + a2*nx2(i) + a3*nx3(i) + a4*nx4(i)
839 n2(i) = a1*ny1(i) + a2*ny2(i) + a3*ny3(i) + a4*ny4(i)
840 n3(i) = a1*nz1(i) + a2*nz2(i) + a3*nz3(i) + a4*nz4(i)
841 ENDDO
842C
843 DO i=1,jlt
844 IF(ix3(i)/=ix4(i))THEN
845 pene(i) = max(p1(i),p2(i),p3(i),p4(i))
846C
847 la1 = one - lb1(i) - lc1(i)
848 la2 = one - lb2(i) - lc2(i)
849 la3 = one - lb3(i) - lc3(i)
850 la4 = one - lb4(i) - lc4(i)
851C
852 h0 = fourth *
853 . (p1(i)*la1 + p2(i)*la2 + p3(i)*la3 + p4(i)*la4)
854 h1(i) = h0 + p1(i) * lb1(i) + p4(i) * lc4(i)
855 h2(i) = h0 + p2(i) * lb2(i) + p1(i) * lc1(i)
856 h3(i) = h0 + p3(i) * lb3(i) + p2(i) * lc2(i)
857 h4(i) = h0 + p4(i) * lb4(i) + p3(i) * lc3(i)
858 h00 = one/max(em20,h1(i) + h2(i) + h3(i) + h4(i))
859 h1(i) = h1(i) * h00
860 h2(i) = h2(i) * h00
861 h3(i) = h3(i) * h00
862 h4(i) = h4(i) * h00
863C
864 ELSE
865 pene(i) = p1(i)
866 n1(i) = nx1(i)
867 n2(i) = ny1(i)
868 n3(i) = nz1(i)
869 h1(i) = lb1(i)
870 h2(i) = lc1(i)
871 h3(i) = one - lb1(i) - lc1(i)
872 h4(i) = zero
873 ENDIF
874 ENDDO
875C---------------------
876C COURBURE FIXE
877C---------------------
878 IF(icurv(1)==1)THEN
879C spherique (que concave pour le moment)
880 na1 = icurv(2)
881 DO i=1,jlt
882 rr = 1.e30
883 a0x = x(1,na1)
884 a0y = x(2,na1)
885 a0z = x(3,na1)
886C
887 rx = x1(i)-a0x
888 ry = y1(i)-a0y
889 rz = z1(i)-a0z
890 rr = min(rr , rx*rx + ry*ry + rz*rz)
891 rx = x2(i)-a0x
892 ry = y2(i)-a0y
893 rz = z2(i)-a0z
894 rr = min(rr , rx*rx + ry*ry + rz*rz)
895 rx = x3(i)-a0x
896 ry = y3(i)-a0y
897 rz = z3(i)-a0z
898 rr = min(rr , rx*rx + ry*ry + rz*rz)
899 IF(ix3(i)/=ix4(i))THEN
900 rx = x4(i)-a0x
901 ry = y4(i)-a0y
902 rz = z4(i)-a0z
903 rr = min(rr , rx*rx + ry*ry + rz*rz)
904 ENDIF
905 rx = xi(i)-a0x
906 ry = yi(i)-a0y
907 rz = zi(i)-a0z
908 rs = sqrt(rx*rx + ry*ry + rz*rz)
909 rr = sqrt(rr)
910 IF(rs-rr+gapv(i)<0.)THEN
911 stif(i) = 0.
912 pene(i) = 0.
913 ELSEIF(rs-rr+gapv(i)<pene(i))THEN
914 pene(i) = rs-rr+gapv(i)
915 ENDIF
916 n1(i) = -rx
917 n2(i) = -ry
918 n3(i) = -rz
919 ENDDO
920 ELSEIF(icurv(1)==2)THEN
921C cylindrique (que concave pour le moment)
922 na1 = icurv(2)
923 na2 = icurv(3)
924 DO i=1,jlt
925 rr = 1.e30
926 a0x = x(1,na1)
927 a0y = x(2,na1)
928 a0z = x(3,na1)
929 anx = x(1,na2)-a0x
930 any = x(2,na2)-a0y
931 anz = x(3,na2)-a0z
932 aan = 1. / (anx*anx + any*any + anz*anz)
933C
934 aax = x1(i)-a0x
935 aay = y1(i)-a0y
936 aaz = z1(i)-a0z
937 aaa = (aax*anx + aay*any + aaz*anz) * aan
938 rx = aax - aaa * anx
939 ry = aay - aaa * any
940 rz = aaz - aaa * anz
941 rr = min(rr , rx*rx + ry*ry + rz*rz)
942 aax = x2(i)-a0x
943 aay = y2(i)-a0y
944 aaz = z2(i)-a0z
945 aaa = (aax*anx + aay*any + aaz*anz) * aan
946 rx = aax - aaa * anx
947 ry = aay - aaa * any
948 rz = aaz - aaa * anz
949 rr = min(rr , rx*rx + ry*ry + rz*rz)
950 aax = x3(i)-a0x
951 aay = y3(i)-a0y
952 aaz = z3(i)-a0z
953 aaa = (aax*anx + aay*any + aaz*anz) * aan
954 rx = aax - aaa * anx
955 ry = aay - aaa * any
956 rz = aaz - aaa * anz
957 rr = min(rr , rx*rx + ry*ry + rz*rz)
958 IF(ix3(i)/=ix4(i))THEN
959 aax = x4(i)-a0x
960 aay = y4(i)-a0y
961 aaz = z4(i)-a0z
962 aaa = (aax*anx + aay*any + aaz*anz) * aan
963 rx = aax - aaa * anx
964 ry = aay - aaa * any
965 rz = aaz - aaa * anz
966 rr = min(rr , rx*rx + ry*ry + rz*rz)
967 ENDIF
968 aax = xi(i)-a0x
969 aay = yi(i)-a0y
970 aaz = zi(i)-a0z
971 aaa = (aax*anx + aay*any + aaz*anz) * aan
972 rx = aax - aaa * anx
973 ry = aay - aaa * any
974 rz = aaz - aaa * anz
975 rs = sqrt(rx*rx + ry*ry + rz*rz)
976 rr = sqrt(rr)
977 IF(rs-rr+gapv(i)<0.)THEN
978 stif(i) = 0.
979 pene(i) = 0.
980 ELSEIF(rs-rr+gapv(i)<pene(i))THEN
981 pene(i) = rs-rr+gapv(i)
982 n1(i) = -rx
983 n2(i) = -ry
984 n3(i) = -rz
985 ENDIF
986 ENDDO
987 ELSEIF(icurv(1) == 3)THEN
988c print *,'ICURV=3 NOT AVAILABLE WITH IMPLICIT OPTION'
989 ENDIF
990C
991 DO i=1,jlt
992 s2 = one/max(em30,sqrt(n1(i)**2 + n2(i)**2 + n3(i)**2))
993 n1(i) = n1(i)*s2
994 n2(i) = n2(i)*s2
995 n3(i) = n3(i)*s2
996 ENDDO
997C
998 DO i=1,jlt
999 dx(i) = dxi(i) - h1(i)*d(1,ix1(i)) - h2(i)*d(1,ix2(i))
1000 . - h3(i)*d(1,ix3(i)) - h4(i)*d(1,ix4(i))
1001 dy(i) = dyi(i) - h1(i)*d(2,ix1(i)) - h2(i)*d(2,ix2(i))
1002 . - h3(i)*d(2,ix3(i)) - h4(i)*d(2,ix4(i))
1003 dz(i) = dzi(i) - h1(i)*d(3,ix1(i)) - h2(i)*d(3,ix2(i))
1004 . - h3(i)*d(3,ix3(i)) - h4(i)*d(3,ix4(i))
1005 dn0(i) = n1(i)*dx(i) + n2(i)*dy(i) + n3(i)*dz(i)
1006 ENDDO
1007C
1008 DO 500 i=1,jlt
1009 vx(i) = vxi(i) - h1(i)*v(1,ix1(i)) - h2(i)*v(1,ix2(i))
1010 . - h3(i)*v(1,ix3(i)) - h4(i)*v(1,ix4(i))
1011 vy(i) = vyi(i) - h1(i)*v(2,ix1(i)) - h2(i)*v(2,ix2(i))
1012 . - h3(i)*v(2,ix3(i)) - h4(i)*v(2,ix4(i))
1013 vz(i) = vzi(i) - h1(i)*v(3,ix1(i)) - h2(i)*v(3,ix2(i))
1014 . - h3(i)*v(3,ix3(i)) - h4(i)*v(3,ix4(i))
1015 vn(i) = n1(i)*vx(i) + n2(i)*vy(i) + n3(i)*vz(i)
1016C correction hourglass
1017 h0 = -.25*(h1(i) - h2(i) + h3(i) - h4(i))
1018 h0 = min(h0,h2(i),h4(i))
1019 h0 = max(h0,-h1(i),-h3(i))
1020 IF(ix3(i)==ix4(i))h0 = zero
1021 h1(i) = h1(i) + h0
1022 h2(i) = h2(i) - h0
1023 h3(i) = h3(i) + h0
1024 h4(i) = h4(i) - h0
1025 500 CONTINUE
1026C---------------------
1027C PENE INITIALE
1028C---------------------
1029 IF(inacti==5)THEN
1030 DO i=1,jlt
1031C REDUCTION DE LA PENE INITIALE
1032 cand_p(index(i))=min(cand_p(index(i)),
1033 . ((one-fiveem2)*cand_p(index(i))+fiveem2*pene(i)) )
1034C SOUSTRACTION DE LA PENE INITIALE A LA PENE ET AU GAP
1035 pene(i)=max(zero,pene(i)-cand_p(index(i)))
1036 IF( pene(i)==zero ) stif(i) = zero
1037 gapv(i)=gapv(i)-cand_p(index(i))
1038 ENDDO
1039 ELSE IF(inacti==6)THEN
1040 DO i=1,jlt
1041C REDUCTION DE LA PENE INITIALE
1042C CAND_P(INDEX(I))=MIN(CAND_P(INDEX(I)),PENE(I))
1043 cand_p(index(i))=min(cand_p(index(i)),
1044 . ( (one-fiveem2)*cand_p(index(i))
1045 . +fiveem2*(pene(i)+fiveem2*(gapv(i)-pene(i)))) )
1046C SOUSTRACTION DE LA PENE INITIALE A LA PENE ET AU GAP
1047 pene(i)=max(zero,pene(i)-cand_p(index(i)))
1048 IF( pene(i)==zero ) stif(i) = zero
1049 gapv(i)=gapv(i)-cand_p(index(i))
1050 ENDDO
1051 ENDIF
1052C
1053 DO i=1,jlt
1054 stif0(i) = stif(i)
1055 ENDDO
1056C-------------------------------------------
1057 IF(imp_int7>=2)THEN
1058 DO i=1,jlt
1059 IF(stiglo<=zero)THEN
1060 stif(i) = half*stif(i)
1061 ELSEIF(stif(i)/=zero)THEN
1062 stif(i) = stiglo
1063 ENDIF
1064 ENDDO
1065 ELSEIF(imp_int7==1)THEN
1066 DO i=1,jlt
1067 fac = gapv(i)/max( em10,( gapv(i)-pene(i) ) )
1068 IF( (gapv(i)-pene(i))/gapv(i) <prec .AND.
1069 . stif(i)>0. ) stif(i) = 0.
1070 IF(stiglo<=zero)THEN
1071 stif(i) = half*stif(i) * fac
1072 ELSEIF(stif(i)/=zero)THEN
1073 stif(i) = stiglo * fac
1074 ENDIF
1075 ENDDO
1076 ELSE
1077 DO i=1,jlt
1078 fac = gapv(i)/max( em10,( gapv(i)-pene(i) ) )
1079 IF( (gapv(i)-pene(i))/gapv(i) <prec .AND.
1080 . stif(i)>0. ) stif(i) = 0.
1081 IF(stiglo<=zero)THEN
1082 stif(i) = half*stif(i) * fac
1083 ELSEIF(stif(i)/=zero)THEN
1084 stif(i) = stiglo * fac
1085 ENDIF
1086 ENDDO
1087 DO i=1,jlt
1088 stif(i) = stif(i) * gapv(i) /
1089 . max((gapv(i) - pene(i)),em10)
1090 ENDDO
1091 ENDIF ! (IMP_INT7>=2)
1092C
1093 fac = abs(scalk)
1094 DO i=1,jlt
1095 stif(i)=stif(i)*fac
1096 fni(i)= -stif(i) * dn0(i)
1097 fxi(i)=n1(i)*fni(i)
1098 fyi(i)=n2(i)*fni(i)
1099 fzi(i)=n3(i)*fni(i)
1100 ENDDO
1101C------------------
1102C TANGENT FORCE CALCULATION
1103C------------------
1104C---------------------------------
1105C NEW FRICTION MODELS
1106C---------------------------------
1107 IF (mfrot==0) THEN
1108C--- Coulomb friction
1109 DO i=1,jlt
1110 xmu(i) = fric
1111 ENDDO
1112 ELSEIF (mfrot==1) THEN
1113C--- Viscous friction
1114 DO i=1,jlt
1115 ie=ce_loc(i)
1116 aa = n1(i)*vx(i) + n2(i)*vy(i) + n3(i)*vz(i)
1117 v2 = (vx(i) - n1(i)*aa)**2
1118 . + (vy(i) - n2(i)*aa)**2
1119 . + (vz(i) - n3(i)*aa)**2
1120 vv = sqrt(max(em30,v2))
1121 ax1 = x(1,irect(3,ie)) - x(1,irect(1,ie))
1122 ay1 = x(2,irect(3,ie)) - x(2,irect(1,ie))
1123 az1 = x(3,irect(3,ie)) - x(3,irect(1,ie))
1124 ax2 = x(1,irect(4,ie)) - x(1,irect(2,ie))
1125 ay2 = x(2,irect(4,ie)) - x(2,irect(2,ie))
1126 az2 = x(3,irect(4,ie)) - x(3,irect(2,ie))
1127 ax = ay1*az2 - az1*ay2
1128 ay = az1*ax2 - ax1*az2
1129 az = ax1*ay2 - ay1*ax2
1130 area = half*sqrt(ax*ax+ay*ay+az*az)
1131 p = -fni(i)/area
1132 xmu(i) = fric + (frot_p(1) + frot_p(4)*p ) * p
1133 . +(frot_p(2) + frot_p(3)*p) * vv + frot_p(5)*v2
1134 ENDDO
1135 ELSEIF(mfrot==2)THEN
1136C--- Loi Darmstad
1137 DO i=1,jlt
1138 ie=ce_loc(i)
1139 aa = n1(i)*vx(i) + n2(i)*vy(i) + n3(i)*vz(i)
1140 v2 = (vx(i) - n1(i)*aa)**2
1141 . + (vy(i) - n2(i)*aa)**2
1142 . + (vz(i) - n3(i)*aa)**2
1143 vv = sqrt(max(em30,v2))
1144 ax1 = x(1,irect(3,ie)) - x(1,irect(1,ie))
1145 ay1 = x(2,irect(3,ie)) - x(2,irect(1,ie))
1146 az1 = x(3,irect(3,ie)) - x(3,irect(1,ie))
1147 ax2 = x(1,irect(4,ie)) - x(1,irect(2,ie))
1148 ay2 = x(2,irect(4,ie)) - x(2,irect(2,ie))
1149 az2 = x(3,irect(4,ie)) - x(3,irect(2,ie))
1150 ax = ay1*az2 - az1*ay2
1151 ay = az1*ax2 - ax1*az2
1152 az = ax1*ay2 - ay1*ax2
1153 area = half*sqrt(ax*ax+ay*ay+az*az)
1154 p = -fni(i)/area
1155 xmu(i) = frot_p(1)*exp(frot_p(2)*vv)*p*p
1156 . + frot_p(3)*exp(frot_p(4)*vv)*p
1157 . + frot_p(5)*exp(frot_p(6)*vv)
1158 ENDDO
1159 ELSEIF (mfrot==3) THEN
1160C--- Loi Renard
1161 DO i=1,jlt
1162 aa = n1(i)*vx(i) + n2(i)*vy(i) + n3(i)*vz(i)
1163 v2 = (vx(i) - n1(i)*aa)**2
1164 . + (vy(i) - n2(i)*aa)**2
1165 . + (vz(i) - n3(i)*aa)**2
1166 vv = sqrt(max(em30,v2))
1167 IF(vv>=0.AND.vv<=frot_p(5)) THEN
1168 dmu = frot_p(3)-frot_p(1)
1169 vv1 = vv / frot_p(5)
1170 xmu(i) = frot_p(1)+ dmu*vv1*(two-vv1)
1171 ELSEIF(vv>frot_p(5).AND.vv<frot_p(6)) THEN
1172 dmu = frot_p(4)-frot_p(3)
1173 vv1 = (vv - frot_p(5))/(frot_p(6)-frot_p(5))
1174 xmu(i) = frot_p(3)+ dmu * (three-two*vv1)*vv1**2
1175 ELSE
1176 dmu = frot_p(2)-frot_p(4)
1177 vv2 = (vv - frot_p(6))**2
1178 xmu(i) = frot_p(2) - dmu / (one + dmu*vv2)
1179 ENDIF
1180 ENDDO
1181 ENDIF
1182C
1183C---------------------------------
1184C INCREMENTAL (STIFFNESS) FORMULATION
1185C---------------------------------
1186 IF (ifq>=10.AND.scalk<zero) THEN
1187 IF (ifq==13) THEN
1188 alpha = max(one,alpha0*dt12)
1189 ELSE
1190 alpha = alpha0
1191 ENDIF
1192 DO i=1,jlt
1193 fx = stif0(i)*dx(i)
1194 fy = stif0(i)*dy(i)
1195 fz = stif0(i)*dz(i)
1196 fx = cand_fx(index(i)) + alpha*fx
1197 fy = cand_fy(index(i)) + alpha*fy
1198 fz = cand_fz(index(i)) + alpha*fz
1199 ftn = fx*n1(i) + fy*n2(i) + fz*n3(i)
1200 fx = fx - ftn*n1(i)
1201 fy = fy - ftn*n2(i)
1202 fz = fz - ftn*n3(i)
1203 ft = fx*fx + fy*fy + fz*fz
1204 ft = max(ft,em30)
1205 fn = fxi(i)**2+fyi(i)**2+fzi(i)**2
1206 beta = min(one,xmu(i)*sqrt(fn/ft))
1207 fxt(i) = fx * beta
1208 fyt(i) = fy * beta
1209 fzt(i) = fz * beta
1210C------- total force
1211 fxi(i)=fxi(i) + fxt(i)
1212 fyi(i)=fyi(i) + fyt(i)
1213 fzi(i)=fzi(i) + fzt(i)
1214 ENDDO
1215 ELSE
1216C---------------------------------
1217C TOTAL (VISCOUS) FORMULATION / NO FRICTION FILTERING
1218C---------------------------------
1219 IF (fric>0) THEN
1220 DO i=1,jlt
1221 vnx = n1(i)*dn0(i)
1222 vny = n2(i)*dn0(i)
1223 vnz = n3(i)*dn0(i)
1224 vx(i) = dx(i) - vnx
1225 vy(i) = dy(i) - vny
1226 vz(i) = dz(i) - vnz
1227 v2 = vx(i)**2 + vy(i)**2 + vz(i)**2
1228 dxt = sqrt(v2)
1229 aa = dxt/max(em30,v2)
1230 t1 = vx(i)*aa
1231 t2 = vy(i)*aa
1232 t3 = vz(i)*aa
1233 ftn = -xmu(i)*stif(i) * dxt
1234 fxt(i) = ftn * t1
1235 fyt(i) = ftn * t2
1236 fzt(i) = ftn * t3
1237C------- total force
1238 fxi(i)=fxi(i) + fxt(i)
1239 fyi(i)=fyi(i) + fyt(i)
1240 fzi(i)=fzi(i) + fzt(i)
1241 ENDDO
1242 ENDIF
1243 ENDIF
1244C
1245C--------main part-------
1246c
1247 DO i=1,jlt
1248 fx=fxi(i)
1249 fy=fyi(i)
1250 fz=fzi(i)
1251 a(1,ix1(i))=a(1,ix1(i))+fx*h1(i)
1252 a(1,ix2(i))=a(1,ix2(i))+fx*h2(i)
1253 a(1,ix3(i))=a(1,ix3(i))+fx*h3(i)
1254 a(1,ix4(i))=a(1,ix4(i))+fx*h4(i)
1255 a(2,ix1(i))=a(2,ix1(i))+fy*h1(i)
1256 a(2,ix2(i))=a(2,ix2(i))+fy*h2(i)
1257 a(2,ix3(i))=a(2,ix3(i))+fy*h3(i)
1258 a(2,ix4(i))=a(2,ix4(i))+fy*h4(i)
1259 a(3,ix1(i))=a(3,ix1(i))+fz*h1(i)
1260 a(3,ix2(i))=a(3,ix2(i))+fz*h2(i)
1261 a(3,ix3(i))=a(3,ix3(i))+fz*h3(i)
1262 a(3,ix4(i))=a(3,ix4(i))+fz*h4(i)
1263 ENDDO
1264C--------secnd part-------
1265
1266 DO i=1,jlt
1267 ig=nsvg(i)
1268 IF(ig>0)THEN
1269 a(1,ig)=a(1,ig)-fxi(i)
1270 a(2,ig)=a(2,ig)-fyi(i)
1271 a(3,ig)=a(3,ig)-fzi(i)
1272 ELSE
1273 nn=-ig
1274 ns=ind_int(nin)%P(nn)
1275 ffi(1,ns)=ffi(1,ns)-fxi(i)
1276 ffi(2,ns)=ffi(2,ns)-fyi(i)
1277 ffi(3,ns)=ffi(3,ns)-fzi(i)
1278 ENDIF
1279 ENDDO
1280C
1281 RETURN
alpha
#define alpha
Definition eval.h:35
area
subroutine area(d1, x, x2, y, y2, eint, stif0)
Definition hm_read_prop32.F:567