OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
i17for3.F File Reference
#include "implicit_f.inc"
#include "comlock.inc"
#include "mvsiz_p.inc"
#include "task_c.inc"
#include "com01_c.inc"
#include "com04_c.inc"
#include "parit_c.inc"
#include "mic_lockon.inc"
#include "mic_lockoff.inc"
#include "scr07_c.inc"
#include "scr11_c.inc"
#include "scr14_c.inc"
#include "scr16_c.inc"
#include "com06_c.inc"
#include "com08_c.inc"
#include "lockon.inc"
#include "lockoff.inc"

Go to the source code of this file.

Functions/Subroutines

subroutine i17for3 (x, v, candn, cande, i_stok, ixs, ixs16, eminxm, neles, nelem, itask, a, itied, nint, eminxs, stifn, fskyi, isky, nme, nse, frotm, frots, km, ks, fric, fsav, fcont, ms, niskyfi, lskyi17, noint, h3d_data)
subroutine i17lll_pena (llt, v, stifn, xx, fric, yy, zz, iii, fskyi, isky, a, x, itied, nint, xxs, yys, zzs, iiis, vit_min, le, les, ie_min, ies_min, itask, frotm, frots, km, ks, fsav, fcont, ms, niskyfi, noint, h3d_data)
subroutine i17mini_pena (llt, r_cs, s_cs, t_cs, ri_s, si_s, ti_s, ni_s, xxs, yys, zzs, xx, yy, zz, r_cm, s_cm, t_cm, nx, ny, nz, r_1s, r_2s, t_1s, t_2s, r_1m, r_2m, r_3m, r_4m, t_1m, t_2m, t_3m, t_4m, icont, area, area_tot, area_el)
subroutine i17abc_pena (i, f, r, t, b1, b2, b3, c1, c2, c3)
subroutine i17borne_pena (llt, r_s, a, b, c, icont, rs)
subroutine i17racine_pena (llt, a, b, c, r1, r2)
subroutine i17surf (llt, r1, r2, r3, r4, t1, t2, t3, t4, area, xx, yy, zz, sm)
subroutine i17vit_pena (llt, nint, v, a, iii, iiis, ni_m, ni_s, nx, ny, nz, vit, icont, rm, tm, rs, ts, sm, les)
subroutine i17lll4_pena (llt, itied, nint, v, a, fric, iii, iiis, ni_m, ni_s, s_cm, s_cs, nx, ny, nz, vit, le, les, icont, rm, tm, rs, ts, xx, yy, zz, stifn, fskyi, isky, frotm, frots, area, area_tot, km, ks, fsav, fcont, ms, area_el, niskyfi, noint, h3d_data)
subroutine i17lll_ass (nskyi17, iskyi17, fskyi17, frots, nmes, lskyi17)

Function/Subroutine Documentation

◆ i17abc_pena()

subroutine i17abc_pena ( integer i,
f,
r,
t,
b1,
b2,
b3,
c1,
c2,
c3 )

Definition at line 934 of file i17for3.F.

936C-----------------------------------------------
937C i=1,4
938C ri=+-1 ti=+-1
939C Ni = 1/4 (1+ro)(1+to)(ro+to-1)
940C Ni = 1/4 (r^2 + ti r^2 t + t^2 + ri ti r t + ri r t^2 - 1)
941C
942C i=6;8
943C ri=0 ti=+-1
944C Ni = 1/2 ( 1 + ti t - r^2 - ti t r^2)
945C
946C
947C i=5;7
948C ri=+-1 ti=0
949C Ni = 1/2 (1 + ri r - t^2 - ri r t^2)
950C-----------------------------------------------
951C f = Somme( fi Ni )
952C
953C f = A1 + A2 r + A3 t + A4 rt + A5 r^2 + A6 t^2 + A7 tr^2 + A8 rt^2
954C f = (A1 + A3 t + A6 t^2) + (A2 + A4 t + A8 t^2)r + (A5 + A7 t)r^2
955C f = B1 + B2 r + B3 r^2
956C f = (A1 + A2 r + A5 r^2) + (A3 + A4 r + A7 r^2)t + (A6 + A8 r)t^2
957C f = C1 + C2 t + C3 t^2
958C
959C A1 = (-f1 - f2 - f3 - f4 + 2 f5 + 2 f6 + 2 f7 + 2 f8)/4
960C A2 = ( - f5 + f7 )/2
961C A3 = ( + f6 - f8)/2
962C A4 = (+f1 - f2 + f3 - f4 )/4
963C A5 = (+f1 + f2 + f3 + f4 - 2 f6 - 2 f8)/4
964C A6 = (+f1 + f2 + f3 + f4 - 2 f5 - 2 f7 )/4
965C A7 = (-f1 + f2 + f3 - f4 - 2 f6 + 2 f8)/4
966C A8 = (-f1 - f2 + f3 + f4 + 2 f5 - 2 f7 )/4
967C-----------------------------------------------
968C I m p l i c i t T y p e s
969C-----------------------------------------------
970#include "implicit_f.inc"
971C-----------------------------------------------
972C G l o b a l P a r a m e t e r s
973C-----------------------------------------------
974#include "mvsiz_p.inc"
975C-----------------------------------------------
976C D u m m y A r g u m e n t s
977C-----------------------------------------------
978 INTEGER I
979 my_real
980 + f(mvsiz,*),r,t,
981 + b1,b2,b3,
982 + c1,c2,c3
983C-----------------------------------------------
984C L o c a l V a r i a b l e s
985C-----------------------------------------------
986 my_real
987 + a1,a2,a3,a4,a5,a6,a7,a8,r2,t2,ff5,ff6,ff7,ff8
988C
989 ff5 = f(i,5) + f(i,5)
990 ff6 = f(i,6) + f(i,6)
991 ff7 = f(i,7) + f(i,7)
992 ff8 = f(i,8) + f(i,8)
993c
994c A1 = (-F(I,1)-F(I,2)-F(I,3)-F(I,4)+FF5+FF6+FF7+FF8)*0.25 ...
995c B1(I) = A1 + ( A3 + A6 * T(I) ) *T(I) ...
996c
997 a1 = (-f(i,1)-f(i,2)-f(i,3)-f(i,4)+ff5+ff6+ff7+ff8)
998 a2 = ( -ff5 +ff7 )
999 a3 = ( +ff6 -ff8)
1000 a4 = (+f(i,1)-f(i,2)+f(i,3)-f(i,4) )
1001 a5 = (+f(i,1)+f(i,2)+f(i,3)+f(i,4) -ff6 -ff8)
1002 a6 = (+f(i,1)+f(i,2)+f(i,3)+f(i,4)-ff5 -ff7 )
1003 a7 = (-f(i,1)+f(i,2)+f(i,3)-f(i,4) -ff6 +ff8)
1004 a8 = (-f(i,1)-f(i,2)+f(i,3)+f(i,4)+ff5 -ff7 )
1005c
1006 b1 = ( a1 + ( a3 + a6 * t ) *t )*fourth
1007 b2 = ( a2 + ( a4 + a8 * t ) *t )*fourth
1008 b3 = ( a5 + a7 * t )*fourth
1009c
1010 c1 = ( a1 + ( a2 + a5 * r ) *r )*fourth
1011 c2 = ( a3 + ( a4 + a7 * r ) *r )*fourth
1012 c3 = ( a6 + a8 * r )*fourth
1013c
1014 RETURN
my_real
#define my_real
Definition cppsort.cpp:32

◆ i17borne_pena()

subroutine i17borne_pena ( integer llt,
r_s,
a,
b,
c,
integer, dimension(mvsiz) icont,
rs )

Definition at line 1021 of file i17for3.F.

1022C-----------------------------------------------
1023C I m p l i c i t T y p e s
1024C-----------------------------------------------
1025#include "implicit_f.inc"
1026C-----------------------------------------------
1027C G l o b a l P a r a m e t e r s
1028C-----------------------------------------------
1029#include "mvsiz_p.inc"
1030C-----------------------------------------------
1031C D u m m y A r g u m e n t s
1032C-----------------------------------------------
1033 INTEGER LLT,ICONT(MVSIZ)
1034 my_real
1035 + r_s(mvsiz),c(mvsiz),b(mvsiz),a(mvsiz),rs(mvsiz)
1036C-----------------------------------------------
1037C L o c a l V a r i a b l e s
1038C-----------------------------------------------
1039 INTEGER I
1040 my_real
1041 + cc(mvsiz),r1(mvsiz),r2(mvsiz),rs1,rs2,unpeps
1042c-----------------------------------------------------------------------
1043c 1: borne r(ou s au 2eme appel) a +-1 sur le main
1044c-----------------------------------------------------------------------
1045c-----------------------------------------------------------------------
1046c 1.1: borne <r1,t> et <r,t1> a +-1 sur le main
1047c-----------------------------------------------------------------------
1048 unpeps = one + em4
1049 DO i=1,llt
1050 rs(i) = c(i) + ( b(i) + a(i)*r_s(i) )*r_s(i)
1051 IF(rs(i)>one)THEN
1052 cc(i) = c(i) - one
1053 ELSEIF(rs(i)<-one)THEN
1054 cc(i) = c(i) + one
1055 ELSE
1056 cc(i) = one
1057 ENDIF
1058 ENDDO
1059c
1060C f = C + B r + A r^2
1061c
1062 CALL i17racine(llt,a,b,cc,r1,r2)
1063c
1064 DO i=1,llt
1065 IF(rs(i)>unpeps.OR.rs(i)<-unpeps)THEN
1066 rs1 = c(i) + ( b(i) + a(i)*r1(i) )*r1(i)
1067 rs2 = c(i) + ( b(i) + a(i)*r2(i) )*r2(i)
1068 IF(rs1>=-unpeps.AND.rs1<=unpeps.AND.
1069 + rs2>=-unpeps.AND.rs2<=unpeps)THEN
1070 IF(abs(rs(i)-r1(i))<abs(rs(i)-r2(i)))THEN
1071 r_s(i) = r1(i)
1072 rs(i) = rs1
1073 ELSE
1074 r_s(i) = r2(i)
1075 rs(i) = rs2
1076 ENDIF
1077 ELSEIF(rs1>=-unpeps.AND.rs1<=unpeps)THEN
1078 r_s(i) = r1(i)
1079 rs(i) = rs1
1080 ELSEIF(rs2>=-unpeps.AND.rs2<=unpeps)THEN
1081 r_s(i) = r2(i)
1082 rs(i) = rs2
1083 ELSE
1084C pas de recouvrement main/second
1085 icont(i)=0
1086 ENDIF
1087 ENDIF
1088 ENDDO
1089c
1090 RETURN
i17racine
subroutine i17racine(llt, a, b, c, r1, r2)
Definition i17lagm.F:1316

◆ i17for3()

subroutine i17for3 ( x,
v,
integer, dimension(*) candn,
integer, dimension(*) cande,
integer i_stok,
integer, dimension(nixs,*) ixs,
integer, dimension(8,*) ixs16,
eminxm,
integer, dimension(*) neles,
integer, dimension(*) nelem,
integer itask,
a,
integer itied,
integer nint,
eminxs,
stifn,
fskyi,
integer, dimension(*) isky,
integer nme,
integer nse,
frotm,
frots,
km,
ks,
fric,
fsav,
fcont,
ms,
integer niskyfi,
integer lskyi17,
integer noint,
type(h3d_database) h3d_data )

Definition at line 41 of file i17for3.F.

