44
45
46
48
49
50
51#include "implicit_f.inc"
52
53
54
55#include "mvsiz_p.inc"
56
57
58
59#include "com01_c.inc"
60#include "com08_c.inc"
61#include "impl1_c.inc"
62
63
64
65 INTEGER JLT, ITIED,NIN,LREM
66 INTEGER IX1(MVSIZ), IX2(), IX3(MVSIZ), IX4(MVSIZ),
67 . NSVG(MVSIZ), INDEX(MVSIZ)
69 . a(3,*), ms(*), v(3,*),x1(*),x2(*),x3(*),x4(*),
70 . y1(*),y2(*),y3(*),y4(*),z1(*),z2(*),z3(*),z4(*),
71 . cand_f(6,*),fric,off(*),scalk,
72 . vxi(mvsiz),vyi(mvsiz),vzi(mvsiz),msi(mvsiz)
74 . nx1(mvsiz), nx2(mvsiz), nx3(mvsiz), nx4(mvsiz),
75 . ny1(mvsiz), ny2(mvsiz), ny3(mvsiz), ny4(mvsiz),
76 . nz1(mvsiz), nz2(mvsiz), nz3(mvsiz), nz4(mvsiz),
77 . lb1(mvsiz), lb2(mvsiz), lb3(mvsiz), lb4(mvsiz),
78 . lc1(mvsiz), lc2(mvsiz), lc3(mvsiz), lc4(mvsiz),
79 . p1(mvsiz), p2(mvsiz), p3(mvsiz), p4(mvsiz), stif(mvsiz),
80 . gapv(mvsiz),ki11(3,3,mvsiz),kj11(3,3,mvsiz),
81 . kk11(3,3,mvsiz),kl11(3,3,mvsiz),ki12(3,3,mvsiz),
82 . kj12(3,3,mvsiz),kk12(3,3,mvsiz),kl12(3,3,mvsiz)
83
84
85
86 INTEGER I, J1, J, K,IG,ISF,NN,NS,JLTF,NE,II
88 . n1(mvsiz), n2(mvsiz), n3(mvsiz), pene(mvsiz),
89 . h1(mvsiz), h2(mvsiz), h3(mvsiz), h4(mvsiz),
90 . vt1(mvsiz), vt2(mvsiz),fni(mvsiz),
91 . nx(mvsiz), ny(mvsiz), nz(mvsiz),
92 . t1x(mvsiz),t1y(mvsiz),t1z(mvsiz),
93 . t2x(mvsiz),t2y(mvsiz),t2z(mvsiz),norminv,
94 . vx(mvsiz), vy(mvsiz), vz(mvsiz), vn(mvsiz),
95 . s2,fac,facf, h0, la1, la2, la3, la4,fact(mvsiz),
96 . d1,d2,d3,d4,a1,a2,a3,a4,kn(4,mvsiz),q(3,3,mvsiz),fac10
98 . prec,q11,q12,q13,q22,q23,q33,h00,vtx,vty,vtz,vt,
99 . kt1,kt2,kt3,kt4,q1,q2
100
101 fric =one
102 fac10 = ten
103 IF (imp_int7==3) THEN
104 DO i=1,jlt
105 d1 = sqrt(p1(i))
106 p1(i) = fourth*gapv(i)
107 d2 = sqrt(p2(i))
108 p2(i) = fourth*gapv(i)
109 d3 = sqrt(p3(i))
110 p3(i) = fourth*gapv(i)
111 d4 = sqrt(p4(i))
112 p4(i) = fourth*gapv(i)
113 ENDDO
114 ELSE
115 DO i=1,jlt
116
117 d1 = sqrt(p1(i))
118 p1(i) =
max(zero, gapv(i) - d1)
119
120 d2 = sqrt(p2(i))
121 p2(i) =
max(zero, gapv(i) - d2)
122
123 d3 = sqrt(p3(i))
124 p3(i) =
max(zero, gapv(i) - d3)
125
126 d4 = sqrt(p4(i))
127 p4(i) =
max(zero, gapv(i) - d4)
128 ENDDO
129 ENDIF
130
131 DO i=1,jlt
132 IF(ix3(i)/=ix4(i))THEN
133 pene(i) =
max(p1(i),p2(i),p3(i),p4(i))
134
135 la1 = one - lb1(i) - lc1(i)
136 la2 = one - lb2(i) - lc2(i)
137 la3 = one - lb3(i) - lc3(i)
138 la4 = one - lb4(i) - lc4(i)
139
140 h0 = fourth *
141 . (p1(i)*la1 + p2(i)*la2 + p3(i)*la3 + p4(i)*la4)
142 h1(i) = h0 + p1(i) * lb1(i) + p4(i) * lc4(i)
143 h2(i) = h0 + p2(i) * lb2(i) + p1(i) * lc1(i)
144 h3(i) = h0 + p3(i) * lb3(i) + p2(i) * lc2(i)
145 h4(i) = h0 + p4(i) * lb4(i) + p3(i) * lc3(i)
146 h00 = one/
max(em20,h1(i) + h2(i) + h3(i) + h4(i))
147 h1(i) = h1(i) * h00
148 h2(i) = h2(i) * h00
149 h3(i) = h3(i) * h00
150 h4(i) = h4(i) * h00
151
152 ELSE
153 pene(i) = p1(i)
154 n1(i) = nx1(i)
155 n2(i) = ny1(i)
156 n3(i) = nz1(i)
157 h1(i) = lb1(i)
158 h2(i) = lc1(i)
159 h3(i) = one - lb1(i) - lc1(i)
160 h4(i) = zero
161 ENDIF
162 ENDDO
163
164
165
166
167
168
169
170
171 DO i=1,jlt
172
173 IF(ix3(i)/=ix4(i))THEN
174 h0 = -fourth*(h1(i) - h2(i) + h3(i) - h4(i))
175 h0 =
min(h0,h2(i),h4(i))
176 h0 =
max(h0,-h1(i),-h3(i))
177 h1(i) = h1(i) + h0
178 h2(i) = h2(i) - h0
179 h3(i) = h3(i) + h0
180 h4(i) = h4(i) - h0
181 ENDIF
182 ENDDO
183
184 DO i=1,jlt
185 ii = index(i)
186 IF(cand_f(1,ii)==zero)THEN
187
188
189
190
191
192
193 ELSE
194
195
196
197 h1(i) = cand_f(4,ii)
198 h2(i) = cand_f(5,ii)
199 h3(i) = cand_f(6,ii)
200 h4(i) = one - h1(i) - h2(i) - h3(i)
201 ENDIF
202 ENDDO
203
204 DO i=1,jlt
205 vx(i) = vxi(i) - h1(i)*v(1,ix1(i)) - h2(i)*v(1,ix2(i))
206 . - h3(i)*v(1,ix3(i)) - h4(i)*v(1,ix4(i))
207 vy(i) = vyi(i) - h1(i)*v(2,ix1(i)) - h2(i)*v(2,ix2(i))
208 . - h3(i)*v(2,ix3(i)) - h4(i)*v(2,ix4(i))
209 vz(i) = vzi(i) - h1(i)*v(3,ix1(i)) - h2(i)*v(3,ix2(i))
210 . - h3(i)*v(3,ix3(i)) - h4(i)*v(3,ix4(i))
211 ENDDO
212
213 DO i=1,jlt
214 t1x(i) = x3(i) - x1(i)
215 t1y(i) = y3(i) - y1(i)
216 t1z(i) = z3(i) - z1(i)
217 norminv = one/sqrt(t1x(i)**2+t1y(i)**2+t1z(i)**2)
218 t1x(i) = t1x(i)*norminv
219 t1y(i) = t1y(i)*norminv
220 t1z(i) = t1z(i)*norminv
221
222 t2x(i) = x4(i) - x2(i)
223 t2y(i) = y4(i) - y2(i)
224 t2z(i) = z4(i) - z2(i)
225
226 nx(i) = t1y(i)*t2z(i) - t1z(i)*t2y(i)
227 ny(i) = t1z(i)*t2x(i) - t1x(i)*t2z(i)
228 nz(i) = t1x(i)*t2y(i) - t1y(i)*t2x(i)
229 norminv = one/sqrt(nx(i)**2+ny(i)**2+nz(i)**2)
230 nx(i) = nx(i)*norminv
231 ny(i) = ny(i)*norminv
232 nz(i) = nz(i)*norminv
233
234 t2x(i) = ny(i)*t1z(i) - nz(i)*t1y(i)
235 t2y(i) = nz(i)*t1x(i) - nx(i)*t1z(i)
236 t2z(i) = nx(i)*t1y(i) - ny(i)*t1x(i)
237
238 vn(i) = vx(i)*nx(i) + vy(i)*ny(i) + vz(i)*nz(i)
239 vt1(i) = vx(i)*t1x(i) + vy(i)*t1y(i) + vz(i)*t1z(i)
240 vt2(i) = vx(i)*t2x(i) + vy(i)*t2y(i) + vz(i)*t2z(i)
241 ENDDO
242
243 DO i=1,jlt
244 IF(pene(i)==zero.AND.cand_f(1,index(i))==zero)THEN
245
246
247
248 vn(i) = zero
249 vt1(i) = zero
250 vt2(i) = zero
251 ENDIF
252 ENDDO
253
254 DO i=1,jlt
255 ii = index(i)
256 fni(i) = cand_f(1,ii) + vn(i) * dt1 * stif(i)
257 ENDDO
258
259 DO 100 i=1,jlt
260 IF(itied==0)THEN
261 IF(cand_f(1,index(i))*fni(i)<zero)THEN
262
263
264
265 fni(i) = zero
266 vn(i) = zero
267 vt1(i) = zero
268 vt2(i) = zero
269 stif(i) = zero
270 ELSE
271
272 ENDIF
273 ELSE
274 stif(i) = stif(i) * abs(vn(i)) * dt1/
max(pene(i),em10)
275 ENDIF
276
277 100 CONTINUE
278
279
280
281 DO i=1,jlt
282 IF (abs(vt1(i))>zero.OR.abs(vt2(i))>zero) THEN
283 q(1,1,i)=t1x(i)
284 q(1,2,i)=t1y(i)
285 q(1,3,i)=t1z(i)
286 q(3,1,i)=nx(i)
287 q(3,2,i)=ny(i)
288 q(3,3,i)=nz(i)
289 q(2,1,i)=t2x(i)
290 q(2,2,i)=t2y(i)
291 q(2,3,i)=t2z(i)
292 fact(i)=fric
293 ELSE
294 fact(i)=zero
295 ENDIF
296 ENDDO
297 IF (scalk<0) THEN
298 isf=1
299 ELSE
300 isf=0
301 ENDIF
302 facf=fac10*abs(scalk)
303 IF (isf==1) THEN
304 DO i=1,jlt
305 IF (vn(i)>zero) THEN
306 fac=stif(i)*facf
307 ELSEIF (vn(i)<zero) THEN
308 fac=stif(i)/facf
309 ELSE
310 fac=stif(i)
311 ENDIF
312 kn(1,i)=fac*h1(i)
313 kn(2,i)=fac*h2(i)
314 kn(3,i)=fac*h3(i)
315 kn(4,i)=fac*h4(i)
316 fact(i)=fac*fact(i)
317 ENDDO
318 ELSE
319 DO i=1,jlt
320 fac=stif(i)*facf
321 kn(1,i)=fac*h1(i)
322 kn(2,i)=fac*h2(i)
323 kn(3,i)=fac*h3(i)
324 kn(4,i)=fac*h4(i)
325 fact(i)=fac*fact(i)
326 ENDDO
327 ENDIF
328 DO i=1,jlt
329 q11=nx(i)*nx(i)
330 q12=nx(i)*ny(i)
331 q13=nx(i)*nz(i)
332 q22=ny(i)*ny(i)
333 q23=ny(i)*nz(i)
334 q33=nz(i)*nz(i)
335 ki11(1,1,i)=kn(1,i)*q11
336 ki11(1,2,i)=kn(1,i)*q12
337 ki11(1,3,i)=kn(1,i)*q13
338 ki11(2,2,i)=kn(1,i)*q22
339 ki11(2,3,i)=kn(1,i)*q23
340 ki11(3,3,i)=kn(1,i)*q33
341 kj11(1,1,i)=kn(2,i)*q11
342 kj11(1,2,i)=kn(2,i)*q12
343 kj11(1,3,i)=kn(2,i)*q13
344 kj11(2,2,i)=kn(2,i)*q22
345 kj11(2,3,i)=kn(2,i)*q23
346 kj11(3,3,i)=kn(2,i)*q33
347 kk11(1,1,i)=kn(3,i)*q11
348 kk11(1,2,i)=kn(3,i)*q12
349 kk11(1,3,i)=kn(3,i)*q13
350 kk11(2,2,i)=kn(3,i)*q22
351 kk11(2,3,i)=kn(3,i)*q23
352 kk11(3,3,i)=kn(3,i)*q33
353 kl11(1,1,i)=kn(4,i)*q11
354 kl11(1,2,i)=kn(4,i)*q12
355 kl11(1,3,i)=kn(4,i)*q13
356 kl11(2,2,i)=kn(4,i)*q22
357 kl11(2,3,i)=kn(4,i)*q23
358 kl11(3,3,i)=kn(4,i)*q33
359 ENDDO
360
361 DO j=1,3
362 DO k=j,3
363 DO i=1,jlt
364 IF (fact(i)>zero) THEN
365 q1 =q(1,j,i)*q(1,k,i)
366 q2 =q(2,j,i)*q(2,k,i)
367 fac=fact(i)*(q1+q2)
368 kt1=fac*h1(i)
369 ki11(j,k,i)=ki11(j,k,i)+kt1
370 kt2=fac*h2(i)
371 kj11(j,k,i)=kj11(j,k,i)+kt2
372 kt3=fac*h3(i)
373 kk11(j,k,i)=kk11(j,k,i)+kt3
374 kt4=fac*h4(i)
375 kl11(j,k,i)=kl11(j,k,i)+kt4
376 ENDIF
377 ENDDO
378 ENDDO
379 ENDDO
380
381 DO j=1,3
382 DO k=j,3
383 DO i=1,jlt
384 ki12(j,k,i)=-ki11(j,k,i)
385 kj12(j,k,i)=-kj11(j,k,i)
386 kk12(j,k,i)=-kk11(j,k,i)
387 kl12(j,k,i)=-kl11(j,k,i)
388 ENDDO
389 ENDDO
390 ENDDO
391 DO j=1,3
392 DO k=j+1,3
393 DO i=1,jlt
394 ki12(k,j,i)=-ki11(j,k,i)
395 kj12(k,j,i)=-kj11(j,k,i)
396 kk12(k,j,i)=-kk11(j,k,i)
397 kl12(k,j,i)=-kl11(j,k,i)
398 ENDDO
399 ENDDO
400 ENDDO
401
402 DO i=1,jlt
403 off(i)=one
404 ENDDO
405 IF (nspmd>1) THEN
407 DO i=1,jlt
408 IF(nsvg(i)<0) THEN
409 nn=-nsvg(i)
411
412 ffi(1,ns)=zero
413 ffi(2,ns)=zero
414 ffi(3,ns)=zero
415 dfi(1,ns)=zero
416 dfi(2,ns)=zero
417 dfi(3,ns)=zero
418 ENDIF
419 ENDDO
420 ELSE
421 jltf = 0
422 DO i=1,jlt
423 IF(nsvg(i)<0) THEN
424 nn=-nsvg(i)
425 jltf = jltf + 1
428 stifs(ne)=stif(i)
429 h_e(1,ne)=h1(i)
430 h_e(2,ne)=h2(i)
431 h_e(3,ne)=h3(i)
432 h_e(4,ne)=h4(i)
433 n_e(1,ne)=nx(i)
434 n_e(2,ne)=ny(i)
435 n_e(3,ne)=nz(i)
436
437 ffi(1,ns)=zero
438 ffi(2,ns)=zero
439 ffi(3,ns)=zero
440 dfi(1,ns)=zero
441 dfi(2,ns)=zero
442 dfi(3,ns)=zero
443 ENDIF
444 ENDDO
445 ENDIF
446 ENDIF
447
448 RETURN
integer, dimension(:), allocatable shf_int
type(int_pointer2), dimension(:), allocatable ind_int