43
44
45
48
49
50
51#include "implicit_f.inc"
52#include "comlock.inc"
53
54
55
56#include "mvsiz_p.inc"
57
58
59
60#include "com01_c.inc"
61#include "scr05_c.inc"
62#include "impl1_c.inc"
63
64
65
66 INTEGER JLT,INACTI,NRTS,NIN,LREM,IDESAC
67 INTEGER CS_LOC(MVSIZ), CM_LOC(MVSIZ),N1(MVSIZ), N2(MVSIZ)
69 . stiglo,a(3,*), ms(*), v(3,*),
70 . gapv(*),penis(2,*), penim(2,*),gap, fric,scalk
72 . hs1(mvsiz), hs2(mvsiz), hm1(mvsiz), hm2(mvsiz),
73 . nx(mvsiz), ny(mvsiz), nz(mvsiz), stif(mvsiz),
74 . ms1(mvsiz),ms2(mvsiz),mm1(mvsiz),mm2(mvsiz), off(mvsiz),
75 . vxs1(mvsiz),vys1(mvsiz),vzs1(mvsiz),vxs2(mvsiz),vys2(mvsiz),
76 . vzs2(mvsiz),vxm1(mvsiz),vym1(mvsiz),vzm1(mvsiz),vxm2(mvsiz),
77 . vym2(mvsiz),vzm2(mvsiz),k1i11(3,3,mvsiz),k1j11(3,3,mvsiz),
78 . k2i11(3,3,mvsiz),k2j11(3,3,mvsiz),k1i12(3,3,mvsiz),
79 . k1j12(3,3,mvsiz),k2i12(3,3,mvsiz),k2j12(3,3,mvsiz)
80
81
82
83 INTEGER I, J1, J , K0,K1S,K, NI,ISF,NN,NS,JLTF,NE,NN1,NN2
85 . vx(mvsiz), vy(mvsiz), vz(mvsiz), vn(mvsiz),pene(mvsiz),
86 . vnx, vny, vnz, aa, vmax,s2,dist,rdist,
87 . v2, fm2, dt1inv, visca, fac, ff,
88 . fx, fy, fz, f2, mas2, facm1, pplus
90 . prec,fact(mvsiz),kn(4,mvsiz),q(3,3,mvsiz),
91 . h1(mvsiz), h2(mvsiz), h3(mvsiz), h4(mvsiz)
93 . q11,q12,q13,q22,q23,q33,h00,vtx,vty,vtz,vt,
94 . kt1,kt2,kt3,kt4,q1,q2,facf
95
96 IF (iresp == 1) THEN
97 prec = fiveem4
98 ELSE
99 prec = em10
100 ENDIF
101
102 DO i=1,jlt
103 s2 = sqrt(nx(i)**2 + ny(i)**2 + nz(i)**2)
104 pene(i) = gapv(i) - s2
105 s2 = one/
max(em30,s2)
106 nx(i) = nx(i)*s2
107 ny(i) = ny(i)*s2
108 nz(i) = nz(i)*s2
109 ENDDO
110
111 IF(inacti==5)THEN
112 DO i=1,jlt
113 IF(cs_loc(i)<=nrts) THEN
114 penis(2,cs_loc(i)) =
max(penis(2,cs_loc(i)),half*pene(i))
115 ELSE
116 ni = cs_loc(i)-nrts
118 END IF
119 penim(2,cm_loc(i)) =
max(penim(2,cm_loc(i)),half*pene(i))
120 ENDDO
121 DO i=1,jlt
122 IF(cs_loc(i)<=nrts) THEN
123 pene(i) = pene(i) - penis(1,cs_loc(i)) - penim(1,cm_loc(i))
124 pene(i) =
max(pene(i),zero)
125 IF(pene(i)==zero)stif(i)=zero
126 gapv(i) = gapv(i) - penis(1,cs_loc(i)) - penim(1,cm_loc(i))
127 ELSE
128 ni = cs_loc(i)-nrts
129 pene(i) = pene(i) -
penfi(nin)%P(1,ni) - penim(1,cm_loc(i))
130 pene(i) =
max(pene(i),zero)
131 IF(pene(i)==zero)stif(i)=zero
132 gapv(i) = gapv(i) -
penfi(nin)%P(1,ni) - penim(1,cm_loc(i))
133 END IF
134 END DO
135 ELSE IF(inacti==6)THEN
136 DO i=1,jlt
137 pplus=half*(pene(i)+fiveem2*(gapv(i)-pene(i)))
138 IF(cs_loc(i)<=nrts) THEN
139 penis(2,cs_loc(i)) =
max(penis(2,cs_loc(i)),pplus)
140 ELSE
141 ni = cs_loc(i)-nrts
143 END IF
144 penim(2,cm_loc(i)) =
max(penim(2,cm_loc(i)),pplus)
145 ENDDO
146 DO i=1,jlt
147 IF(cs_loc(i)<=nrts) THEN
148 pene(i) = pene(i) - penis(1,cs_loc(i)) - penim(1,cm_loc(i))
149 pene(i) =
max(pene(i),zero)
150 IF(pene(i)==zero)stif(i)=zero
151 gapv(i) = gapv(i) - penis(1,cs_loc(i)) - penim(1,cm_loc(i))
152 ELSE
153 ni = cs_loc(i)-nrts
154 pene(i) = pene(i) -
penfi(nin)%P(1,ni) - penim(1,cm_loc(i))
155 pene(i) =
max(pene(i),zero)
156 IF(pene(i)==zero)stif(i)=zero
157 gapv(i) = gapv(i) -
penfi(nin)%P(1,ni) - penim(1,cm_loc(i))
158 END IF
159 END DO
160 ENDIF
161
162 DO i=1,jlt
163 gapv(i) = zep9*gapv(i)
164 vx(i) = hs1(i)*vxs1(i) + hs2(i)*vxs2(i)
165 . - hm1(i)*vxm1(i) - hm2(i)*vxm2(i)
166 vy(i) = hs1(i)*vys1(i) + hs2(i)*vys2(i)
167 . - hm1(i)*vym1(i) - hm2(i)*vym2(i)
168 vz(i) = hs1(i)*vzs1(i) + hs2(i)*vzs2(i)
169 . - hm1(i)*vzm1(i) - hm2(i)*vzm2(i)
170 vn(i) = nx(i)*vx(i) + ny(i)*vy(i) + nz(i)*vz(i)
171 h1(i) = hs1(i)*hm1(i)
172 h2(i) = hs1(i)*hm2(i)
173 h3(i) = hs2(i)*hm1(i)
174 h4(i) = hs2(i)*hm2(i)
175 ENDDO
176
177
178 IF(imp_int7>=2)THEN
179 DO i=1,jlt
180 stif(i) = half*stif(i)
181 ENDDO
182 ELSEIF(imp_int7==1)THEN
183 DO i=1,jlt
184 fac = gapv(i)/
max( em10,( gapv(i)-pene(i) ) )
185 IF(( (gapv(i)-pene(i))/gapv(i) )<prec .AND.
186 . stif(i)>zero ) THEN
187 stif(i) = zero
188 pene(i)= zero
189 idesac = 1
190 ELSE
191 stif(i) = half*stif(i) * fac
192 ENDIF
193 ENDDO
194 ELSE
195 DO i=1,jlt
196 fac = gapv(i)/
max( em10,( gapv(i)-pene(i) ) )
197 IF(( (gapv(i)-pene(i))/gapv(i) )<prec .AND.
198 . stif(i)>zero ) THEN
199 stif(i) = zero
200 pene(i)= zero
201 idesac = 1
202 ELSE
203 stif(i) = half*stif(i) * fac
204 ENDIF
205 ENDDO
206 DO i=1,jlt
207 stif(i) = stif(i) * gapv(i) /
208 .
max((gapv(i) - pene(i)),em10)
209 ENDDO
210
211 END IF
212 IF(idesac>0) RETURN
213
214
215
216
217 DO i=1,jlt
218 vtx = vx(i) -vn(i)*nx(i)
219 vty = vy(i) -vn(i)*ny(i)
220 vtz = vz(i) -vn(i)*nz(i)
221 vt = vtx*vtx+vty*vty+vtz*vtz
222 IF (vt>em20) THEN
223 s2=one/sqrt(vt)
224 q(1,1,i)=vtx*s2
225 q(1,2,i)=vty*s2
226 q(1,3,i)=vtz*s2
227 q(3,1,i)=nx(i)
228 q(3,2,i)=ny(i)
229 q(3,3,i)=nz(i)
230 q(2,1,i)=q(3,2,i)*q(1,3,i)-q(3,3,i)*q(1,2,i)
231 q(2,2,i)=q(3,3,i)*q(1,1,i)-q(3,1,i)*q(1,3,i)
232 q(2,3,i)=q(3,1,i)*q(1,2,i)-q(3,2,i)*q(1,1,i)
233 fact(i)=fric
234 ELSE
235 fact(i)=zero
236 ENDIF
237 ENDDO
238 IF (scalk<0) THEN
239 isf=1
240 ELSE
241 isf=0
242 ENDIF
243 facf=abs(scalk)
244 IF (isf==1) THEN
245 DO i=1,jlt
246 IF (vn(i)>zero) THEN
247 fac=stif(i)*facf
248
249 ELSEIF (vn(i)<zero) THEN
250 fac=stif(i)/facf
251
252 ELSE
253 fac=stif(i)
254 ENDIF
255 kn(1,i)=fac*h1(i)
256 kn(2,i)=fac*h2(i)
257 kn(3,i)=fac*h3(i)
258 kn(4,i)=fac*h4(i)
259 fact(i)=fac*fact(i)
260 ENDDO
261 ELSE
262 DO i=1,jlt
263 fac=stif(i)*facf
264 kn(1,i)=fac*h1(i)
265 kn(2,i)=fac*h2(i)
266 kn(3,i)=fac*h3(i)
267 kn(4,i)=fac*h4(i)
268 fact(i)=fac*fact(i)
269 ENDDO
270 ENDIF
271 DO i=1,jlt
272 q11=nx(i)*nx(i)
273 q12=nx(i)*ny(i)
274 q13=nx(i)*nz(i)
275 q22=ny(i)*ny(i)
276 q23=ny(i)*nz(i)
277 q33=nz(i)*nz(i)
278 k1i11(1,1,i)=kn(1,i)*q11
279 k1i11(1,2,i)=kn(1,i)*q12
280 k1i11(1,3,i)=kn(1,i)*q13
281 k1i11(2,2,i)=kn(1,i)*q22
282 k1i11(2,3,i)=kn(1,i)*q23
283 k1i11(3,3,i)=kn(1,i)*q33
284 k1j11(1,1,i)=kn(2,i)*q11
285 k1j11(1,2,i)=kn(2,i)*q12
286 k1j11(1,3,i)=kn(2,i)*q13
287 k1j11(2,2,i)=kn(2,i)*q22
288 k1j11(2,3,i)=kn(2,i)*q23
289 k1j11(3,3,i)=kn(2,i)*q33
290 k2i11(1,1,i)=kn(3,i)*q11
291 k2i11(1,2,i)=kn(3,i)*q12
292 k2i11(1,3,i)=kn(3,i)*q13
293 k2i11(2,2,i)=kn(3,i)*q22
294 k2i11(2,3,i)=kn(3,i)*q23
295 k2i11(3,3,i)=kn(3,i)*q33
296 k2j11(1,1,i)=kn(4,i)*q11
297 k2j11(1,2,i)=kn(4,i)*q12
298 k2j11(1,3,i)=kn(4,i)*q13
299 k2j11(2,2,i)=kn(4,i)*q22
300 k2j11(2,3,i)=kn(4,i)*q23
301 k2j11(3,3,i)=kn(4,i)*q33
302 ENDDO
303
304 DO j=1,3
305 DO k=j,3
306 DO i=1,jlt
307 IF (fact(i)>zero) THEN
308 q1 =q(1,j,i)*q(1,k,i)
309 q2 =q(2,j,i)*q(2,k,i)
310 fac=fact(i)*(q1+q2)
311 kt1=fac*h1(i)
312 k1i11(j,k,i)=k1i11(j,k,i)+kt1
313 kt2=fac*h2(i)
314 k1j11(j,k,i)=k1j11(j,k,i)+kt2
315 kt3=fac*h3(i)
316 k2i11(j,k,i)=k2i11(j,k,i)+kt3
317 kt4=fac*h4(i)
318 k2j11(j,k,i)=k2j11(j,k,i)+kt4
319 ENDIF
320 ENDDO
321 ENDDO
322 ENDDO
323
324 DO j=1,3
325 DO k=j,3
326 DO i=1,jlt
327 k1i12(j,k,i)=-k1i11(j,k,i)
328 k1j12(j,k,i)=-k1j11(j,k,i)
329 k2i12(j,k,i)=-k2i11(j,k,i)
330 k2j12(j,k,i)=-k2j11(j,k,i)
331 ENDDO
332 ENDDO
333 ENDDO
334 DO j=1,3
335 DO k=j+1,3
336 DO i=1,jlt
337 k1i12(k,j,i)=-k1i11(j,k,i)
338 k1j12(k,j,i)=-k1j11(j,k,i)
339 k2i12(k,j,i)=-k2i11(j,k,i)
340 k2j12(k,j,i)=-k2j11(j,k,i)
341 ENDDO
342 ENDDO
343 ENDDO
344
345 DO i=1,jlt
346 off(i)=one
347 ENDDO
348
349 IF (nspmd>1) THEN
351 DO i=1,jlt
352 IF(cs_loc(i)>nrts) THEN
353 nn=cs_loc(i)-nrts
355
356 nn1 = ns
357 ffi(1,nn1)=zero
358 ffi(2,nn1)=zero
359 ffi(3,nn1)=zero
360 dfi(1,nn1)=zero
361 dfi(2,nn1)=zero
362 dfi(3,nn1)=zero
363 nn2 = nn1 + 1
364 ffi(1,nn2)=zero
365 ffi(2,nn2)=zero
366 ffi(3,nn2)=zero
367 dfi(1,nn2)=zero
368 dfi(2,nn2)=zero
369 dfi(3,nn2)=zero
370 ENDIF
371 ENDDO
372 ELSE
373 jltf = 0
374 DO i=1,jlt
375 IF(cs_loc(i)>nrts) THEN
376 jltf = jltf + 1
378 nn=cs_loc(i)-nrts
380 stifs(ne)=stif(i)
381 h_e(1,ne)=hs1(i)
382 h_e(2,ne)=hs2(i)
383 h_e(3,ne)=hm1(i)
384 h_e(4,ne)=hm2(i)
385 n_e(1,ne)=nx(i)
386 n_e(2,ne)=ny(i)
387 n_e(3,ne)=nz(i)
388
389 nn1 = ns
390 ffi(1,nn1)=zero
391 ffi(2,nn1)=zero
392 ffi(3,nn1)=zero
393 dfi(1,nn1)=zero
394 dfi(2,nn1)=zero
395 dfi(3,nn1)=zero
396 nn2 = nn1 + 1
397 ffi(1,nn2)=zero
398 ffi(2,nn2)=zero
399 ffi(3,nn2)=zero
400 dfi(1,nn2)=zero
401 dfi(2,nn2)=zero
402 dfi(3,nn2)=zero
403 ENDIF
404 ENDDO
405 ENDIF
406 ENDIF
407
408 RETURN
integer, dimension(:), allocatable shf_int
type(int_pointer2), dimension(:), allocatable ind_int
type(real_pointer2), dimension(:), allocatable penfi