49C-----------------------------------------------
50C M o d u l e s
51C-----------------------------------------------
52 USE tri7box
53 USE icontact_mod
54 USE message_mod
55 USE h3d_mod
56 USE anim_mod
57C-----------------------------------------------
58C I m p l i c i t T y p e s
59C-----------------------------------------------
60#include "implicit_f.inc"
61#include "comlock.inc"
62C-----------------------------------------------
63C G l o b a l P a r a m e t e r s
64C-----------------------------------------------
65#include "mvsiz_p.inc"
66C-----------------------------------------------
67C C o m m o n B l o c k s
68C-----------------------------------------------
69#include "task_c.inc"
70#include "com01_c.inc"
71#include "com04_c.inc"
72#include "parit_c.inc"
73 COMMON /i17globi/ie_min,ies_min
74 COMMON /i17globr/vit_min
75C-----------------------------------------------
76C D u m m y A r g u m e n t s
77C-----------------------------------------------
78 INTEGER I_STOK,ITASK,ITIED,NINT,NME,NSE,NISKYFI,LSKYI17,NOINT,
79 . CANDN(*),CANDE(*), ISKY(*),
80 . IXS(NIXS,*),IXS16(8,*),NELES(*) ,NELEM(*)
81C REAL
82 my_real
83 . x(3,*),v(3,*),eminxm(6,*),eminxs(6,*),a(3,*),stifn(*),
84 . fskyi ,frotm(7,*),frots(7,*),km(*),ks(*),fric(*), fsav(*),
85 . fcont(3,*), ms(*)
86 TYPE(H3D_DATABASE) :: H3D_DATA
87C-----------------------------------------------
88C L o c a l V a r i a b l e s
89C-----------------------------------------------
90 INTEGER I,J,K,IK,IE,IS,IC,NK,III(MVSIZ,17),LLT,NFT,LE,FIRST,LAST,
91 . I16,LES,IES,IIIS(MVSIZ,16),LEL(MVSIZ),LESL(MVSIZ),
92 . IE_MIN,IES_MIN, IERR1, IERR2
93 my_real
94 . xx(mvsiz,17),yy(mvsiz,17),zz(mvsiz,17),
95 . xxs(mvsiz,16),yys(mvsiz,16),zzs(mvsiz,16),
96 . vit_min
97C-----------------------------------------------------------------------
98C boucle sur les candidats au contact
99C-----------------------------------------------------------------------
100 vit_min = zero
101 ie_min = 99999999
102 ies_min = 99999999
103 nskyi17 = 0
104 nrskyi17 = 0
105 IF(iparit/=0.AND.itask==0)THEN
106 ierr1=0
107 ierr2=0
108 ALLOCATE(iskyi17(lskyi17),stat=ierr2)
109 ierr1 = ierr1+ierr2
110 ALLOCATE(fskyi17(lskyi17,4),stat=ierr2)
111 ierr1 = ierr1+ierr2
112 IF(nspmd>1)THEN
113 ALLOCATE(irskyi17(lskyi17),stat=ierr2)
114 ierr1 = ierr1+ierr2
115 ALLOCATE(frskyi17(4,lskyi17),stat=ierr2)
116 ierr1 = ierr1+ierr2
117 END IF
118 IF(ierr1/=0)THEN
119 CALL ancmsg(msgid=88,anmode=aninfo,
120 . i1=noint)
121 CALL arret(2)
122 END IF
123 END IF
124C-----------------------------------------------------------------------
125 CALL my_barrier
126C-----------------------------------------------------------------------
127 first = 1 + i_stok * itask / nthread
128 last = i_stok*(itask+1) / nthread
129 llt = 0
130 nft=llt+1
131 DO ic=first,last
132 le = cande(ic)
133 les = candn(ic)
134 ie = nelem(le)
135 IF(les<=nse)THEN ! candidat interne en SPMD
136 ies = neles(les)
137C-----------------------------------------------------------------------
138C test si dans la boite
139C-----------------------------------------------------------------------
140 IF(le >0.AND.les>0.AND.
141 . eminxs(4,les)>eminxm(1,le).AND.
142 . eminxs(5,les)>eminxm(2,le).AND.
143 . eminxs(6,les)>eminxm(3,le).AND.
144 . eminxs(1,les)<eminxm(4,le).AND.
145 . eminxs(2,les)<eminxm(5,le).AND.
146 . eminxs(3,les)<eminxm(6,le))THEN
147 llt = llt+1
148 lel(llt)=le
149 lesl(llt)=les
150 DO k=1,4
151C III 1,2,3,4,5,6,7,8, 9,10,11,12,13,14,15,16
152 iii(llt,k) =ixs(k+1,ie)
153 iii(llt,k+4) =ixs(k+5,ie)
154 iii(llt,k+8) =ixs16(k,ie-numels8-numels10-numels20)
155 iii(llt,k+12)=ixs16(k+4,ie-numels8-numels10-numels20)
156C IIIS 1,2,3,4,9,10,11,12, 5,6,7,8,13,14,15,16
157 iiis(llt,k) =ixs(k+1,ies)
158 iiis(llt,k+8) =ixs(k+5,ies)
159 iiis(llt,k+4)=ixs16(k,ies-numels8-numels10-numels20)
160 iiis(llt,k+12)=ixs16(k+4,ies-numels8-numels10-numels20)
161 ENDDO
162 DO k=1,16
163 i = iii(llt,k)
164 xx(llt,k)=x(1,i)
165 yy(llt,k)=x(2,i)
166 zz(llt,k)=x(3,i)
167 i = iiis(llt,k)
168 xxs(llt,k)=x(1,i)
169 yys(llt,k)=x(2,i)
170 zzs(llt,k)=x(3,i)
171 ENDDO
172 END IF
173 ELSE ! complement frontiere SPMD
174 ies = les-nse
175C-----------------------------------------------------------------------
176C test si dans la boite
177C-----------------------------------------------------------------------
178 IF(le >0.AND.les>0.AND.
179 . eminxfi(nint)%P(4,ies)>eminxm(1,le).AND.
180 . eminxfi(nint)%P(5,ies)>eminxm(2,le).AND.
181 . eminxfi(nint)%P(6,ies)>eminxm(3,le).AND.
182 . eminxfi(nint)%P(1,ies)<eminxm(4,le).AND.
183 . eminxfi(nint)%P(2,ies)<eminxm(5,le).AND.
184 . eminxfi(nint)%P(3,ies)<eminxm(6,le))THEN
185 llt = llt+1
186 lel(llt)=le
187 lesl(llt)=-ies ! reperage en negatif des candidats remote SPMD
188#include "mic_lockon.inc"
189 nsvfi(nint)%P(ies) = -abs(nsvfi(nint)%P(ies))
190#include "mic_lockoff.inc"
191 DO k=1,4
192C III 1,2,3,4,5,6,7,8, 9,10,11,12,13,14,15,16
193 iii(llt,k) =ixs(k+1,ie)
194 iii(llt,k+4) =ixs(k+5,ie)
195 iii(llt,k+8) =ixs16(k,ie-numels8-numels10-numels20)
196 iii(llt,k+12)=ixs16(k+4,ie-numels8-numels10-numels20)
197C IIIS 1,2,3,4,9,10,11,12, 5,6,7,8,13,14,15,16
198 iiis(llt,k) =k
199 iiis(llt,k+8) =k+4
200 iiis(llt,k+4) =k+8
201 iiis(llt,k+12)=k+12
202 ENDDO
203 DO k=1,16
204 i = iii(llt,k)
205 xx(llt,k)=x(1,i)
206 yy(llt,k)=x(2,i)
207 zz(llt,k)=x(3,i)
208 i = iiis(llt,k) ! I indique le no de noeud de 1 a 16
209 xxs(llt,k)=xfi17(nint)%P(1,i,ies)
210 yys(llt,k)=xfi17(nint)%P(2,i,ies)
211 zzs(llt,k)=xfi17(nint)%P(3,i,ies)
212 ENDDO
213 END IF
214 END IF
215C-----------------------------------------------------------------------
216C calcul de F par paquet de mvsiz
217C-----------------------------------------------------------------------
218 IF(llt==mvsiz-1)THEN
219 CALL i17lll_pena(
220 1 llt ,v ,stifn ,xx ,fric ,
221 2 yy ,zz ,iii ,fskyi ,isky ,
222 3 a ,x ,itied ,nint ,
223 4 xxs ,yys ,zzs ,iiis ,vit_min ,
224 5 lel ,lesl ,ie_min ,ies_min ,itask ,
225 6 frotm ,frots ,km ,ks ,fsav ,
226 7 fcont ,ms ,niskyfi ,noint ,h3d_data )
227 nft=llt+1
228 llt = 0
229 ENDIF
230 ENDDO
231C-----------------------------------------------------------------------
232C calcul de F pour dernier paquet
233C-----------------------------------------------------------------------
234 IF(llt/=0) CALL i17lll_pena(
235 1 llt ,v ,stifn ,xx ,fric ,
236 2 yy ,zz ,iii ,fskyi ,isky ,
237 3 a ,x ,itied ,nint ,
238 4 xxs ,yys ,zzs ,iiis ,vit_min ,
239 5 lel ,lesl ,ie_min ,ies_min ,itask ,
240 6 frotm ,frots ,km ,ks ,fsav ,
241 7 fcont ,ms ,niskyfi ,noint ,h3d_data )
242C-----------------------------------------------------------------------
243C parallel Arithmetic for assembly (frots)
244C-----------------------------------------------------------------------
245 IF(iparit/=0) THEN
246 CALL my_barrier
247 IF(itask==0)THEN
248 IF(nspmd>1)THEN
249 CALL spmd_i17frots_pon(
250 1 nskyi17 ,iskyi17,fskyi17,nrskyi17,irskyi17,
251 2 frskyi17,nint ,lskyi17,noint)
252 DEALLOCATE(irskyi17)
253 DEALLOCATE(frskyi17)
254 END IF
255C
256 IF(nskyi17>0)THEN
257 CALL i17lll_ass(nskyi17 ,iskyi17 ,fskyi17, frots, nse,
258 2 lskyi17 )
259 END IF
260 DEALLOCATE(iskyi17)
261 DEALLOCATE(fskyi17)
262 END IF
263 END IF
264C
265 RETURN
i17lll_ass
subroutine i17lll_ass(nskyi17, iskyi17, fskyi17, frots, nmes, lskyi17)
Definition i17for3.F:2054
i17lll_pena
subroutine i17lll_pena(llt, v, stifn, xx, fric, yy, zz, iii, fskyi, isky, a, x, itied, nint, xxs, yys, zzs, iiis, vit_min, le, les, ie_min, ies_min, itask, frotm, frots, km, ks, fsav, fcont, ms, niskyfi, noint, h3d_data)
Definition i17for3.F:288
anim_mod
Definition anim_mod.F:96
h3d_mod
Definition h3d_mod.F:308
icontact_mod
Definition icontact_mod.F:34
icontact_mod::nskyi17
integer nskyi17
Definition icontact_mod.F:42
icontact_mod::irskyi17
integer, dimension(:), allocatable irskyi17
Definition icontact_mod.F:43
icontact_mod::nrskyi17
integer nrskyi17
Definition icontact_mod.F:42
icontact_mod::iskyi17
integer, dimension(:), allocatable iskyi17
Definition icontact_mod.F:43
message_mod
Definition message_mod.F:1249
tri7box::eminxfi
type(real_pointer2), dimension(:), allocatable eminxfi
Definition tri7box.F:467
tri7box::xfi17
type(real_pointer3), dimension(:), allocatable xfi17
Definition tri7box.F:470
tri7box::nsvfi
type(int_pointer), dimension(:), allocatable nsvfi
Definition tri7box.F:431
spmd_i17frots_pon
subroutine spmd_i17frots_pon(nskyi17, iskyi17, fskyi17, nrskyi17, irskyi17, frskyi17, nin, lskyi17, noint)
Definition spmd_i17frots_pon.F:41
ancmsg
subroutine ancmsg(msgid, msgtype, anmode, i1, i2, i3, i4, i5, i6, i7, i8, i9, i10, i11, i12, i13, i14, i15, i16, i17, i18, i19, i20, r1, r2, r3, r4, r5, r6, r7, r8, r9, c1, c2, c3, c4, c5, c6, c7, c8, c9, prmode)
Definition message.F:889
arret
subroutine arret(nn)
Definition arret.F:87
my_barrier
subroutine my_barrier
Definition machine.F:31

◆ i17lll4_pena()

subroutine i17lll4_pena ( integer llt,
integer itied,
integer nint,
v,
a,
fric,
integer, dimension(mvsiz,17) iii,
integer, dimension(mvsiz,16) iiis,
ni_m,
ni_s,
s_cm,
s_cs,
nx,
ny,
nz,
vit,
integer, dimension(*) le,
integer, dimension(*) les,
integer, dimension(mvsiz) icont,
rm,
tm,
rs,
ts,
xx,
yy,
zz,
stifn,
fskyi,
integer, dimension(*) isky,
frotm,
frots,
area,
area_tot,
km,
ks,
fsav,
fcont,
ms,
area_el,
integer niskyfi,
integer noint,
type(h3d_database) h3d_data )

Definition at line 1498 of file i17for3.F.

