45
46
47
48#include "implicit_f.inc"
49
50
51
52#include "mvsiz_p.inc"
53
54
55
56 INTEGER, INTENT(INOUT) :: LFT
57 INTEGER, INTENT(INOUT) :: LLT
58 INTEGER, INTENT(INOUT) :: NFT
59 INTEGER IRECT(4,*),MSR(*),NSV(*),IRTL(*),IRTLO(*),NPC(*),LOLD(*)
60 INTEGER ,IFRICV
61
63 . x(3,*),cst(2,*),fric0(3,*),es(*),em(*),tf(*),vnt(*)
65 . fric,sfric,ascalf,ascalv,stiff
66 my_real,
DIMENSION(MVSIZ),
INTENT(IN) :: n1,n2,n3,ssc,ttc,xface
67 my_real,
DIMENSION(MVSIZ),
INTENT(IN) :: xp,yp,zp
68 my_real,
DIMENSION(MVSIZ),
INTENT(IN) :: h1,h2,h3,h4,fni
69 my_real,
DIMENSION(MVSIZ),
INTENT(INOUT) :: fxi,fyi,fzi
70 my_real,
DIMENSION(MVSIZ),
INTENT(INOUT) :: fx1,fx2,fx3,fx4
71 my_real,
DIMENSION(MVSIZ),
INTENT(INOUT) :: fy1,fy2,fy3,fy4
72 my_real,
DIMENSION(MVSIZ),
INTENT(INOUT) :: fz1,fz2,fz3,fz4
73
74
75
76
77
78
79 INTEGER I, IL, JJ, NN, L, J3, J2, J1, I3, I2, I1
81 . h(4), xx1(4), xx2(4), xx3(4)
83 . ss0,tt0,xc0,yc0,zc0,sp,sm,tp,tm,ft,ftx,fty,ftz,
84 . ansx,ansy,ansz,fmax,fti,fn,tn1,tn2,tn3,tn
86 . kfric(llt),kfricv(llt),xx(llt)
87
88
89
90 IF (ifricf == 0) THEN
91 DO i=lft,llt
92 kfric(i) = fric
93 ENDDO
94 ELSE
95 DO i=lft,llt
96 xx(i) = fni(i)*ascalf
97 ENDDO
98 CALL ninterp(ifricf,npc,tf,llt,xx,kfric)
99 DO i=lft,llt
100 kfric(i) = fric*kfric(i)
101 ENDDO
102 ENDIF
103 IF (ifricv /= 0) THEN
104 DO i=lft,llt
105 xx(i) = vnt(i)*ascalv
106 ENDDO
107 CALL ninterp(ifricv,npc,tf,llt,xx,kfricv)
108 kfric(i) = kfric(i)*kfricv(i)
109 ENDIF
110
111 DO i=lft,llt
112 il = i+nft
113 IF (lold(i) == 0) THEN
114 fxi(i) = zero
115 fyi(i) = zero
116 fzi(i) = zero
117 fric0(1,il) = zero
118 fric0(2,il) = zero
119 fric0(3,il) = zero
120 ELSE
121
122
123
124 ss0 = cst(1,il)
125 tt0 = cst(2,il)
126
127 DO jj=1,4
128 nn=msr(irect(jj,iabs(irtlo(il))))
129 xx1(jj) = x(1,nn)
130 xx2(jj) = x(2,nn)
131 xx3(jj) = x(3,nn)
132 ENDDO
133
134 sp = one + ss0
135 sm = one - ss0
136 tp = fourth*(one + tt0)
137 tm = fourth*(one - tt0)
138 h(1)=tm*sm
139 h(2)=tm*sp
140 h(3)=tp*sp
141 h(4)=tp*sm
142
143 xc0 = zero
144 yc0 = zero
145 zc0 = zero
146 DO jj=1,4
147 xc0 = xc0 + h(jj)*xx1(jj)
148 yc0 = yc0 + h(jj)*xx2(jj)
149 zc0 = zc0 + h(jj)*xx3(jj)
150 ENDDO
151
152 ansx = xp(i) - xc0
153 ansy = yp(i) - yc0
154 ansz = zp(i) - zc0
155
156 fmax = sfric -
min(kfric(i)*fni(i),zero)
157
158 fxi(i) = fric0(1,il) + ansx*stiff
159 fyi(i) = fric0(2,il) + ansy*stiff
160 fzi(i) = fric0(3,il) + ansz*stiff
161
162 fn = fxi(i)*n1(i)+fyi(i)*n2(i)+fzi(i)*n3(i)
163 ftx = fxi(i) - n1(i)*fn
164 fty = fyi(i) - n2(i)*fn
165 ftz = fzi(i) - n3(i)*fn
166 ft = sqrt(ftx*ftx + fty*fty + ftz*ftz)
167 IF (ft /= zero) THEN
168 tn1 = ftx/ft
169 tn2 = fty/ft
170 tn3 = ftz/ft
171 ELSE
172 tn3 = zero
173 tn = sqrt(n1(i)*n1(i)+n2(i)*n2(i))
174 IF(tn/=zero)THEN
175 tn2 =-n1(i)/tn
176 tn1 = n2(i)/tn
177 ELSE
178 tn2
179 tn1 = one
180 ENDIF
181 ENDIF
182
183 IF (ft > fmax) THEN
184
185
186
187 fxi(i) = tn1*fmax
188 fyi(i) = tn2*fmax
189 fzi(i) = tn3*fmax
190 irtlo(il) = irtl(il)*xface(i)
191 cst(1,il) = ssc(i)
192 cst(2,il) = ttc(i)
193 ELSE
194
195
196
197 fxi(i) = ftx
198 fyi(i) = fty
199 fzi(i) = ftz
200 ENDIF
201 fric0(1,il) = fxi(i)
202 fric0(2,il) = fyi(i)
203 fric0(3,il) = fzi(i)
204
205 ENDIF
206 ENDDO
207
208
209
210 DO i=lft,llt
211 il=i+nft
212 l =irtl(il)
213 fx1(i)=fxi(i)*h1(i)
214 fy1(i)=fyi(i)*h1(i)
215 fz1(i)=fzi(i)*h1(i)
216
217 fx2(i)=fxi(i)*h2(i)
218 fy2(i)=fyi(i)*h2(i)
219 fz2(i)=fzi(i)*h2(i)
220
221 fx3(i)=fxi(i)*h3(i)
222 fy3(i)=fyi(i)*h3(i)
223 fz3(i)=fzi(i)*h3(i)
224
225 fx4(i)=fxi(i)*h4(i)
226 fy4(i)=fyi(i)*h4(i)
227 fz4(i)=fzi(i)*h4(i)
228
229
230
231 j3=3*irect(1,l)
232 j2=j3-1
233 j1=j2-1
234 em(j1)=em(j1)+fx1
235 em(j2)=em(j2)+fy1(i)
236 em(j3)=em(j3)+fz1(i)
237
238 j3=3*irect(2,l)
239 j2=j3-1
240 j1=j2-1
241 em(j1)=em(j1)+fx2(i)
242 em(j2)=em(j2)+fy2(i)
243 em(j3)=em(j3)+fz2(i)
244
245 j3=3*irect(3,l)
246 j2=j3-1
247 j1=j2-1
248 em(j1)=em(j1)+fx3(i)
249 em(j2)=em(j2)+fy3(i)
250 em(j3)=em(j3)+fz3
251
252 j3=3*irect(4,l)
253 j2=j3-1
254 j1=j2-1
255 em(j1)=em(j1)+fx4(i)
256 em(j2)=em(j2)+fy4(i)
257 em(j3)=em(j3)+fz4(i)
258
259
260
261 i3=3*il
262 i2=i3-1
263 i1=i2-1
264 es(i1)=es(i1)-fxi(i)
265 es(i2)=es(i2)-fyi(i)
266 es(i3)=es(i3)-fzi(i)
267 ENDDO
268
269
270 RETURN
subroutine ninterp(ifunc, npc, pld, npoint, xv, yv)