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"
#include "comlock.inc"
#include "com08_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, 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, ISF, 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)
594C-----------------------------------------------
595C
596 DO i=1,jlt
597 vx(i) = vxi(i) - h1(i)*v(1,ix1(i)) - h2(i)*v(1,ix2(i))
598 . - h3(i)*v(1,ix3(i)) - h4(i)*v(1,ix4(i))
599 vy(i) = vyi(i) - h1(i)*v(2,ix1(i)) - h2(i)*v(2,ix2(i))
600 . - h3(i)*v(2,ix3(i)) - h4(i)*v(2,ix4(i))
601 vz(i) = vzi(i) - h1(i)*v(3,ix1(i)) - h2(i)*v(3,ix2(i))
602 . - h3(i)*v(3,ix3(i)) - h4(i)*v(3,ix4(i))
603 vn(i) = n1(i)*vx(i) + n2(i)*vy(i) + n3(i)*vz(i)
604 dx(i) = dxi(i) - h1(i)*d(1,ix1(i)) - h2(i)*d(1,ix2(i))
605 . - h3(i)*d(1,ix3(i)) - h4(i)*d(1,ix4(i))
606 dy(i) = dyi(i) - h1(i)*d(2,ix1(i)) - h2(i)*d(2,ix2(i))
607 . - h3(i)*d(2,ix3(i)) - h4(i)*d(2,ix4(i))
608 dz(i) = dzi(i) - h1(i)*d(3,ix1(i)) - h2(i)*d(3,ix2(i))
609 . - h3(i)*d(3,ix3(i)) - h4(i)*d(3,ix4(i))
610 dn(i) = n1(i)*dx(i) + n2(i)*dy(i) + n3(i)*dz(i)
611 dni(i) = n1(i)*dxi(i) + n2(i)*dyi(i) + n3(i)*dzi(i)
612 ENDDO
613C-------------------------------------------
614 DO i=1,jlt
615 vtx = vx(i) -vn(i)*n1(i)
616 vty = vy(i) -vn(i)*n2(i)
617 vtz = vz(i) -vn(i)*n3(i)
618 vt = vtx*vtx+vty*vty+vtz*vtz
619 IF (vt>em20) THEN
620 s2=one/sqrt(vt)
621 t1(i)=vtx*s2
622 t2(i)=vty*s2
623 t3(i)=vtz*s2
624 fact(i)=fric
625 ELSE
626 fact(i)=zero
627 t1(i)=one
628 t2(i)=zero
629 t3(i)=zero
630 ENDIF
631 ENDDO
632 DO i=1,jlt
633 dt(i) = t1(i)*dx(i) + t2(i)*dy(i) + t3(i)*dz(i)
634 dti(i) = t1(i)*dxi(i) + t2(i)*dyi(i) + t3(i)*dzi(i)
635 ENDDO
636 IF (scalk<0) THEN
637 isf=1
638 ELSE
639 isf=0
640 ENDIF
641 facf=abs(scalk)
642 IF (isf==1) THEN
643 DO i=1,jlt
644 IF (vn(i)>zero) THEN
645 facn(i)=stif(i)*facf
646 ELSEIF (vn(i)<zero) THEN
647 facn(i)=stif(i)/facf
648 ELSE
649 facn(i)=stif(i)
650 ENDIF
651 fact(i)=facn(i)*fact(i)
652 ENDDO
653 ELSE
654 DO i=1,jlt
655 facn(i)=stif(i)*facf
656 fact(i)=facn(i)*fact(i)
657 ENDDO
658 ENDIF
659C-------- part nml --------
660 DO i=1,jlt
661 fn = -facn(i)*dni(i)
662 fx=fn*n1(i)
663 fy=fn*n2(i)
664 fz=fn*n3(i)
665 IF (fact(i)/=zero) THEN
666 ft = -fact(i)*dti(i)
667 fx = fx + ft*t1(i)
668 fy = fy + ft*t2(i)
669 fz = fz + ft*t3(i)
670 ENDIF
671 a(1,ix1(i))=a(1,ix1(i))+fx*h1(i)
672 a(1,ix2(i))=a(1,ix2(i))+fx*h2(i)
673 a(1,ix3(i))=a(1,ix3(i))+fx*h3(i)
674 a(1,ix4(i))=a(1,ix4(i))+fx*h4(i)
675 a(2,ix1(i))=a(2,ix1(i))+fy*h1(i)
676 a(2,ix2(i))=a(2,ix2(i))+fy*h2(i)
677 a(2,ix3(i))=a(2,ix3(i))+fy*h3(i)
678 a(2,ix4(i))=a(2,ix4(i))+fy*h4(i)
679 a(3,ix1(i))=a(3,ix1(i))+fz*h1(i)
680 a(3,ix2(i))=a(3,ix2(i))+fz*h2(i)
681 a(3,ix3(i))=a(3,ix3(i))+fz*h3(i)
682 a(3,ix4(i))=a(3,ix4(i))+fz*h4(i)
683 ENDDO
684C-------- part nsl --------
685 DO i=1,jlt
686 fni = facn(i)*dn(i)
687 ni = index(i)
688 ffi(1,ni)=ffi(1,ni)+fni*n1(i)
689 ffi(2,ni)=ffi(2,ni)+fni*n2(i)
690 ffi(3,ni)=ffi(3,ni)+fni*n3(i)
691 IF (fact(i)/=zero) THEN
692 fti = fact(i)*dt(i)
693 ffi(1,ni)=ffi(1,ni)+fti*t1(i)
694 ffi(2,ni)=ffi(2,ni)+fti*t2(i)
695 ffi(3,ni)=ffi(3,ni)+fti*t3(i)
696 ENDIF
697 ENDDO
698C
699 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, J, K, 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
107C-----------------------------------------------
108C dev021
109 IF (iresp==1) THEN
110 prec = fiveem4
111 ELSE
112 prec = em10
113 ENDIF
114C--------------------------------------------------------
115C case of mixed packets
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 spherical (only concave for now)
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 cylindrical (only concave for now)
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 ---- Without friction first ---
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 with friction
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 708 of file i7keg3.F.

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