1507C-----------------------------------------------
1508C M o d u l e s
1509C-----------------------------------------------
1510 USE tri7box
1511 USE icontact_mod
1512 USE message_mod
1513 USE h3d_mod
1514 USE anim_mod
1515C-----------------------------------------------
1516C I m p l i c i t T y p e s
1517C-----------------------------------------------
1518#include "implicit_f.inc"
1519#include "comlock.inc"
1520C-----------------------------------------------
1521C G l o b a l P a r a m e t e r s
1522C-----------------------------------------------
1523#include "mvsiz_p.inc"
1524C-----------------------------------------------
1525C C o m m o n B l o c k s
1526C-----------------------------------------------
1527#include "scr07_c.inc"
1528#include "scr11_c.inc"
1529#include "scr14_c.inc"
1530#include "scr16_c.inc"
1531#include "com06_c.inc"
1532#include "com08_c.inc"
1533#include "parit_c.inc"
1534C-----------------------------------------------
1535C D u m m y A r g u m e n t s
1536C-----------------------------------------------
1537 INTEGER LLT,ITIED,NINT,NISKYFI,NOINT
1538 INTEGER III(MVSIZ,17),IIIS(MVSIZ,16),ICONT(MVSIZ), ISKY(*),
1539 . LE(*) ,LES(*)
1540C REAL
1541 my_real
1542 . v(3,*),a(3,*),vit(*), ms(*)
1543 my_real
1544 . rm(mvsiz) ,rs(mvsiz) ,tm(mvsiz) ,ts(mvsiz) ,
1545 . xx(mvsiz,17),yy(mvsiz,17),zz(mvsiz,17),s_cm(mvsiz),
1546 . s_cs(mvsiz), fsav(*),fcont(3,*),
1547 . ni_m(mvsiz,*) ,ni_s(mvsiz,*),nx(mvsiz) ,ny(mvsiz) ,nz(mvsiz),
1548 . stifn(*),fskyi(lskyi,nfskyi),frotm(7,*),frots(7,*),
1549 . area(mvsiz),area_tot(mvsiz),area_el(mvsiz),km(2,*),ks(2,*),
1550 . fric(*)
1551 TYPE(H3D_DATABASE) :: H3D_DATA
1552C-----------------------------------------------
1553C L o c a l V a r i a b l e s
1554C-----------------------------------------------
1555 INTEGER I,J,IK,NK,I1,I2,I3,I4,IAD,NN,LEM,J1,J2,NISKYL,IES,IIS,
1556 . NISKYFIL
1557 my_real
1558 . vx,vy,vz,vn,aa,surf,ep2,xe,ye,ze,xa,ya,za,xb,yb,zb,xc,yc,zc,
1559 . f,stif,pene,sigtmx,sigtmy,sigtmz,sigtsx,sigtsy,sigtsz,bb,
1560 . fx,fy,fz,ffx,ffy,ffz,ff,ff2,muf2,htsqrtpi,mas4,ep,vis,
1561 . rhom,rhos,stifv,dt1inv,ffvx,ffvy,ffvz,fftx,ffty,fftz,muf,
1562 . fsav1,fsav2,fsav3,fsav4,fsav5,fsav6,econvt,econtt,beta,
1563 . ks1,ks2
1564c-----------------------------------------------------------------------
1565cfrottement
1566c-----------------------------------------------------------------------
1567c0- formulation classique
1568c
1569c Ki = 8/3 E* (pi/A)^1/2 Ai/pi
1570c
1571c Fti = min[(Fti + Ki * Vti * dt) , mu*Fni]
1572c
1573c => limitation en cas de roulement glissement ou roulement adh rence:
1574c si Fti est m moris e sur l' l ment main ou second, Fti est perdue
1575c quand on passe d'un l ment sur son voisin
1576c
1577c1- on consid re une facette main contact e par plusieurs seconds
1578c
1579c
1580c Fti = min[(Sigt*Ai + Ki * Vti * dt) , mu*Fni]
1581c
1582c Sigt = Somme(Fti)/Somme(Ai)
1583c
1584c ou (estimation de A en i): Sigt = Somme(Fti/A)
1585c
1586c * Sigt vecteur contrainte de frottement
1587c
1588c
1589c2- sym trie main second
1590c
1591c Fti = min[((SigtM+SigtS)*Ai/2 + Ki * Vti * dt) , mu*Fni]
1592c
1593c SigtM = SommeM(Fti)/SommeM(Ai)
1594c SigtS = SommeS(Fti)/SommeS(Ai)
1595C-----------------------------------------------
1596 htsqrtpi = eight / (three*sqrt(pi))
1597
1598 IF(dt1>zero)THEN
1599 dt1inv = one/dt1
1600 ELSE
1601 dt1inv =zero
1602 ENDIF
1603
1604C-----------------------------------------------
1605 fsav1 = zero
1606 fsav2 = zero
1607 fsav3 = zero
1608 fsav4 = zero
1609 fsav5 = zero
1610 fsav6 = zero
1611 econvt = zero
1612 econtt = zero
1613 DO i=1,llt
1614 lem = le(i)
1615 IF(icont(i)/=0 )THEN
1616c---------------------------------------------
1617c contact de Hertz
1618c---------------------------------------------
1619c
1620c ellipse de contact 2a x 2b => surface A = a b pi
1621c r: rayon quivalent de l'ellipse r^2 = a b
1622c rayon de courbure 1/R = 1/R1 + 1/R2 + 1/R3 + 1/R4
1623c k1 = (1-v1^2)/E1 k2 = (1-v2^2)/E2 k = k1+k2 = 1/E*
1624c kk1 = (1-v1^2)/piE1 kk2 = (1-v2^2)/piE2 kk = kk1+kk2
1625c E* = 1/(k1+k2)
1626c F: force
1627c pm: pression max au centre de l'ellipse)
1628c d: penetration max
1629c------------------------------------------
1630cd'apres http://fr.wikibooks.org/wiki/Tribologie_-_Contacts_localis%C3%A9s
1631c
1632c pm = 3/2 F/A (1)
1633c
1634c p(x,y) = pm [1 - x^2/a^2 - y^2/b^2]^1/3 (2)
1635c
1636c sur la fronti re p=0 => x^2/a^2 + y^2/b^2 = 1
1637c
1638cd'apres http://barreau.matthieu.free.fr/cours/Liaison-pivot/pages/roulements-1.html
1639c
1640c
1641c F = pm^3 pi^3 R^2 k^2 / 6 (4)
1642c
1643c (1)+(4)
1644c F = 4 A / (3 pi R k) [A/pi]^1/2
1645c
1646c-----------------------
1647c F = 4 A r / (3 pi R k) (5)
1648c F = 4 r^3 / (3 R k) (6)
1649c-----------------------
1650c--------------------------------------------------------------
1651cp n tration, force, rigidit
1652c--------------------------------------------------------------
1653c
1654c d = R - (R^2 - r^2)^1/2
1655c d = R[1 - (1 - r^2/R^2)^1/2]
1656c-----------------------
1657c d ~= 1/2 r^2/R (8)
1658c-----------------------
1659c F = 8 r/ (3 k) d (9)
1660c-----------------------
1661c F = K d (10)
1662c-----------------------
1663c K = 8 r/ (3 k) (11)
1664c-------------------------------------------
1665c---------------------------------------------
1666c contact multiple (i=1,n)
1667c si on connait la courbure
1668c---------------------------------------------
1669c pour chaque contact on evalue ri^2 , Ri et di
1670c
1671c Rmax = ri^2 / 2 di
1672c Ki = 4/3 E* (2/ (Ri di))^1/2 ri^2
1673c Kmin = 8/3 E* ri
1674c Fi = Ki di
1675c
1676c Rmax = Ai / 2pi di
1677c Ki = 4/3 E* (2/(Ri di))^1/2 Ai/pi
1678c Kmin = 8/3 E* (Ai/pi)^1/2
1679c Fi = max(Ki,Kmin) di
1680c---------------------------------------------
1681c contact multiple (i=1,n)
1682c si on connait la surface totale
1683c---------------------------------------------
1684c A: surface totale du contact (estimee)
1685c
1686c Ri = A / 2pi di
1687c Ki = 8/3 E* (pi/A)^1/2 Ai/pi
1688c Fi = Ki di
1689C-----------------------------------------------
1690C epaisseur main
1691C-----------------------------------------------
1692 xe = zero
1693 ye = zero
1694 ze = zero
1695 DO ik=1,4
1696 xe = xe + xx(i,ik) *ni_m(i,ik)
1697 + + xx(i,ik+8) *ni_m(i,ik+4)
1698 + - xx(i,ik+4) *ni_m(i,ik)
1699 + - xx(i,ik+12)*ni_m(i,ik+4)
1700 ye = ye + yy(i,ik) *ni_m(i,ik)
1701 + + yy(i,ik+8) *ni_m(i,ik+4)
1702 + - yy(i,ik+4) *ni_m(i,ik)
1703 + - yy(i,ik+12)*ni_m(i,ik+4)
1704 ze = ze + zz(i,ik) *ni_m(i,ik)
1705 + + zz(i,ik+8) *ni_m(i,ik+4)
1706 + - zz(i,ik+4) *ni_m(i,ik)
1707 + - zz(i,ik+12)*ni_m(i,ik+4)
1708 ENDDO
1709 ep = abs(xe*nx(i) + ye*ny(i) + ze*nz(i))
1710 pene = half*(one-abs(s_cm(i))) * ep
1711 if(pene<zero)then
1712 stop 111
1713 endif
1714C-----------------------------------------------
1715C vitesse relative second/main
1716C-----------------------------------------------
1717 vx = zero
1718 vy = zero
1719 vz = zero
1720 IF(les(i)>0)THEN
1721 DO ik=1,8
1722 vx = vx + (v(1,iiis(i,ik)))*ni_s(i,ik)
1723 . - (v(1,iii(i,ik))) *ni_m(i,ik)
1724 vy = vy + (v(2,iiis(i,ik)))*ni_s(i,ik)
1725 . - (v(2,iii(i,ik))) *ni_m(i,ik)
1726 vz = vz + (v(3,iiis(i,ik)))*ni_s(i,ik)
1727 . - (v(3,iii(i,ik))) *ni_m(i,ik)
1728 ENDDO
1729 ks1 = ks(1,les(i))
1730 ks2 = ks(2,les(i))
1731C
1732 sigtmx = frotm(5,lem)
1733 sigtmy = frotm(6,lem)
1734 sigtmz = frotm(7,lem)
1735 sigtsx = frots(5,les(i))
1736 sigtsy = frots(6,les(i))
1737 sigtsz = frots(7,les(i))
1738 ELSE ! partie complementaire SPMD
1739 ies = abs(les(i)) ! no local
1740 DO ik=1,8
1741 iis = iiis(i,ik)
1742 vx = vx + vfi17(nint)%P(1,iis,ies)*ni_s(i,ik)
1743 . - v(1,iii(i,ik)) *ni_m(i,ik)
1744 vy = vy + vfi17(nint)%P(2,iis,ies)*ni_s(i,ik)
1745 . - v(2,iii(i,ik)) *ni_m(i,ik)
1746 vz = vz + vfi17(nint)%P(3,iis,ies)*ni_s(i,ik)
1747 . - v(3,iii(i,ik)) *ni_m(i,ik)
1748 ENDDO
1749 ks1 = ksfi(nint)%P(1,ies)
1750 ks2 = ksfi(nint)%P(2,ies)
1751C
1752 sigtmx = frotm(5,lem)
1753 sigtmy = frotm(6,lem)
1754 sigtmz = frotm(7,lem)
1755 sigtsx = frotsfi(nint)%P(5,ies)
1756 sigtsy = frotsfi(nint)%P(6,ies)
1757 sigtsz = frotsfi(nint)%P(7,ies)
1758 END IF
1759 vn = nx(i)*vx + ny(i)*vy + nz(i)*vz
1760C-----------------------------------------------
1761C rigidite Hertz
1762C-----------------------------------------------
1763c Ki = 8 E* Ai / 3(A*pi)^1/2
1764c HTSQRTPI = EIGHT / (THREE*SQRT(PI))
1765 stif = htsqrtpi*area(i) /
1766 . ((km(1,lem)+ks1)*sqrt(area_tot(i)))
1767 f = stif*pene
1768C-----------------------------------------------
1769C Viscosite critique
1770C-----------------------------------------------
1771 rhom = km(2,lem)
1772 rhos = ks2
1773 mas4 = (rhom+rhos) * area(i) * min(ep,sqrt(area(i)))
1774 beta=0.5
1775 vis = beta*sqrt(stif*mas4)
1776 f = f + vis * vn
1777 stifv = vis*dt1inv
1778C-----------------------------------------------
1779C force normale
1780C-----------------------------------------------
1781 fx = nx(i) * f
1782 fy = ny(i) * f
1783 fz = nz(i) * f
1784C-----------------------------------------------
1785C vitesse tangente
1786C-----------------------------------------------
1787c VX = VX - VN*NX(I)
1788c VY = VY - VN*NY(I)
1789c VZ = VZ - VN*NZ(I)
1790C-----------------------------------------------
1791C frottement
1792C calcul avec vitesse totale
1793C et projection de la force
1794C-----------------------------------------------
1795 aa = half*area(i)
1796 bb = stif*dt1
1797 ffx = aa*(sigtmx+sigtsx)+bb*vx
1798 ffy = aa*(sigtmy+sigtsy)+bb*vy
1799 ffz = aa*(sigtmz+sigtsz)+bb*vz
1800 ffvx = vis*vx
1801 ffvy = vis*vy
1802 ffvz = vis*vz
1803C-----------------------------------------------
1804C force tangente
1805C-----------------------------------------------
1806 ff = nx(i)*ffx + ny(i)*ffy + nz(i)*ffz
1807 ffx = ffx - ff*nx(i)
1808 ffy = ffy - ff*ny(i)
1809 ffz = ffz - ff*nz(i)
1810 ff = nx(i)*ffvx + ny(i)*ffvy + nz(i)*ffvz
1811 ffvx = ffvx - ff*nx(i)
1812 ffvy = ffvy - ff*ny(i)
1813 ffvz = ffvz - ff*nz(i)
1814 fftx = ffx + ffvx
1815 ffty = ffy + ffvy
1816 fftz = ffz + ffvz
1817C-----------------------------------------------
1818C frottement de Coulomb (penalite viscoelastique)
1819C-----------------------------------------------
1820 ff2 = fftx*fftx
1821 . + ffty*ffty
1822 . + fftz*fftz
1823 muf = fric(1)*max(zero,f)
1824 muf2 = muf*muf
1825 IF(ff2 > muf2)THEN
1826c reduction de la contrainte visqueuse
1827c puis de la contrainte elastique (si necessaire)
1828 aa = muf/sqrt(ff2)
1829 fftx = fftx*aa
1830 ffty = ffty*aa
1831 fftz = fftz*aa
1832 ff2 = ffx*ffx + ffy*ffy + ffz*ffz
1833 IF(ff2 > muf2)THEN
1834 aa = muf / sqrt(ff2)
1835 ffx = ffx*aa
1836 ffy = ffy*aa
1837 ffz = ffz*aa
1838 ENDIF
1839 ENDIF
1840
1841#include "lockon.inc"
1842 frotm(1,lem) = frotm(1,lem) + ffx
1843 frotm(2,lem) = frotm(2,lem) + ffy
1844 frotm(3,lem) = frotm(3,lem) + ffz
1845 frotm(4,lem) = frotm(4,lem) + area(i)
1846 IF(les(i)>0)THEN
1847 IF(iparit == 0)THEN
1848 frots(1,les(i)) = frots(1,les(i)) + ffx
1849 frots(2,les(i)) = frots(2,les(i)) + ffy
1850 frots(3,les(i)) = frots(3,les(i)) + ffz
1851 frots(4,les(i)) = frots(4,les(i)) + area(i)
1852 ELSE
1853 nskyi17 = nskyi17+1
1854 iskyi17(nskyi17) = les(i)
1855 fskyi17(nskyi17,1) = ffx
1856 fskyi17(nskyi17,2) = ffy
1857 fskyi17(nskyi17,3) = ffz
1858 fskyi17(nskyi17,4) = area(i)
1859 END IF
1860 ELSE ! partie complementaire SPMD
1861 IF(iparit == 0)THEN
1862 frotsfi(nint)%P(1,ies) = frotsfi(nint)%P(1,ies) + ffx
1863 frotsfi(nint)%P(2,ies) = frotsfi(nint)%P(2,ies) + ffy
1864 frotsfi(nint)%P(3,ies) = frotsfi(nint)%P(3,ies) + ffz
1865 frotsfi(nint)%P(4,ies) = frotsfi(nint)%P(4,ies) + area(i)
1866 ELSE
1867 nrskyi17=nrskyi17+1
1868 irskyi17(nrskyi17) = ies
1869 frskyi17(1,nrskyi17) = ffx
1870 frskyi17(2,nrskyi17) = ffy
1871 frskyi17(3,nrskyi17) = ffz
1872 frskyi17(4,nrskyi17) = area(i)
1873 END IF
1874 END IF
1875#include "lockoff.inc"
1876C-----------------------------------------------
1877C TH
1878C-----------------------------------------------
1879 fsav1 =fsav1 + fx
1880 fsav2 =fsav2 + fy
1881 fsav3 =fsav3 + fz
1882 fsav4 =fsav4 + fftx
1883 fsav5 =fsav5 + ffty
1884 fsav6 =fsav6 + fftz
1885C-----------------------------------------------
1886C force totale
1887C-----------------------------------------------
1888 fx = fx + fftx
1889 fy = fy + ffty
1890 fz = fz + fftz
1891 stif = stif + stifv
1892 econt = econt
1893 . + vx*fx+vy*fy+vz*fz
1894
1895 IF(les(i)>0)THEN
1896 IF(iparit == 0)THEN
1897#include "lockon.inc"
1898 DO ik=1,8
1899 j1 = iii(i,ik)
1900 a(1,j1) = a(1,j1) + ni_m(i,ik)*fx
1901 a(2,j1) = a(2,j1) + ni_m(i,ik)*fy
1902 a(3,j1) = a(3,j1) + ni_m(i,ik)*fz
1903 stifn(j1) = stifn(j1) + abs(ni_m(i,ik)*stif)
1904
1905 j2 = iiis(i,ik)
1906 a(1,j2) = a(1,j2) - ni_s(i,ik)*fx
1907 a(2,j2) = a(2,j2) - ni_s(i,ik)*fy
1908 a(3,j2) = a(3,j2) - ni_s(i,ik)*fz
1909 stifn(j2) = stifn(j2) + abs(ni_s(i,ik)*stif)
1910 ENDDO
1911#include "lockoff.inc"
1912 ELSE
1913#include "lockon.inc"
1914 niskyl = nisky
1915 nisky = nisky + 16
1916#include "lockoff.inc"
1917 IF (niskyl > lskyi) THEN
1918 CALL ancmsg(msgid=26,anmode=aninfo,
1919 . i1=noint)
1920 CALL arret(2)
1921 ENDIF
1922 DO ik=1,8
1923 niskyl = niskyl + 1
1924 fskyi(niskyl,1)=ni_m(i,ik)*fx
1925 fskyi(niskyl,2)=ni_m(i,ik)*fy
1926 fskyi(niskyl,3)=ni_m(i,ik)*fz
1927 fskyi(niskyl,4)=abs(ni_m(i,ik)*stif)
1928 isky(niskyl) = iii(i,ik)
1929 niskyl = niskyl + 1
1930 fskyi(niskyl,1)= - ni_s(i,ik)*fx
1931 fskyi(niskyl,2)= - ni_s(i,ik)*fy
1932 fskyi(niskyl,3)= - ni_s(i,ik)*fz
1933 fskyi(niskyl,4)= abs(ni_s(i,ik)*stif)
1934 isky(niskyl) = iiis(i,ik)
1935 ENDDO
1936 ENDIF
1937C
1938 IF(anim_v(4)+outp_v(4)+h3d_data%N_VECT_CONT >0.AND.
1939 . ((tt>=tanim .AND. tt<=tanim_stop).OR.tt>=toutp.OR.(tt>=h3d_data%TH3D.AND.tt<=h3d_data%TH3D_STOP).OR.
1940 . (manim>=4.AND.manim<=15).OR. h3d_data%MH3D /= 0))THEN
1941#include "lockon.inc"
1942 DO ik=1,8
1943 j1 = iii(i,ik)
1944 fcont(1,j1) =fcont(1,j1) + ni_m(i,ik)*fx
1945 fcont(2,j1) =fcont(2,j1) + ni_m(i,ik)*fy
1946 fcont(3,j1) =fcont(3,j1) + ni_m(i,ik)*fz
1947
1948 j2 = iiis(i,ik)
1949 fcont(1,j2) =fcont(1,j2) - ni_s(i,ik)*fx
1950 fcont(2,j2) =fcont(2,j2) - ni_s(i,ik)*fy
1951 fcont(3,j2) =fcont(3,j2) - ni_s(i,ik)*fz
1952 ENDDO
1953#include "lockoff.inc"
1954 ENDIF
1955 ELSE ! cas complementaire uniquement en SPMD
1956 IF(iparit == 0)THEN
1957#include "lockon.inc"
1958 DO ik=1,8
1959 j1 = iii(i,ik)
1960 a(1,j1) = a(1,j1) + ni_m(i,ik)*fx
1961 a(2,j1) = a(2,j1) + ni_m(i,ik)*fy
1962 a(3,j1) = a(3,j1) + ni_m(i,ik)*fz
1963 stifn(j1) = stifn(j1) + abs(ni_m(i,ik)*stif)
1964C
1965 j2 = iiis(i,ik)
1966 afi17(nint)%P(1,j2,ies) = afi17(nint)%P(1,j2,ies)
1967 - - ni_s(i,ik)*fx
1968 afi17(nint)%P(2,j2,ies) = afi17(nint)%P(2,j2,ies)
1969 - - ni_s(i,ik)*fy
1970 afi17(nint)%P(3,j2,ies) = afi17(nint)%P(3,j2,ies)
1971 - - ni_s(i,ik)*fz
1972 stnfi17(nint)%P(j2,ies) = stnfi17(nint)%P(j2,ies)
1973 + + abs(ni_s(i,ik)*stif)
1974 ENDDO
1975#include "lockoff.inc"
1976 ELSE
1977#include "lockon.inc"
1978 niskyl = nisky
1979 nisky = nisky + 8
1980 niskyfi=niskyfi+1
1981 niskyfil=niskyfi
1982#include "lockoff.inc"
1983 iskyfi(nint)%P(niskyfil) = ies
1984C
1985 IF (niskyl > lskyi) THEN
1986 CALL ancmsg(msgid=26,anmode=aninfo,
1987 . i1=noint)
1988 CALL arret(2)
1989 ENDIF
1990 IF (niskyfil > nlskyfi(nint)) THEN
1991 CALL ancmsg(msgid=26,anmode=aninfo,
1992 . i1=noint)
1993 CALL arret(2)
1994 ENDIF
1995C
1996 DO ik=1,8
1997 niskyl = niskyl + 1
1998 fskyi(niskyl,1)=ni_m(i,ik)*fx
1999 fskyi(niskyl,2)=ni_m(i,ik)*fy
2000 fskyi(niskyl,3)=ni_m(i,ik)*fz
2001 fskyi(niskyl,4)=abs(ni_m(i,ik)*stif)
2002 isky(niskyl) = iii(i,ik)
2003C
2004 fskyfi(nint)%P(1+(ik-1)*5,niskyfil)= - ni_s(i,ik)*fx
2005 fskyfi(nint)%P(2+(ik-1)*5,niskyfil)= - ni_s(i,ik)*fy
2006 fskyfi(nint)%P(3+(ik-1)*5,niskyfil)= - ni_s(i,ik)*fz
2007 fskyfi(nint)%P(4+(ik-1)*5,niskyfil)= abs(ni_s(i,ik)*stif)
2008 fskyfi(nint)%P(5+(ik-1)*5,niskyfil)= iiis(i,ik)
2009 ENDDO
2010 ENDIF
2011C
2012 IF(anim_v(4)+outp_v(4)+h3d_data%N_VECT_CONT >0.AND.
2013 . ((tt>=tanim .AND. tt<=tanim_stop).OR.tt>=toutp.OR.(tt>=h3d_data%TH3D.AND.tt<=h3d_data%TH3D_STOP).OR.
2014 . (manim>=4.AND.manim<=15).OR. h3d_data%MH3D /= 0))THEN
2015#include "lockon.inc"
2016 DO ik=1,8
2017 j1 = iii(i,ik)
2018 fcont(1,j1) =fcont(1,j1) + ni_m(i,ik)*fx
2019 fcont(2,j1) =fcont(2,j1) + ni_m(i,ik)*fy
2020 fcont(3,j1) =fcont(3,j1) + ni_m(i,ik)*fz
2021 ENDDO
2022#include "lockoff.inc"
2023 ENDIF
2024 END IF
2025 ENDIF
2026 ENDDO
2027C
2028C-----------------------------------------------
2029C TH
2030C-----------------------------------------------
2031#include "lockon.inc"
2032 fsav(1) = fsav(1) + fsav1*dt12
2033 fsav(2) = fsav(2) + fsav2*dt12
2034 fsav(3) = fsav(3) + fsav3*dt12
2035 fsav(4) = fsav(4) + fsav4*dt12
2036 fsav(5) = fsav(5) + fsav5*dt12
2037 fsav(6) = fsav(6) + fsav6*dt12
2038 econtv = econtv + dt1*econvt
2039 econt = econt + dt1*econtt
2040#include "lockoff.inc"
2041C
2042 RETURN
area
subroutine area(d1, x, x2, y, y2, eint, stif0)
Definition hm_read_prop32.F:567
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
tri7box::stnfi17
type(real_pointer2), dimension(:), allocatable stnfi17
Definition tri7box.F:467
tri7box::afi17
type(real_pointer3), dimension(:), allocatable afi17
Definition tri7box.F:470
tri7box::frotsfi
type(real_pointer2), dimension(:), allocatable frotsfi
Definition tri7box.F:467
tri7box::vfi17
type(real_pointer3), dimension(:), allocatable vfi17
Definition tri7box.F:470
tri7box::ksfi
type(real_pointer2), dimension(:), allocatable ksfi
Definition tri7box.F:467
tri7box::fskyfi
type(real_pointer2), dimension(:), allocatable fskyfi
Definition tri7box.F:459
tri7box::iskyfi
type(int_pointer), dimension(:), allocatable iskyfi
Definition tri7box.F:480
tri7box::nlskyfi
integer, dimension(:), allocatable nlskyfi
Definition tri7box.F:512

