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"

Go to the source code of this file.

Functions/Subroutines

subroutine i17for3 (output, 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 (output, 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 (output, 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 (output, 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 941 of file i17for3.F.

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

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

◆ i17for3()

subroutine i17for3 ( type(output_), intent(inout) output,
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 42 of file i17for3.F.

50C-----------------------------------------------
51C M o d u l e s
52C-----------------------------------------------
53 USE tri7box
54 USE icontact_mod
55 USE message_mod
56 USE h3d_mod
57 USE output_mod, ONLY : output_
58 use element_mod , only : nixs
59C-----------------------------------------------
60C I m p l i c i t T y p e s
61C-----------------------------------------------
62#include "implicit_f.inc"
63#include "comlock.inc"
64C-----------------------------------------------
65C G l o b a l P a r a m e t e r s
66C-----------------------------------------------
67#include "mvsiz_p.inc"
68C-----------------------------------------------
69C C o m m o n B l o c k s
70C-----------------------------------------------
71#include "task_c.inc"
72#include "com01_c.inc"
73#include "com04_c.inc"
74#include "parit_c.inc"
75 COMMON /i17globi/ie_min,ies_min
76 COMMON /i17globr/vit_min
77C-----------------------------------------------
78C D u m m y A r g u m e n t s
79C-----------------------------------------------
80 TYPE(OUTPUT_), INTENT(INOUT) :: OUTPUT
81 INTEGER I_STOK,ITASK,ITIED,NINT,NME,NSE,NISKYFI,LSKYI17,NOINT,
82 . CANDN(*),CANDE(*), ISKY(*),
83 . IXS(NIXS,*),IXS16(8,*),NELES(*) ,NELEM(*)
84C REAL
85 my_real
86 . x(3,*),v(3,*),eminxm(6,*),eminxs(6,*),a(3,*),stifn(*),
87 . fskyi ,frotm(7,*),frots(7,*),km(*),ks(*),fric(*), fsav(*),
88 . fcont(3,*), ms(*)
89 TYPE(H3D_DATABASE) :: H3D_DATA
90C-----------------------------------------------
91C L o c a l V a r i a b l e s
92C-----------------------------------------------
93 INTEGER I,J,K,IK,IE,IS,IC,NK,III(MVSIZ,17),LLT,NFT,LE,FIRST,LAST,
94 . I16,LES,IES,IIIS(MVSIZ,16),LEL(MVSIZ),LESL(MVSIZ),
95 . IE_MIN,IES_MIN, IERR1, IERR2
96 my_real
97 . xx(mvsiz,17),yy(mvsiz,17),zz(mvsiz,17),
98 . xxs(mvsiz,16),yys(mvsiz,16),zzs(mvsiz,16),
99 . vit_min
100C-----------------------------------------------------------------------
101C loop over contact candidates
102C-----------------------------------------------------------------------
103 vit_min = zero
104 ie_min = 99999999
105 ies_min = 99999999
106 nskyi17 = 0
107 nrskyi17 = 0
108 IF(iparit/=0.AND.itask==0)THEN
109 ierr1=0
110 ierr2=0
111 ALLOCATE(iskyi17(lskyi17),stat=ierr2)
112 ierr1 = ierr1+ierr2
113 ALLOCATE(fskyi17(lskyi17,4),stat=ierr2)
114 ierr1 = ierr1+ierr2
115 IF(nspmd>1)THEN
116 ALLOCATE(irskyi17(lskyi17),stat=ierr2)
117 ierr1 = ierr1+ierr2
118 ALLOCATE(frskyi17(4,lskyi17),stat=ierr2)
119 ierr1 = ierr1+ierr2
120 END IF
121 IF(ierr1/=0)THEN
122 CALL ancmsg(msgid=88,anmode=aninfo,
123 . i1=noint)
124 CALL arret(2)
125 END IF
126 END IF
127C-----------------------------------------------------------------------
128 CALL my_barrier
129C-----------------------------------------------------------------------
130 first = 1 + i_stok * itask / nthread
131 last = i_stok*(itask+1) / nthread
132 llt = 0
133 nft=llt+1
134 DO ic=first,last
135 le = cande(ic)
136 les = candn(ic)
137 ie = nelem(le)
138 IF(les<=nse)THEN ! candidat interne en SPMD
139 ies = neles(les)
140C-----------------------------------------------------------------------
141C test if inside the box
142C-----------------------------------------------------------------------
143 IF(le >0.AND.les>0.AND.
144 . eminxs(4,les)>eminxm(1,le).AND.
145 . eminxs(5,les)>eminxm(2,le).AND.
146 . eminxs(6,les)>eminxm(3,le).AND.
147 . eminxs(1,les)<eminxm(4,le).AND.
148 . eminxs(2,les)<eminxm(5,le).AND.
149 . eminxs(3,les)<eminxm(6,le))THEN
150 llt = llt+1
151 lel(llt)=le
152 lesl(llt)=les
153 DO k=1,4
154C III 1,2,3,4,5,6,7,8, 9,10,11,12,13,14,15,16
155 iii(llt,k) =ixs(k+1,ie)
156 iii(llt,k+4) =ixs(k+5,ie)
157 iii(llt,k+8) =ixs16(k,ie-numels8-numels10-numels20)
158 iii(llt,k+12)=ixs16(k+4,ie-numels8-numels10-numels20)
159C IIIS 1,2,3,4,9,10,11,12, 5,6,7,8,13,14,15,16
160 iiis(llt,k) =ixs(k+1,ies)
161 iiis(llt,k+8) =ixs(k+5,ies)
162 iiis(llt,k+4)=ixs16(k,ies-numels8-numels10-numels20)
163 iiis(llt,k+12)=ixs16(k+4,ies-numels8-numels10-numels20)
164 ENDDO
165 DO k=1,16
166 i = iii(llt,k)
167 xx(llt,k)=x(1,i)
168 yy(llt,k)=x(2,i)
169 zz(llt,k)=x(3,i)
170 i = iiis(llt,k)
171 xxs(llt,k)=x(1,i)
172 yys(llt,k)=x(2,i)
173 zzs(llt,k)=x(3,i)
174 ENDDO
175 END IF
176 ELSE ! Complement SPMD front
177 ies = les-nse
178C-----------------------------------------------------------------------
179C test if inside the box
180C-----------------------------------------------------------------------
181 IF(le >0.AND.les>0.AND.
182 . eminxfi(nint)%P(4,ies)>eminxm(1,le).AND.
183 . eminxfi(nint)%P(5,ies)>eminxm(2,le).AND.
184 . eminxfi(nint)%P(6,ies)>eminxm(3,le).AND.
185 . eminxfi(nint)%P(1,ies)<eminxm(4,le).AND.
186 . eminxfi(nint)%P(2,ies)<eminxm(5,le).AND.
187 . eminxfi(nint)%P(3,ies)<eminxm(6,le))THEN
188 llt = llt+1
189 lel(llt)=le
190 lesl(llt)=-ies ! reperage en negatif des candidats remote SPMD
191#include "mic_lockon.inc"
192 nsvfi(nint)%P(ies) = -abs(nsvfi(nint)%P(ies))
193#include "mic_lockoff.inc"
194 DO k=1,4
195C III 1,2,3,4,5,6,7,8, 9,10,11,12,13,14,15,16
196 iii(llt,k) =ixs(k+1,ie)
197 iii(llt,k+4) =ixs(k+5,ie)
198 iii(llt,k+8) =ixs16(k,ie-numels8-numels10-numels20)
199 iii(llt,k+12)=ixs16(k+4,ie-numels8-numels10-numels20)
200C IIIS 1,2,3,4,9,10,11,12, 5,6,7,8,13,14,15,16
201 iiis(llt,k) =k
202 iiis(llt,k+8) =k+4
203 iiis(llt,k+4) =k+8
204 iiis(llt,k+12)=k+12
205 ENDDO
206 DO k=1,16
207 i = iii(llt,k)
208 xx(llt,k)=x(1,i)
209 yy(llt,k)=x(2,i)
210 zz(llt,k)=x(3,i)
211 i = iiis(llt,k) ! I indicates the node number from 1 to 16
212 xxs(llt,k)=xfi17(nint)%P(1,i,ies)
213 yys(llt,k)=xfi17(nint)%P(2,i,ies)
214 zzs(llt,k)=xfi17(nint)%P(3,i,ies)
215 ENDDO
216 END IF
217 END IF
218C-----------------------------------------------------------------------
219C calculation of F by mvsiz packet
220C-----------------------------------------------------------------------
221 IF(llt==mvsiz-1)THEN
222 CALL i17lll_pena(output,
223 1 llt ,v ,stifn ,xx ,fric ,
224 2 yy ,zz ,iii ,fskyi ,isky ,
225 3 a ,x ,itied ,nint ,
226 4 xxs ,yys ,zzs ,iiis ,vit_min ,
227 5 lel ,lesl ,ie_min ,ies_min ,itask ,
228 6 frotm ,frots ,km ,ks ,fsav ,
229 7 fcont ,ms ,niskyfi ,noint ,h3d_data )
230 nft=llt+1
231 llt = 0
232 ENDIF
233 ENDDO
234C-----------------------------------------------------------------------
235C calculation of F for last packet
236C-----------------------------------------------------------------------
237 IF(llt/=0) CALL i17lll_pena(output,
238 1 llt ,v ,stifn ,xx ,fric ,
239 2 yy ,zz ,iii ,fskyi ,isky ,
240 3 a ,x ,itied ,nint ,
241 4 xxs ,yys ,zzs ,iiis ,vit_min ,
242 5 lel ,lesl ,ie_min ,ies_min ,itask ,
243 6 frotm ,frots ,km ,ks ,fsav ,
244 7 fcont ,ms ,niskyfi ,noint ,h3d_data )
245C-----------------------------------------------------------------------
246C parallel Arithmetic for assembly (frots)
247C-----------------------------------------------------------------------
248 IF(iparit/=0) THEN
249 CALL my_barrier
250 IF(itask==0)THEN
251 IF(nspmd>1)THEN
252 CALL spmd_i17frots_pon(
253 1 nskyi17 ,iskyi17,fskyi17,nrskyi17,irskyi17,
254 2 frskyi17,nint ,lskyi17,noint)
255 DEALLOCATE(irskyi17)
256 DEALLOCATE(frskyi17)
257 END IF
258C
259 IF(nskyi17>0)THEN
260 CALL i17lll_ass(nskyi17 ,iskyi17 ,fskyi17, frots, nse,
261 2 lskyi17 )
262 END IF
263 DEALLOCATE(iskyi17)
264 DEALLOCATE(fskyi17)
265 END IF
266 END IF
267C
268 RETURN
i17lll_pena
subroutine i17lll_pena(output, 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:291
i17lll_ass
subroutine i17lll_ass(nskyi17, iskyi17, fskyi17, frots, nmes, lskyi17)
Definition i17for3.F:2062
h3d_mod
Definition h3d_mod.F:311
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:1257
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:895
arret
subroutine arret(nn)
Definition arret.F:86
my_barrier
subroutine my_barrier
Definition machine.F:31

◆ i17lll4_pena()

subroutine i17lll4_pena ( type(output_), intent(inout) output,
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 1505 of file i17for3.F.

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

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

◆ i17lll_pena()

subroutine i17lll_pena ( type(output_), intent(inout) output,
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 283 of file i17for3.F.

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

◆ i17mini_pena()

subroutine i17mini_pena ( type(output_), intent(inout) output,
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 501 of file i17for3.F.

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

◆ i17racine_pena()

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

Definition at line 1104 of file i17for3.F.

1105C-----------------------------------------------
1106C
1107C calculation des racines r1,r2 bornees a +- 1
1108C
1109C de a x^2 + b x + c = 0
1110C
1111C The routine returns -1,+1 if there are no roots
1112C-----------------------------------------------
1113C I m p l i c i t T y p e s
1114C-----------------------------------------------
1115#include "implicit_f.inc"
1116C-----------------------------------------------
1117C G l o b a l P a r a m e t e r s
1118C-----------------------------------------------
1119#include "mvsiz_p.inc"
1120C-----------------------------------------------
1121C D u m m y A r g u m e n t s
1122C-----------------------------------------------
1123 INTEGER LLT
1124 my_real
1125 . c(mvsiz),b(mvsiz),a(mvsiz),r1(mvsiz),r2(mvsiz)
1126C-----------------------------------------------
1127C L o c a l V a r i a b l e s
1128C-----------------------------------------------
1129 INTEGER I
1130 my_real
1131 . d
1132C
1133 DO i=1,llt
1134 d = b(i)*b(i) - four*a(i)*c(i)
1135 IF(abs(a(i))<em20)THEN
1136 IF(abs(b(i))<em20)THEN
1137 r1(i) = -one
1138 r2(i) = one
1139 ELSE
1140 r1(i) = -c(i) / b(i)
1141 r2(i) = b(i) / abs(b(i))
1142 ENDIF
1143 ELSEIF(d<zero)THEN
1144 r1(i) = -one
1145 r2(i) = one
1146 ELSE
1147 d = sqrt( d )
1148 r1(i) = (-b(i) - d) / (two*a(i))
1149 r2(i) = (-b(i) + d) / (two*a(i))
1150 ENDIF
1151 ENDDO
1152c-----------------------------------------------------------------------
1153c bound r (or t) at +-1 on the second
1154c-----------------------------------------------------------------------
1155 RETURN

◆ i17surf()

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

Definition at line 1162 of file i17for3.F.

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

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