◆ i17lll_ass()

subroutine i17lll_ass ( integer nskyi17,
integer, dimension(*) iskyi17,
fskyi17,
frots,
integer nmes,
integer lskyi17 )

Definition at line 2052 of file i17for3.F.

2054C-----------------------------------------------
2055C I m p l i c i t T y p e s
2056C-----------------------------------------------
2057#include "implicit_f.inc"
2058C-----------------------------------------------
2059C D u m m y A r g u m e n t s
2060C-----------------------------------------------
2061 INTEGER NSKYI17, NMES, LSKYI17, ISKYI17(*)
2062C REAL
2063 my_real
2064 . fskyi17(lskyi17,4),frots(7,*)
2065C-----------------------------------------------
2066C L o c a l V a r i a b l e s
2067C-----------------------------------------------
2068 INTEGER I, J, JJ1, JJ2, K, KK, L, N, NN, IDIFF,
2069 . ADSKYI(0:NMES+1)
2070 my_real
2071 . ff, fskyit(nskyi17) , bid
2072C-----------------------------------------------
2073 bid = zero
2074 DO n = 1, nmes+1
2075 adskyi(n) = 0
2076 END DO
2077 DO i=1,nskyi17
2078 n = iskyi17(i)+1
2079 adskyi(n) = adskyi(n)+1
2080 END DO
2081 adskyi(0) = 1
2082 adskyi(1) = 1
2083 DO n = 1, nmes
2084 nn = n+1
2085 adskyi(nn) = adskyi(nn) + adskyi(n)
2086 END DO
2087C
2088 DO i=1,nskyi17
2089 n = iskyi17(i)
2090 j = adskyi(n)
2091 iskyi17(i) = j
2092 adskyi(n) = adskyi(n) + 1
2093 END DO
2094C
2095 DO l = 1, 4
2096 DO i=1,nskyi17
2097 j = iskyi17(i)
2098 fskyit(j) = fskyi17(i,l)
2099 END DO
2100 DO i=1,nskyi17
2101 fskyi17(i,l) = fskyit(i)
2102 END DO
2103 ENDDO
2104C
2105 DO i=1,nskyi17
2106 iskyi17(i) = 0
2107 END DO
2108C
2109 DO n = 1, nmes
2110 nn = n-1
2111 idiff = adskyi(n)-adskyi(nn)
2112 IF(idiff==1)THEN
2113 k = adskyi(nn)
2114 frots(1,n) = frots(1,n) + fskyi17(k,1)
2115 frots(2,n) = frots(2,n) + fskyi17(k,2)
2116 frots(3,n) = frots(3,n) + fskyi17(k,3)
2117 frots(4,n) = frots(4,n) + fskyi17(k,4)
2118 ELSEIF(idiff>1.AND.idiff<20)THEN
2119 jj1 = adskyi(nn)
2120 jj2 = adskyi(n)-1
2121 DO k=jj1,jj2-1
2122 DO kk=k+1,jj2
2123 DO l=1,4
2124 IF(fskyi17(kk,l)>fskyi17(k,l))THEN
2125 ff = fskyi17(kk,l)
2126 fskyi17(kk,l) = fskyi17(k,l)
2127 fskyi17(k,l) = ff
2128 END IF
2129 END DO
2130 END DO
2131 END DO
2132 DO k=jj1,jj2
2133 frots(1,n)= frots(1,n)+ fskyi17(k,1)
2134 frots(2,n)= frots(2,n)+ fskyi17(k,2)
2135 frots(3,n)= frots(3,n)+ fskyi17(k,3)
2136 frots(4,n)= frots(4,n)+ fskyi17(k,4)
2137 ENDDO
2138 ELSEIF(idiff>=20)THEN
2139C qsort necessaire si > 20
2140 jj1 = adskyi(nn)
2141 jj2 = adskyi(n)-1
2142 CALL ass2sort(fskyi17,jj1,jj2,fskyit,4)
2143 DO k=jj1,jj2
2144 frots(1,n)= frots(1,n)+ fskyi17(k,1)
2145 frots(2,n)= frots(2,n)+ fskyi17(k,2)
2146 frots(3,n)= frots(3,n)+ fskyi17(k,3)
2147 frots(4,n)= frots(4,n)+ fskyi17(k,4)
2148 ENDDO
2149 END IF
2150 END DO
2151C
2152 RETURN
ass2sort
subroutine ass2sort(fskyi, jj1, jj2, fskyt, nfsk)
Definition ass2sort.F:34

◆ i17lll_pena()

subroutine i17lll_pena ( integer llt,
v,
stifn,
xx,
fric,
yy,
zz,
integer, dimension(mvsiz,17) iii,
fskyi,
integer, dimension(*) isky,
a,
x,
integer itied,
integer nint,
xxs,
yys,
zzs,
integer, dimension(mvsiz,16) iiis,
vit_min,
integer, dimension(*) le,
integer, dimension(*) les,
integer ie_min,
integer ies_min,
integer itask,
frotm,
frots,
km,
ks,
fsav,
fcont,
ms,
integer niskyfi,
integer noint,
type(h3d_database) h3d_data )

Definition at line 280 of file i17for3.F.

288C-----------------------------------------------
289C M o d u l e s
290C-----------------------------------------------
291 USE h3d_mod
292 USE anim_mod
293C-----------------------------------------------
294C I m p l i c i t T y p e s
295C-----------------------------------------------
296#include "implicit_f.inc"
297#include "comlock.inc"
298C-----------------------------------------------
299C G l o b a l P a r a m e t e r s
300C-----------------------------------------------
301#include "mvsiz_p.inc"
302C-----------------------------------------------
303C D u m m y A r g u m e n t s
304C-----------------------------------------------
305 INTEGER LLT,ITIED,NINT ,IE_MIN,IES_MIN,NISKYFI,NOINT
306 INTEGER ITASK,
307 . III(MVSIZ,17),IIIS(MVSIZ,16),LE(*) ,LES(*), ISKY(*)
308C REAL
309 my_real
310 . v(3,*),a(3,*),km(*),ks(*),fric(*), ms(*)
311 my_real
312 . xx(mvsiz,17),yy(mvsiz,17),zz(mvsiz,17),x(3,*),
313 . xxs(mvsiz,16) ,yys(mvsiz,16) ,zzs(mvsiz,16) ,vit_min,
314 . stifn(*) ,fskyi ,frotm(7,*),frots(7,*),fsav(*) ,fcont(3,*)
315 TYPE(H3D_DATABASE) :: H3D_DATA
316C-----------------------------------------------
317C L o c a l V a r i a b l e s
318C-----------------------------------------------
319 INTEGER I,J,IK,NK,I1,I2,I3,I4,IAD,ICON,
320 + ICONT(MVSIZ,4)
321 my_real
322 . vx,vy,vz,vn,aa,vv,pene
323 my_real
324 . r_cm(mvsiz),t_cm(mvsiz),s_cm(mvsiz),si_s(mvsiz,8),
325 . ri_s(mvsiz,8),ti_s(mvsiz,8),
326 . nx(mvsiz), ny(mvsiz), nz(mvsiz),
327 . ni_m(mvsiz,17),ni_s(mvsiz,8) ,
328 . r_1s(mvsiz) ,r_2s(mvsiz) ,t_1s(mvsiz) ,t_2s(mvsiz),
329 . r_1m(mvsiz) ,r_2m(mvsiz) ,t_1m(mvsiz) ,t_2m(mvsiz),
330 . r_3m(mvsiz) ,r_4m(mvsiz) ,t_3m(mvsiz) ,t_4m(mvsiz),
331 . r_cs(mvsiz) ,s_cs(mvsiz) ,t_cs(mvsiz) ,vit(mvsiz),
332 . r_s(mvsiz) ,t_s(mvsiz) ,area(mvsiz) ,area_tot(mvsiz),
333 . area_el(mvsiz)
334C-----------------------------------------------
335C calcul de r,s,t
336C-----------------------------------------------
337C noeud centrale face 1 2 3 4
338C-----------------------------------------------
339 DO i=1,llt
340 xx(i,17) = half *(xxs(i,5) +xxs(i,6) +xxs(i,7) +xxs(i,8))
341 . - fourth*(xxs(i,1) +xxs(i,2) +xxs(i,3) +xxs(i,4))
342 yy(i,17) = half *(yys(i,5) +yys(i,6) +yys(i,7) +yys(i,8))
343 . - fourth*(yys(i,1) +yys(i,2) +yys(i,3) +yys(i,4))
344 zz(i,17) = half *(zzs(i,5) +zzs(i,6) +zzs(i,7) +zzs(i,8))
345 . - fourth*(zzs(i,1) +zzs(i,2) +zzs(i,3) +zzs(i,4))
346 ENDDO
347 CALL i17rst(llt ,ri_s(1,1),si_s(1,1),ti_s(1,1),ni_m ,
348 2 xx ,yy ,zz )
349C-----------------------------------------------
350C noeud centrale face 5 6 7 8
351C-----------------------------------------------
352 DO i=1,llt
353 xx(i,17) = half *(xxs(i,13)+xxs(i,14)+xxs(i,15)+xxs(i,16))
354 . - fourth*(xxs(i,9) +xxs(i,10)+xxs(i,11)+xxs(i,12))
355 yy(i,17) = half *(yys(i,13)+yys(i,14)+yys(i,15)+yys(i,16))
356 . - fourth*(yys(i,9) +yys(i,10)+yys(i,11)+yys(i,12))
357 zz(i,17) = half *(zzs(i,13)+zzs(i,14)+zzs(i,15)+zzs(i,16))
358 . - fourth*(zzs(i,9) +zzs(i,10)+zzs(i,11)+zzs(i,12))
359 ENDDO
360 CALL i17rst(llt ,ri_s(1,2),si_s(1,2),ti_s(1,2),ni_m ,
361 2 xx ,yy ,zz )
362C-----------------------------------------------
363C choix face 1 2 3 4 face 5 6 7 8 cote second
364C-----------------------------------------------
365 DO i=1,llt
366 IF(abs(si_s(i,1))<=abs(si_s(i,2)))THEN
367C face 1 2 3 4
368 s_cs(i) = -one
369 ELSE
370C face 5 6 7 8
371 s_cs(i) = one
372 iiis(i,1) = iiis(i,12)
373 iiis(i,2) = iiis(i,11)
374 iiis(i,3) = iiis(i,10)
375 iiis(i,4) = iiis(i, 9)
376 iiis(i,5) = iiis(i,15)
377 iiis(i,6) = iiis(i,14)
378 iiis(i,7) = iiis(i,13)
379 iiis(i,8) = iiis(i,16)
380C
381 xxs(i,1) = xxs(i,12)
382 xxs(i,2) = xxs(i,11)
383 xxs(i,3) = xxs(i,10)
384 xxs(i,4) = xxs(i, 9)
385 xxs(i,5) = xxs(i,15)
386 xxs(i,6) = xxs(i,14)
387 xxs(i,7) = xxs(i,13)
388 xxs(i,8) = xxs(i,16)
389C
390 yys(i,1) = yys(i,12)
391 yys(i,2) = yys(i,11)
392 yys(i,3) = yys(i,10)
393 yys(i,4) = yys(i, 9)
394 yys(i,5) = yys(i,15)
395 yys(i,6) = yys(i,14)
396 yys(i,7) = yys(i,13)
397 yys(i,8) = yys(i,16)
398C
399 zzs(i,1) = zzs(i,12)
400 zzs(i,2) = zzs(i,11)
401 zzs(i,3) = zzs(i,10)
402 zzs(i,4) = zzs(i, 9)
403 zzs(i,5) = zzs(i,15)
404 zzs(i,6) = zzs(i,14)
405 zzs(i,7) = zzs(i,13)
406 zzs(i,8) = zzs(i,16)
407C
408 ENDIF
409 ENDDO
410C-----------------------------------------------
411C calcul de SI_S=s(relatif element main)
412c des 8 noeuds de la face de l'element second
413C-----------------------------------------------
414 DO j=1,8
415 DO i=1,llt
416 xx(i,17) = xxs(i,j)
417 yy(i,17) = yys(i,j)
418 zz(i,17) = zzs(i,j)
419 ENDDO
420 CALL i17rst(llt ,ri_s(1,j),si_s(1,j),ti_s(1,j),ni_m ,
421 2 xx ,yy ,zz )
422 ENDDO
423C-----------------------------------------------
424C calcul du point de distance min
425C dSI_S/dr = 0 dSI_S/dt = 0 (r,t de l'element second)
426C-----------------------------------------------
427c
428c il est possible de calculer la surface effective
429c de contact pour ponderer la force
430c => pas de discontinuit si contact sur 2 elements
431c
432c a faire dans I17MINI
433c
434c possible aussi de calculer la pression de Hertz
435c
436 CALL i17mini_pena(
437 1 llt ,r_cs ,s_cs ,t_cs ,ri_s ,si_s ,
438 2 ti_s ,ni_s ,xxs ,yys ,zzs ,xx ,
439 3 yy ,zz ,r_cm ,s_cm ,t_cm ,nx ,
440 4 ny ,nz ,r_1s ,r_2s ,t_1s ,t_2s ,
441 5 r_1m ,r_2m ,r_3m ,r_4m ,t_1m ,t_2m ,
442 6 t_3m ,t_4m ,icont ,area ,area_tot,area_el)
443C-----------------------------------------------
444C choix face 1 2 3 4 face 5 6 7 8 cote main
445C-----------------------------------------------
446 DO i=1,llt
447 IF(s_cm(i)<=0.)THEN
448C face 1 2 3 4
449 iii(i,5) = iii(i,9)
450 iii(i,6) = iii(i,10)
451 iii(i,7) = iii(i,11)
452 iii(i,8) = iii(i,12)
453 ELSE
454C face 5 6 7 8
455 iii(i,1) = iii(i,5)
456 iii(i,2) = iii(i,6)
457 iii(i,3) = iii(i,7)
458 iii(i,4) = iii(i,8)
459 iii(i,5) = iii(i,13)
460 iii(i,6) = iii(i,14)
461 iii(i,7) = iii(i,15)
462 iii(i,8) = iii(i,16)
463 ENDIF
464 ENDDO
465 CALL i17vit_pena(
466 1 llt ,nint ,v ,a ,iii ,iiis ,
467 2 ni_m ,ni_s ,nx ,ny ,nz ,vit ,
468 3 icont(1,1),r_cm ,t_cm ,r_cs ,t_cs ,s_cm ,
469 4 les )
470C-----------------------------------------------------------------------
471C calcul du contact
472C-----------------------------------------------------------------------
473 CALL i17lll4_pena(
474 1 llt ,itied ,nint ,v ,a ,fric ,
475 3 iii ,iiis ,ni_m ,ni_s ,s_cm ,s_cs ,
476 4 nx ,ny ,nz ,vit ,le ,les ,
477 5 icont(1,1),r_cm ,t_cm ,r_cs ,t_cs ,xx ,
478 6 yy ,zz ,stifn ,fskyi ,isky ,
479 7 frotm ,frots ,area ,area_tot,km ,ks ,
480 8 fsav ,fcont ,ms ,area_el ,niskyfi ,noint ,
481 9 h3d_data )
482 RETURN
i17lll4_pena
subroutine i17lll4_pena(llt, itied, nint, v, a, fric, iii, iiis, ni_m, ni_s, s_cm, s_cs, nx, ny, nz, vit, le, les, icont, rm, tm, rs, ts, xx, yy, zz, stifn, fskyi, isky, frotm, frots, area, area_tot, km, ks, fsav, fcont, ms, area_el, niskyfi, noint, h3d_data)
Definition i17for3.F:1507
i17mini_pena
subroutine i17mini_pena(llt, r_cs, s_cs, t_cs, ri_s, si_s, ti_s, ni_s, xxs, yys, zzs, xx, yy, zz, r_cm, s_cm, t_cm, nx, ny, nz, r_1s, r_2s, t_1s, t_2s, r_1m, r_2m, r_3m, r_4m, t_1m, t_2m, t_3m, t_4m, icont, area, area_tot, area_el)
Definition i17for3.F:502
i17vit_pena
subroutine i17vit_pena(llt, nint, v, a, iii, iiis, ni_m, ni_s, nx, ny, nz, vit, icont, rm, tm, rs, ts, sm, les)
Definition i17for3.F:1405
i17rst
subroutine i17rst(llt, r, s, t, ni, xx, yy, zz)
Definition i17lagm.F:772

◆ i17mini_pena()

subroutine i17mini_pena ( integer llt,
r_cs,
s_cs,
t_cs,
ri_s,
si_s,
ti_s,
ni_s,
xxs,
yys,
zzs,
xx,
yy,
zz,
r_cm,
s_cm,
t_cm,
nx,
ny,
nz,
r_1s,
r_2s,
t_1s,
t_2s,
r_1m,
r_2m,
r_3m,
r_4m,
t_1m,
t_2m,
t_3m,
t_4m,
integer, dimension(mvsiz,4) icont,
area,
area_tot,
area_el )

Definition at line 495 of file i17for3.F.

502C-----------------------------------------------
503C I m p l i c i t T y p e s
504C-----------------------------------------------
505#include "implicit_f.inc"
506C-----------------------------------------------
507C G l o b a l P a r a m e t e r s
508C-----------------------------------------------
509#include "mvsiz_p.inc"
510C-----------------------------------------------
511C D u m m y A r g u m e n t s
512C-----------------------------------------------
513 INTEGER LLT,
514 + ICONT(MVSIZ,4)
515C REAL
516 my_real
517 + si_s(mvsiz,*),ni_s(mvsiz,*),ri_s(mvsiz,*),ti_s(mvsiz,*),
518 + xxs(mvsiz,*) ,yys(mvsiz,*) ,zzs(mvsiz,*) ,
519 + xx(mvsiz,*) ,yy(mvsiz,*) ,zz(mvsiz,*) ,
520 + r_cm(mvsiz) ,s_cm(mvsiz) ,t_cm(mvsiz),
521 + r_cs(mvsiz) ,s_cs(mvsiz) ,t_cs(mvsiz),
522 + r_1s(mvsiz) ,r_2s(mvsiz) ,t_1s(mvsiz) ,t_2s(mvsiz),
523 + r_1m(mvsiz) ,r_2m(mvsiz) ,t_1m(mvsiz) ,t_2m(mvsiz),
524 + r_3m(mvsiz) ,r_4m(mvsiz) ,t_3m(mvsiz) ,t_4m(mvsiz),
525 + nx(mvsiz) ,ny(mvsiz) ,nz(mvsiz) ,area(mvsiz) ,area_tot(mvsiz),
526 + area_el(mvsiz)
527C-----------------------------------------------
528C L o c a l V a r i a b l e s
529C-----------------------------------------------
530 INTEGER I,ITER,NITERMAX,IR,IT
531 my_real
532 + a1(mvsiz),a2(mvsiz),a3(mvsiz),a4(mvsiz),a5(mvsiz),
533 + b1(mvsiz),b2(mvsiz),b3(mvsiz),b4(mvsiz),b5(mvsiz),
534 + c1(mvsiz),c2(mvsiz),c3(mvsiz),azero(mvsiz),
535 + f1,f2,f3,f4,f5,f6,f7,f8,
536 + bb1,bb2,bb3,
537 + cc1,cc2,cc3,dd1,dd2,dd3,dd,d,
538 + a0,ab,ba,a4r,b4t,a5t,b5r,eps,
539 + xa,ya,za,xb,yb,zb,xc,yc,zc,aaa,unpeps,
540 + p,rm,tm,sm,pp,rr,tt,aa,rt(9),as(9)
541 DATA rt/-1.,-0.75,-0.5,-0.25,0.0,0.25,0.5,0.75,1./
542 DATA as/0.0625,0.125,0.125,0.125,0.125,0.125,0.125,0.125,0.0625/
543C-----------------------------------------------
544C ro = r ri to = t ti
545C
546C i=1,4
547C ri=+-1 ti=+-1
548C Ni = 1/4 (1+ro)(1+to)(ro+to-1)
549C Ni = 1/4 (r^2+r^2to+t^2+roto+rot^2-1)
550C Ni = 1/4 (r^2(1+to)+t^2(1+ro)+roto-1)
551C dNi/dr = ri/4 (1+to)(2ro+to)
552C dNi/dt = ti/4 (1+ro)(2to+ro)
553C dNi/dr = 1/4 (2 r + ri ti t + 2ti rt + ri t^2)
554C dNi/dt = 1/4 (2 t + ri ti r + 2ri rt + ti r^2)
555C
556C i=6;8
557C ri=0 ti=+-1
558C Ni = 1/2 (1-r^2)(1+to)
559C dNi/dr = -r (1+to)
560C dNi/dt = ti/2 (1-r^2)
561C dNi/dr = 1/4 (-4 r - 4ti r t)
562C dNi/dt = 1/4 (2ti - 2ti r^2 )
563C
564C
565C i=5;7
566C ri=+-1 ti=0
567C Ni = 1/2 (1-t^2)(1+ro)
568C dNi/dr = ri/2 (1-t^2)
569C dNi/dt = -t (1+ro)
570C dNi/dr = 1/4 (2ri - 2ri t^2 )
571C dNi/dt = 1/4 (-4 t - 4ri rt )
572C-----------------------------------------------
573C df/dr = Somme( fi dNi/dr )
574C df/dt = Somme( fi dNi/dt )
575C
576C df/dr = A1 + A2 r + A3 t + A4 rt + A5 t^2
577C df/dt = B1 + B2 r + B3 t + B4 rt + B5 r^2
578C
579C A1 = ( -f5 + f7 )/2
580C A2 = ( f1 + f2 + f3 + f4 - 2 f6 - 2 f8)/2
581C A3 = ( f1 - f2 + f3 - f4 )/4
582C A4 = (-f1 + f2 + f3 - f4 - 2 f6 + 2 f8)/2
583C A5 = (-f1 - f2 + f3 + f4 + 2 f5 - 2 f7 )/4
584C
585C B1 = ( f6 - f8)/2
586C B2 = ( f1 - f2 + f3 - f4 )/4
587C B3 = ( f1 + f2 + f3 + f4 - 2 f5 - 2 f7 )/2
588C B4 = (-f1 - f2 + f3 + f4 + 2 f5 - 2 f7 )/2
589C B5 = (-f1 + f2 + f3 - f4 - 2 f6 + 2 f8)/4
590C
591C df/dr = A1 + A2 r + A3 t + A4 rt + A5 t^2 = 0
592C df/dt = B1 + B2 r + B3 t + B4 rt + B5 r^2 = 0
593C r = -(A1 + A3 t + A5 t^2 ) / (A2 + A4 t)
594C t = -(B1 + B2 r + B5 r^2 ) / (B3 + B4 r)
595c
596c
597c r
598C ^
599c |
600C . |7
601C 4 O-----------------O-----------------O 3
602C | . | |
603C | . | |
604C | | C |
605C | . r+------+df/dt=0 |
606C | . | |df/dr=0 |
607C | | | |
608C | | | |
609C | | | |6
610C 8 0 +------+----------0----> t
611C | t |
612C | |
613C | |
614C | |
615C | |
616C | |
617C | |
618C |' |
619C O-----------------O-----------------O
620C 1 5 2
621C
622C-----------------------------------------------
623C
624C second. main
625C (1): <R1S , TCS> <R1M , T1M>
626C (2): <R2S , TCS> <R2M , T2M>
627C (3): <RCS , T1S> <R3M , T3M>
628C (4): <RCS , T2S> <R4M , T4M>
629C
630C
631c O
632c / \
633c / _____
634c / / \ \
635c / / \ +r2s
636c / / \ / \
637c / / r2s'+ \
638c / / /|\ |<-----Area_tot
639c rm O / / | \ |
640C ^ / / / | \| rs
641c |/t1s/ / | / ^
642C / + / (2)| / \ /
643C O----------------/O---o------------oO/ O
644C | / | /(3)\:::/:::::::| / \
645C | / |/:::::\:/:::::::/| / \
646C | /rcm+-------+rcs,tcs/ | / \
647C | / /|::::::/|\:::::/ |/ \
648C | / /:|:::::/:|:\:::/ / \
649C | O /::|::::/::|::\:/ /| \
650C | \ /:::|:::/:::|:::o(4) | \
651C | \::::|::/::::|::/t2s | \
652C O |\:::+-/-----+-/-------O----> tm \
653C | | \:::/::::tcm/ \| \
654C | | \:/:::::::/ \ O
655C | \‍(1)o:::::::/ |\ /
656C | \ / \:::::/ | \ /
657C | + \:::/ | \ /
658C | r1s\___\_/ | \ /
659C | \ | \ /
660C |main \ | \ /
661C O-----------------O----O------------O \ /
662C \ O
663C \ / \
664C \ / \
665C :::: Area \ / v
666C :::: commune \ / ts
667C :::: main Sec \ /
668C \second/
669C \ /
670C \ /
671C \ /
672C O
673c
674c Area_tot = pi/4 | (r1s,r2s) x (t1s,t2s) |
675c Area ~= pi/4 | ((1),(2)) x ((3),(4)) |
676C-----------------------------------------------
677C
678c
679 nitermax = 5
680C
681 DO i=1,llt
682 azero(i) = zero
683 d = si_s(i,1)*si_s(i,1)+si_s(i,2)*si_s(i,2)
684 + + si_s(i,3)*si_s(i,3)+si_s(i,4)*si_s(i,4)
685 + + si_s(i,5)*si_s(i,5)+si_s(i,6)*si_s(i,6)
686 + + si_s(i,7)*si_s(i,7)+si_s(i,8)*si_s(i,8)
687 d = 1./max(em20,sqrt(d))
688 f1 = d * si_s(i,1)
689 f2 = d * si_s(i,2)
690 f3 = d * si_s(i,3)
691 f4 = d * si_s(i,4)
692 f5 = d * si_s(i,5)
693 f6 = d * si_s(i,6)
694 f7 = d * si_s(i,7)
695 f8 = d * si_s(i,8)
696 a0 = ( f1 + f2 + f3 + f4 )*half
697 ab = f5 - f7
698 ba = f6 - f8
699C
700 a1(i) = -ab * half
701 a2(i) = a0 - f6 - f8
702 a3(i) = ( f1 - f2 + f3 - f4 )*fourth
703 a4(i) = (-f1 + f2 + f3 - f4 )*half - ba
704 a5(i) = (-f1 - f2 + f3 + f4 )*fourth - a1(i)
705C
706 b1(i) = ba*half
707 b2(i) = a3(i)
708 b3(i) = a0 - f5 - f7
709 b4(i) = (-f1 - f2 + f3 + f4 )*half + ab
710 b5(i) = (-f1 + f2 + f3 - f4 )*fourth - b1(i)
711c
712 r_cs(i) = zero
713 t_cs(i) = zero
714 ENDDO
715c------------------------------------------------
716c newton iter : lineari_sation en r et t
717c------------------------------------------------
718C fr = df/dr = A1 + A2 r + A3 t + A4 t r + A5 t t = 0
719C ft = df/dt = B1 + B2 r + B3 t + B4 r t + B5 r r = 0
720c
721C fr = fr_ + dfr/dr dr + dfr/dt dt
722C ft = ft_ + dft/dr dr + dft/dt dt
723c
724C fr = fr_ + dfr/dr (r-r_) + dfr/dt (t-t_) = 0
725C ft = ft_ + dft/dr (r-r_) + dft/dt (t-t_) = 0
726c
727C dfr/dr = A2 + A4 t_ = C2
728C dfr/dt = A3 + A4 r_ + 2 A5 t_ = C3
729C dft/dr = B2 + B4 t_ + 2 B5 r_ = D2
730C dft/dt = B3 + B4 r_ = D3
731c
732c C1 = A1 - A4 r_ t_ - A5 t_^2
733c C2 = A2 + A4 t_
734c C3 = A3 + A4 r_ + 2 A5 t_
735c
736c D1 = B1 - B4 r_ t_ - B5 r_^2
737c D2 = B2 + B4 t_ + 2 B5 r_
738c D3 = B3 + B4 r_
739c
740C fr = C1 + C2 r + C3 t = 0
741C ft = D1 + D2 r + D3 t = 0
742c
743C r = (C3 D1 - D3 C1) / (D3 C2 - C3 D2)
744C t = (D2 C1 - C2 D1) / (D3 C2 - C3 D2)
745c------------------------------------------------
746 DO iter=1,nitermax
747 DO i=1,llt
748 a4r = a4(i) * r_cs(i)
749 a5t = a5(i) * t_cs(i)
750 b4t = b4(i) * t_cs(i)
751 b5r = b5(i) * r_cs(i)
752 cc1 = a1(i) -(a4r + a5t) * t_cs(i)
753 cc2 = a2(i) + a4(i) * t_cs(i)
754 cc3 = a3(i) + a4r + a5t + a5t
755 dd1 = b1(i) -(b4t + b5r )* r_cs(i)
756 dd2 = b2(i) + b4t + b5r + b5r
757 dd3 = b3(i) + b4(i) * r_cs(i)
758 d = dd3 * cc2 - cc3 * dd2
759 IF(abs(d)<em20)THEN
760 cc2 = cc2 + em10
761 dd3 = dd3 + em10
762 d = dd3 * cc2 - cc3 * dd2
763 ENDIF
764 r_cs(i) = (cc3 * dd1 - dd3 * cc1) / d
765 t_cs(i) = (dd2 * cc1 - cc2 * dd1) / d
766 r_cs(i) = max(-one,min(one,r_cs(i)))
767 t_cs(i) = max(-one,min(one,t_cs(i)))
768c dfdr = A1 + A2 * r + A3 * t + A4 * t * r + A5 * t * t
769c dfdt = B1 + B2 * r + B3 * t + B4 * r * t + B5 * r * r
770 ENDDO
771 ENDDO
772C-----------------------------------------------------------------------
773C calcul de 4 points de l'iso S_M = +- 1
774C
775C second. main
776C 1: <R_1S , T_CS> <R_1M , T_1M>
777C 2: <R_2S , T_CS> <R_2M , T_2M>
778C 3: <R_CS , T_1S> <R_3M , T_3M>
779C 4: <R_CS , T_2S> <R_4M , T_4M>
780C
781C on borne r et t a +- 1 sur l'element second. :
782C -1 < R_1S < 1 et -1 < R_2S < 1
783C -1 < T_1S < 1 et -1 < T_2S < 1
784c
785C on borne r et t a +- 1 sur l'element main :
786C -1 < R_M < 1 et -1 < T_M < 1
787C-----------------------------------------------------------------------
788 CALL i17abc(llt ,si_s,r_cs,t_cs,
789 + b1 ,b2 ,b3 ,c1 ,c2 ,c3 )
790C
791 DO i=1,llt
792 s_cm(i) = c1(i) + (c2(i) + c3(i)*t_cs(i))*t_cs(i)
793 s_cm(i) = b1(i) + (b2(i) + b3(i)*r_cs(i))*r_cs(i)
794 ENDDO
795c
796 DO i=1,llt
797 IF(s_cm(i)>zero)THEN
798 b1(i) = b1(i) - one
799 c1(i) = c1(i) - one
800 ELSE
801 b1(i) = b1(i) + one
802 c1(i) = c1(i) + one
803 ENDIF
804 ENDDO
805c
806C f = B1 + B2 r + B3 r^2
807c
808 CALL i17racine_pena(llt,b3,b2,b1,r_1s,r_2s)
809C f = C1 + C2 t + C3 t^2
810C
811 CALL i17racine_pena(llt,c3,c2,c1,t_1s,t_2s)
812C-----------------------------------------------
813C calcul de la surface totale second
814C-----------------------------------------------
815 CALL i17surf(llt ,r_1s ,r_2s ,r_cs ,r_cs ,
816 2 t_cs ,t_cs ,t_1s ,t_2s ,area_tot,
817 3 xxs ,yys ,zzs ,azero )
818c-----------------------------------------------------------------------
819c borne surface totale la surface de l'second
820c-----------------------------------------------------------------------
821 DO i=1,llt
822 xa = xxs(i,7)-xxs(i,5)
823 ya = yys(i,7)-yys(i,5)
824 za = zzs(i,7)-zzs(i,5)
825 xb = xxs(i,8)-xxs(i,6)
826 yb = yys(i,8)-yys(i,6)
827 zb = zzs(i,8)-zzs(i,6)
828 xc = ya*zb - za*yb
829 yc = za*xb - xa*zb
830 zc = xa*yb - ya*xb
831c AAA = PI*FOURTH*ABS(XC*NX(I) + YC*NY(I) + ZC*NZ(I))
832 aaa = pi*fourth*sqrt(xc*xc+yc*yc+zc*zc)
833 area_el(i) = aaa
834 area_tot(i) = min(area_tot(i),aaa)
835c IF(AREA_TOT(I) == ZERO)AREA_TOT(I) = AAA
836 ENDDO
837c-----------------------------------------------------------------------
838c decoupe de l'second 9x9
839c calcul surface commune
840c-----------------------------------------------------------------------
841 unpeps = one+em02
842 DO i=1,llt
843 area(i) = zero
844 pp = zero
845 rr = zero
846 tt = zero
847 aa = zero
848 DO ir=1,9
849 DO it=1,9
850 CALL i17abc_pena(i,ri_s,rt(ir),rt(it),
851 + bb1 ,bb2 ,bb3 ,cc1 ,cc2 ,cc3 )
852 rm = bb1 + (bb2 + bb3*rt(ir))*rt(ir)
853 IF(rm >= -unpeps.and.rm <= unpeps)THEN
854 CALL i17abc_pena(i,ti_s,rt(ir),rt(it),
855 + bb1 ,bb2 ,bb3 ,cc1 ,cc2 ,cc3 )
856 tm = cc1 + (cc2 + cc3*rt(it))*rt(it)
857 IF(tm >= -unpeps.and.tm <= unpeps)THEN
858 CALL i17abc_pena(i,si_s,rt(ir),rt(it),
859 + bb1 ,bb2 ,bb3 ,cc1 ,cc2 ,cc3 )
860 sm = cc1 + (cc2 + cc3*rt(it))*rt(it)
861 p = one - abs(sm)
862 IF(p > zero)THEN
863 aa = aa + as(ir)*as(it)
864 p = p * as(ir) * as(it)
865 pp = pp + p
866 rr = rr + p*rt(ir)
867 tt = tt + p*rt(it)
868 ENDIF
869 ENDIF
870 ENDIF
871 ENDDO
872 ENDDO
873 icont(i,1) = 0
874 IF(aa > zero)THEN
875 r_cs(i) = rr / pp
876 t_cs(i) = tt / pp
877 CALL i17abc_pena(i,ri_s,r_cs(i),t_cs(i),
878 + bb1 ,bb2 ,bb3 ,cc1 ,cc2 ,cc3 )
879 r_cm(i) = bb1 + (bb2 + bb3*r_cs(i))*r_cs(i)
880 CALL i17abc_pena(i,ti_s,r_cs(i),t_cs(i),
881 + bb1 ,bb2 ,bb3 ,cc1 ,cc2 ,cc3 )
882 t_cm(i) = cc1 + (cc2 + cc3*t_cs(i))*t_cs(i)
883 CALL i17abc_pena(i,si_s,r_cs(i),t_cs(i),
884 + bb1 ,bb2 ,bb3 ,cc1 ,cc2 ,cc3 )
885 s_cm(i) = cc1 + (cc2 + cc3*t_cs(i))*t_cs(i)
886 IF(abs(s_cm(i)) < one)THEN
887 icont(i,1) = 1
888 area(i) = min(area_el(i) * aa,area_tot(i))
889 ENDIF
890 ENDIF
891 ENDDO
892c-----------------------------------------------------------------------
893c borne surface totale la surface du main
894c-----------------------------------------------------------------------
895 DO i=1,llt
896 IF(s_cm(i) < zero)THEN
897c main 1234
898 xa = xx(i,11)-xx(i,9)
899 ya = yy(i,11)-yy(i,9)
900 za = zz(i,11)-zz(i,9)
901 xb = xx(i,12)-xx(i,10)
902 yb = yy(i,12)-yy(i,10)
903 zb = zz(i,12)-zz(i,10)
904 ELSE
905c main 5678
906 xa = xx(i,15)-xx(i,13)
907 ya = yy(i,15)-yy(i,13)
908 za = zz(i,15)-zz(i,13)
909 xb = xx(i,16)-xx(i,14)
910 yb = yy(i,16)-yy(i,14)
911 zb = zz(i,16)-zz(i,14)
912 ENDIF
913 xc = ya*zb - za*yb
914 yc = za*xb - xa*zb
915 zc = xa*yb - ya*xb
916 aaa = pi*fourth*sqrt(xc*xc+yc*yc+zc*zc)
917 area_el(i) = min(area_el(i),aaa)
918 area_tot(i) = min(area_tot(i),aaa)
919 IF(area_tot(i) == zero)area_tot(i) = area_el(i)
920 ENDDO
921C-----------------------------------------------
922C calcul de Ni(r,s,t); dNi/dr; dNi/dt; normale second(-n)
923C-----------------------------------------------
924 CALL i17norm(llt ,r_cs ,s_cs ,t_cs ,
925 2 nx ,ny ,nz ,xxs ,yys ,zzs )
926c-----------------------------------------------------------------------
927 RETURN
i17racine_pena
subroutine i17racine_pena(llt, a, b, c, r1, r2)
Definition i17for3.F:1098
i17abc_pena
subroutine i17abc_pena(i, f, r, t, b1, b2, b3, c1, c2, c3)
Definition i17for3.F:936
i17surf
subroutine i17surf(llt, r1, r2, r3, r4, t1, t2, t3, t4, area, xx, yy, zz, sm)
Definition i17for3.F:1159
i17norm
subroutine i17norm(llt, rr, ss, tt, nx, ny, nz, xx, yy, zz)
Definition i17lagm.F:1551
i17abc
subroutine i17abc(llt, f, r, t, b1, b2, b3, c1, c2, c3)
Definition i17lagm.F:1459

◆ i17racine_pena()

subroutine i17racine_pena ( integer llt,
a,
b,
c,
r1,
r2 )

Definition at line 1097 of file i17for3.F.

1098C-----------------------------------------------
1099C
1100C calcul des racines r1,r2 bornees a +- 1
1101C
1102C de a x^2 + b x + c = 0
1103C
1104C la routine retourne -1,+1 s'il n'y a pas de racines
1105C-----------------------------------------------
1106C I m p l i c i t T y p e s
1107C-----------------------------------------------
1108#include "implicit_f.inc"
1109C-----------------------------------------------
1110C G l o b a l P a r a m e t e r s
1111C-----------------------------------------------
1112#include "mvsiz_p.inc"
1113C-----------------------------------------------
1114C D u m m y A r g u m e n t s
1115C-----------------------------------------------
1116 INTEGER LLT
1117 my_real
1118 . c(mvsiz),b(mvsiz),a(mvsiz),r1(mvsiz),r2(mvsiz)
1119C-----------------------------------------------
1120C L o c a l V a r i a b l e s
1121C-----------------------------------------------
1122 INTEGER I
1123 my_real
1124 . d
1125C
1126 DO i=1,llt
1127 d = b(i)*b(i) - four*a(i)*c(i)
1128 IF(abs(a(i))<em20)THEN
1129 IF(abs(b(i))<em20)THEN
1130 r1(i) = -one
1131 r2(i) = one
1132 ELSE
1133 r1(i) = -c(i) / b(i)
1134 r2(i) = b(i) / abs(b(i))
1135 ENDIF
1136 ELSEIF(d<zero)THEN
1137 r1(i) = -one
1138 r2(i) = one
1139 ELSE
1140 d = sqrt( d )
1141 r1(i) = (-b(i) - d) / (two*a(i))
1142 r2(i) = (-b(i) + d) / (two*a(i))
1143 ENDIF
1144 ENDDO
1145c-----------------------------------------------------------------------
1146c borne r (ou t) a +-1 sur l'second
1147c-----------------------------------------------------------------------
1148 RETURN

◆ i17surf()

subroutine i17surf ( integer llt,
r1,
r2,
r3,
r4,
t1,
t2,
t3,
t4,
area,
xx,
yy,
zz,
sm )

Definition at line 1155 of file i17for3.F.

1159C-----------------------------------------------
1160C I m p l i c i t T y p e s
1161C-----------------------------------------------
1162#include "implicit_f.inc"
1163C-----------------------------------------------
1164C G l o b a l P a r a m e t e r s
1165C-----------------------------------------------
1166#include "mvsiz_p.inc"
1167C-----------------------------------------------
1168C D u m m y A r g u m e n t s
1169C-----------------------------------------------
1170 INTEGER LLT
1171 my_real
1172 . r1(mvsiz),r2(mvsiz),r3(mvsiz),r4(mvsiz),
1173 . t1(mvsiz),t2(mvsiz),t3(mvsiz),t4(mvsiz),
1174 . area(mvsiz),xx(mvsiz,*) ,yy(mvsiz,*),zz(mvsiz,*),sm(mvsiz)
1175C-----------------------------------------------
1176C L o c a l V a r i a b l e s
1177C-----------------------------------------------
1178 INTEGER I,K
1179 my_real
1180 . u_m_r,u_p_r,u_m_s,u_p_s,u_m_t,u_p_t,
1181 . ums_umt,ums_upt,ups_umt,ups_upt,
1182 . umr_ums,umr_ups,upr_ums,upr_ups,
1183 . umt_umr,umt_upr,upt_umr,upt_upr,
1184 . a,b,r05,s05,t05,r,t,ni(8),
1185 . x1,x2,x3,y1,y2,y3,z1,z2,z3,pis4
1186C-----------------------------------------------------------------------
1187C calcul de la surface
1188C-----------------------------------------------------------------------
1189 DO i=1,llt
1190C
1191C-----------------------------------------------------------------------
1192 r = r1(i)
1193 t = t1(i)
1194C-----------------------------------------------------------------------
1195 r05 = half*r
1196 t05 = half*t
1197C
1198 u_m_r = half - r05
1199 u_p_r = half + r05
1200 u_m_t = half - t05
1201 u_p_t = half + t05
1202C
1203 ni(1) = u_m_t * u_m_r * (-r-t-one)
1204 ni(2) = u_p_t * u_m_r * (-r+t-one)
1205 ni(3) = u_p_t * u_p_r * ( r+t-one)
1206 ni(4) = u_m_t * u_p_r * ( r-t-one)
1207 a = (one-r*r)
1208 ni(6) = a * u_p_t
1209 ni(8) = a * u_m_t
1210 b = (one-t*t)
1211 ni(5) = b * u_m_r
1212 ni(7) = b * u_p_r
1213 x1 = zero
1214 y1 = zero
1215 z1 = zero
1216 x2 = zero
1217 y2 = zero
1218 z2 = zero
1219C------------------------------------
1220 IF(sm(i) == zero)THEN
1221c secondary 1234 ou 5678
1222 DO k=1,8
1223 x1 = x1-ni(k)*xx(i,k)
1224 y1 = y1-ni(k)*yy(i,k)
1225 z1 = z1-ni(k)*zz(i,k)
1226 ENDDO
1227 ELSEIF(sm(i) < zero)THEN
1228c main 1234
1229 DO k=1,4
1230 x1 = x1-ni(k)*xx(i,k)-ni(k+4)*xx(i,k+8)
1231 y1 = y1-ni(k)*yy(i,k)-ni(k+4)*yy(i,k+8)
1232 z1 = z1-ni(k)*zz(i,k)-ni(k+4)*zz(i,k+8)
1233 ENDDO
1234 ELSEIF(sm(i) > zero)THEN
1235c main 5678
1236 DO k=1,4
1237 x1 = x1-ni(k)*xx(i,k+4)-ni(k+4)*xx(i,k+12)
1238 y1 = y1-ni(k)*yy(i,k+4)-ni(k+4)*yy(i,k+12)
1239 z1 = z1-ni(k)*zz(i,k+4)-ni(k+4)*zz(i,k+12)
1240 ENDDO
1241 ENDIF
1242C------------------------------------
1243C-----------------------------------------------------------------------
1244 r = r2(i)
1245 t = t2(i)
1246C-----------------------------------------------------------------------
1247 r05 = half*r
1248 t05 = half*t
1249C
1250 u_m_r = half - r05
1251 u_p_r = half + r05
1252 u_m_t = half - t05
1253 u_p_t = half + t05
1254C
1255 ni(1) = u_m_t * u_m_r * (-r-t-one)
1256 ni(2) = u_p_t * u_m_r * (-r+t-one)
1257 ni(3) = u_p_t * u_p_r * ( r+t-one)
1258 ni(4) = u_m_t * u_p_r * ( r-t-one)
1259 a = (one-r*r)
1260 ni(6) = a * u_p_t
1261 ni(8) = a * u_m_t
1262 b = (one-t*t)
1263 ni(5) = b * u_m_r
1264 ni(7) = b * u_p_r
1265C------------------------------------
1266 IF(sm(i) == zero)THEN
1267c secondary 1234 ou 5678
1268 DO k=1,8
1269 x1 = x1+ni(k)*xx(i,k)
1270 y1 = y1+ni(k)*yy(i,k)
1271 z1 = z1+ni(k)*zz(i,k)
1272 ENDDO
1273 ELSEIF(sm(i) < zero)THEN
1274c main 1234
1275 DO k=1,4
1276 x1 = x1+ni(k)*xx(i,k)+ni(k+4)*xx(i,k+8)
1277 y1 = y1+ni(k)*yy(i,k)+ni(k+4)*yy(i,k+8)
1278 z1 = z1+ni(k)*zz(i,k)+ni(k+4)*zz(i,k+8)
1279 ENDDO
1280 ELSEIF(sm(i) > zero)THEN
1281c main 5678
1282 DO k=1,4
1283 x1 = x1+ni(k)*xx(i,k+4)+ni(k+4)*xx(i,k+12)
1284 y1 = y1+ni(k)*yy(i,k+4)+ni(k+4)*yy(i,k+12)
1285 z1 = z1+ni(k)*zz(i,k+4)+ni(k+4)*zz(i,k+12)
1286 ENDDO
1287 ENDIF
1288C------------------------------------
1289C-----------------------------------------------------------------------
1290 r = r3(i)
1291 t = t3(i)
1292C-----------------------------------------------------------------------
1293 r05 = half*r
1294 t05 = half*t
1295C
1296 u_m_r = half - r05
1297 u_p_r = half + r05
1298 u_m_t = half - t05
1299 u_p_t = half + t05
1300C
1301 ni(1) = u_m_t * u_m_r * (-r-t-one)
1302 ni(2) = u_p_t * u_m_r * (-r+t-one)
1303 ni(3) = u_p_t * u_p_r * ( r+t-one)
1304 ni(4) = u_m_t * u_p_r * ( r-t-one)
1305 a = (one-r*r)
1306 ni(6) = a * u_p_t
1307 ni(8) = a * u_m_t
1308 b = (one-t*t)
1309 ni(5) = b * u_m_r
1310 ni(7) = b * u_p_r
1311C------------------------------------
1312 IF(sm(i) == zero)THEN
1313c secondary 1234 ou 5678
1314 DO k=1,8
1315 x2 = x2-ni(k)*xx(i,k)
1316 y2 = y2-ni(k)*yy(i,k)
1317 z2 = z2-ni(k)*zz(i,k)
1318 ENDDO
1319 ELSEIF(sm(i) < zero)THEN
1320c main 1234
1321 DO k=1,4
1322 x2 = x2-ni(k)*xx(i,k)-ni(k+4)*xx(i,k+8)
1323 y2 = y2-ni(k)*yy(i,k)-ni(k+4)*yy(i,k+8)
1324 z2 = z2-ni(k)*zz(i,k)-ni(k+4)*zz(i,k+8)
1325 ENDDO
1326 ELSEIF(sm(i) > zero)THEN
1327c main 5678
1328 DO k=1,4
1329 x2 = x2-ni(k)*xx(i,k+4)-ni(k+4)*xx(i,k+12)
1330 y2 = y2-ni(k)*yy(i,k+4)-ni(k+4)*yy(i,k+12)
1331 z2 = z2-ni(k)*zz(i,k+4)-ni(k+4)*zz(i,k+12)
1332 ENDDO
1333 ENDIF
1334C------------------------------------
1335C-----------------------------------------------------------------------
1336 r = r4(i)
1337 t = t4(i)
1338C-----------------------------------------------------------------------
1339 r05 = half*r
1340 t05 = half*t
1341C
1342 u_m_r = half - r05
1343 u_p_r = half + r05
1344 u_m_t = half - t05
1345 u_p_t = half + t05
1346C
1347 ni(1) = u_m_t * u_m_r * (-r-t-one)
1348 ni(2) = u_p_t * u_m_r * (-r+t-one)
1349 ni(3) = u_p_t * u_p_r * ( r+t-one)
1350 ni(4) = u_m_t * u_p_r * ( r-t-one)
1351 a = (one-r*r)
1352 ni(6) = a * u_p_t
1353 ni(8) = a * u_m_t
1354 b = (one-t*t)
1355 ni(5) = b * u_m_r
1356 ni(7) = b * u_p_r
1357C------------------------------------
1358 IF(sm(i) == zero)THEN
1359c secondary 1234 ou 5678
1360 DO k=1,8
1361 x2 = x2+ni(k)*xx(i,k)
1362 y2 = y2+ni(k)*yy(i,k)
1363 z2 = z2+ni(k)*zz(i,k)
1364 ENDDO
1365 ELSEIF(sm(i) < zero)THEN
1366c main 1234
1367 DO k=1,4
1368 x2 = x2+ni(k)*xx(i,k)+ni(k+4)*xx(i,k+8)
1369 y2 = y2+ni(k)*yy(i,k)+ni(k+4)*yy(i,k+8)
1370 z2 = z2+ni(k)*zz(i,k)+ni(k+4)*zz(i,k+8)
1371 ENDDO
1372 ELSEIF(sm(i) > zero)THEN
1373c main 5678
1374 DO k=1,4
1375 x2 = x2+ni(k)*xx(i,k+4)+ni(k+4)*xx(i,k+12)
1376 y2 = y2+ni(k)*yy(i,k+4)+ni(k+4)*yy(i,k+12)
1377 z2 = z2+ni(k)*zz(i,k+4)+ni(k+4)*zz(i,k+12)
1378 ENDDO
1379 ENDIF
1380C------------------------------------
1381C Surface de l'ellipse
1382C------------------------------------
1383 x3 = y1*z2 - z1*y2
1384 y3 = z1*x2 - x1*z2
1385 z3 = x1*y2 - y1*x2
1386 area(i) = pi*fourth*sqrt(x3*x3+y3*y3+z3*z3)
1387 ENDDO
1388C-----------------------------------------------
1389 RETURN

◆ i17vit_pena()

subroutine i17vit_pena ( integer llt,
integer nint,
v,
a,
integer, dimension(mvsiz,17) iii,
integer, dimension(mvsiz,16) iiis,
ni_m,
ni_s,
nx,
ny,
nz,
vit,
integer, dimension(mvsiz,4) icont,
rm,
tm,
rs,
ts,
sm,
integer, dimension(mvsiz) les )

Definition at line 1400 of file i17for3.F.

1405C-----------------------------------------------
1406C M o d u l e s
1407C-----------------------------------------------
1408 USE tri7box
1409C-----------------------------------------------
1410C I m p l i c i t T y p e s
1411C-----------------------------------------------
1412#include "implicit_f.inc"
1413C-----------------------------------------------
1414C G l o b a l P a r a m e t e r s
1415C-----------------------------------------------
1416#include "mvsiz_p.inc"
1417C-----------------------------------------------
1418C D u m m y A r g u m e n t s
1419C-----------------------------------------------
1420 INTEGER LLT,NINT
1421 INTEGER III(MVSIZ,17),IIIS(MVSIZ,16),
1422 + ICONT(MVSIZ,4), LES(MVSIZ)
1423C REAL
1424 my_real
1425 . v(3,*),a(3,*),vit(*)
1426 my_real
1427 . rm(mvsiz) ,rs(mvsiz) ,tm(mvsiz) ,ts(mvsiz) ,sm(mvsiz),
1428 . ni_m(mvsiz,*) ,ni_s(mvsiz,*),nx(mvsiz) ,ny(mvsiz) ,nz(mvsiz)
1429C-----------------------------------------------
1430C L o c a l V a r i a b l e s
1431C-----------------------------------------------
1432 INTEGER I,J,IK,NK,I1,I2,I3,I4,IES,IIS
1433 my_real
1434 . vx,vy,vz,vn,aa
1435C-----------------------------------------------------------------------
1436C calcul de Ni sur face 1 2 3 4 ou 5 6 7 8 du main (s=+-1)
1437C-----------------------------------------------------------------------
1438 CALL i17ni(llt,rm ,tm ,ni_m )
1439C-----------------------------------------------------------------------
1440C calcul de Ni sur face 1 2 3 4 ou 5 6 7 8 de l'second (s=+-1)
1441C-----------------------------------------------------------------------
1442 CALL i17ni(llt,rs ,ts ,ni_s )
1443C-----------------------------------------------
1444C calcul de [L]
1445C-----------------------------------------------
1446 DO i=1,llt
1447C-----------------------------------------------
1448C test si contact
1449C-----------------------------------------------
1450 IF(icont(i,1) /= 0)THEN
1451 vx = zero
1452 vy = zero
1453 vz = zero
1454 IF(les(i)>0)THEN ! test pour SPMD
1455 DO ik=1,8
1456 vx = vx - (v(1,iii(i,ik)))*ni_m(i,ik)
1457 vy = vy - (v(2,iii(i,ik)))*ni_m(i,ik)
1458 vz = vz - (v(3,iii(i,ik)))*ni_m(i,ik)
1459 vx = vx + (v(1,iiis(i,ik)))*ni_s(i,ik)
1460 vy = vy + (v(2,iiis(i,ik)))*ni_s(i,ik)
1461 vz = vz + (v(3,iiis(i,ik)))*ni_s(i,ik)
1462 ENDDO
1463 ELSE ! partie complementaire SPMD
1464 ies = abs(les(i)) ! no local
1465 DO ik=1,8
1466 vx = vx - v(1,iii(i,ik))*ni_m(i,ik)
1467 vy = vy - v(2,iii(i,ik))*ni_m(i,ik)
1468 vz = vz - v(3,iii(i,ik))*ni_m(i,ik)
1469 iis = iiis(i,ik)
1470 vx = vx + vfi17(nint)%P(1,iis,ies)*ni_s(i,ik)
1471 vy = vy + vfi17(nint)%P(2,iis,ies)*ni_s(i,ik)
1472 vz = vz + vfi17(nint)%P(3,iis,ies)*ni_s(i,ik)
1473 ENDDO
1474 END IF
1475C
1476 vn = nx(i)*vx + ny(i)*vy + nz(i)*vz
1477 vit(i)=sm(i)*vn
1478 ENDIF
1479 ENDDO
1480C
1481 RETURN
i17ni
subroutine i17ni(llt, r, t, ni)
Definition i17lagm.F:1252