40
41
42
43
44#include "implicit_f.inc"
45
46
47
48#include "mvsiz_p.inc"
49#include "sms_c.inc"
50
51
52
53 INTEGER JLT,JLT_NEW,IGAP,INTFRIC
54 INTEGER CAND_S(*),CAND_M(*),
55 . N1(*),N2(*),M1(*),M2(*), NSMS(*),INDEX(*),
56 . IPARTFRICSI(MVSIZ),IPARTFRICMI(MVSIZ)
57 my_real ,
INTENT(IN) :: drad, dgapload
59 . h1s(*),h2s(*),h1m(*),h2m(*),nx(*),ny(*),nz(*),stif(*),
60 . xxs1(*) ,xxs2(*) ,xys1(*) ,xys2(*) ,
61 . xzs1(*) ,xzs2(*) ,xxm1(*) ,xxm2(*) ,xym1(*),
62 . xym2(*) ,xzm1(*) ,xzm2(*) ,vxs1(*) ,vxs2(*),
63 . vys1(*) ,vys2(*) ,vzs1(*) ,vzs2(*) ,vxm1(*),
64 . vxm2(*) ,vym1(*) ,vym2(*) ,vzm1(*) ,vzm2(*),
65 . ms1(*) ,ms2(*) ,mm1(*) ,mm2(*), gapv(*)
66
67
68
69 INTEGER I
71 . pene2(mvsiz),
72 . xs12,ys12,zs12,xm12,ym12,zm12,xa,xb,
73 . xs2,xm2,xsm,xs2m2,ys2,ym2,ysm,ys2m2,zs2,zm2,zsm,zs2m2,
74 . xx,yy,zz,als,alm,det,
75 . gap2,drad2
76
77 jlt_new = 0
78 drad2 = drad*drad
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128 DO i=1,jlt
129 xs12 = xxs2(i)-xxs1(i)
130 ys12 = xys2(i)-xys1(i)
131 zs12 = xzs2(i)-xzs1(i)
132 xs2 = xs12*xs12 + ys12*ys12 + zs12*zs12
133 xm12 = xxm2(i)-xxm1(i)
134 ym12 = xym2(i)-xym1(i)
135 zm12 = xzm2(i)-xzm1(i)
136 xm2 = xm12*xm12 + ym12*ym12 + zm12*zm12
137 xsm = - (xs12*xm12 + ys12*ym12 + zs12*zm12)
138 xs2m2 = xxm2(i)-xxs2(i)
139 ys2m2 = xym2(i)-xys2(i)
140 zs2m2 = xzm2(i)-xzs2(i)
141
142 xa = xs12*xs2m2 + ys12*ys2m2 + zs12*zs2m2
143 xb = -xm12*xs2m2 - ym12*ys2m2 - zm12*zs2m2
144 det = xm2*xs2 - xsm*xsm
146
147 h1m(i) = (xa*xsm-xb*xs2) / det
148
151 h1m(i)=
min(one,
max(zero,h1m(i)))
152 h1s(i) = -(xa + h1m(i)*xsm) / xs2
153 h1s(i)=
min(one,
max(zero,h1s(i)))
154 h1m(i) = -(xb + h1s(i)*xsm) / xm2
155 h1m(i)=
min(one,
max(zero,h1m(i)))
156
157 h2s(i) = one -h1s(i)
158 h2m(i) = one -h1m(i)
159
160
161
162 nx(i) = h1s(i)*xxs1(i) + h2s(i)*xxs2(i)
163 . - h1m(i)*xxm1(i) - h2m(i)*xxm2(i)
164 ny(i) = h1s(i)*xys1(i) + h2s(i)*xys2(i)
165 . - h1m(i)*xym1(i) - h2m(i)*xym2(i)
166 nz(i) = h1s(i)*xzs1(i) + h2s(i)*xzs2(i)
167 . - h1m(i)*xzm1(i) - h2m(i)*xzm2(i)
168 gap2 = (gapv(i)+dgapload)*(gapv(i)+dgapload)
169 pene2(i) =
max(gap2,drad2) - nx(i)*nx(i) - ny(i)*ny(i) - nz(i)*nz(i)
170 pene2(i) =
max(zero,pene2(i))
171
172 ENDDO
173
174 IF(idtmins/=2)THEN
175 IF(intfric == 0) THEN
176 DO i=1,jlt
177 IF(pene2(i)/=zero.AND.stif(i)/=zero)THEN
178 jlt_new = jlt_new + 1
179 cand_s(jlt_new) = cand_s
180 cand_m(jlt_new) = cand_m(i)
181 n1(jlt_new) = n1(i)
182 n2(jlt_new) = n2(i)
183 m1(jlt_new) = m1(i)
184 m2(jlt_new) = m2(i)
185 h1s(jlt_new) = h1s(i)
186 h2s(jlt_new) = h2s(i)
187 h1m(jlt_new) = h1m(i)
188
189 nx(jlt_new) = nx(i)
190 ny(jlt_new) = ny(i)
191 nz(jlt_new) = nz(i)
192 stif(jlt_new)= stif(i)
193 gapv(jlt_new)= gapv(i)
194 vxs1(jlt_new) = vxs1(i)
195 vys1(jlt_new) = vys1(i)
196 vzs1(jlt_new) = vzs1(i)
197 vxs2(jlt_new) = vxs2(i)
198 vys2(jlt_new) = vys2(i)
199 vzs2(jlt_new) = vzs2(i)
200 vxm1(jlt_new) = vxm1(i)
201 vym1(jlt_new) = vym1(i)
202 vzm1(jlt_new) = vzm1(i)
203 vxm2(jlt_new) = vxm2(i)
204 vym2(jlt_new) = vym2(i)
205 vzm2(jlt_new) = vzm2(i)
206 ms1(jlt_new) = ms1(i)
207 ms2(jlt_new) = ms2(i)
208 mm1(jlt_new) = mm1(i)
209 mm2(jlt_new) = mm2(i)
210 index(jlt_new) = index(i)
211 ENDIF
212 ENDDO
213 ELSE
214 DO i=1,jlt
215 IF(pene2(i)/=zero.AND.stif(i)/=zero)THEN
216 jlt_new = jlt_new + 1
217 cand_s(jlt_new) = cand_s(i)
218 cand_m(jlt_new) = cand_m(i)
219 n1(jlt_new) = n1(i)
220 n2(jlt_new) = n2(i)
221 m1(jlt_new) = m1(i)
222 m2(jlt_new) = m2(i)
223 h1s(jlt_new) = h1s(i)
224 h2s(jlt_new) = h2s(i)
225 h1m(jlt_new) = h1m(i)
226 h2m(jlt_new) = h2m(i)
227 nx(jlt_new) = nx(i)
228 ny(jlt_new) = ny(i)
229 nz(jlt_new) = nz(i)
230 stif(jlt_new)= stif(i)
231 gapv(jlt_new)= gapv(i)
232 vxs1(jlt_new) = vxs1(i)
233 vys1(jlt_new) = vys1(i)
234 vzs1(jlt_new) = vzs1(i)
235 vxs2(jlt_new) = vxs2(i)
236 vys2(jlt_new) = vys2(i)
237 vzs2(jlt_new) = vzs2(i)
238 vxm1(jlt_new) = vxm1(i)
239 vym1(jlt_new) = vym1(i)
240 vzm1(jlt_new) = vzm1(i)
241 vxm2(jlt_new) = vxm2(i)
242 vym2(jlt_new) = vym2(i)
243 vzm2(jlt_new) = vzm2(i)
244 ms1(jlt_new) = ms1(i)
245 ms2(jlt_new) = ms2(i)
246 mm1(jlt_new) = mm1(i)
247 mm2(jlt_new) = mm2(i)
248 ipartfricsi(jlt_new)=ipartfricsi(i)
249 ipartfricmi(jlt_new)=ipartfricmi(i)
250 index(jlt_new) = index(i)
251 ENDIF
252 ENDDO
253 ENDIF
254 ELSE
255 IF(intfric == 0) THEN
256 DO i=1,jlt
257 IF(pene2(i)/=zero.AND.stif(i)/=zero)THEN
258 jlt_new = jlt_new + 1
259 cand_s(jlt_new) = cand_s(i)
260 cand_m(jlt_new) = cand_m(i)
261 n1(jlt_new) = n1(i)
262 n2(jlt_new) = n2(i)
263 m1(jlt_new) = m1(i)
264 m2(jlt_new) = m2(i)
265 h1s(jlt_new) = h1s(i)
266 h2s(jlt_new) = h2s(i)
267 h1m(jlt_new) = h1m(i)
268 h2m(jlt_new) = h2m(i)
269 nx(jlt_new) = nx(i)
270 ny(jlt_new) = ny(i)
271 nz(jlt_new) = nz(i)
272 stif(jlt_new)= stif(i)
273 gapv(jlt_new)= gapv(i)
274 vxs1(jlt_new) = vxs1(i)
275 vys1(jlt_new) = vys1(i)
276 vzs1(jlt_new) = vzs1(i)
277 vxs2(jlt_new) = vxs2(i)
278 vys2(jlt_new) = vys2(i)
279 vzs2(jlt_new) = vzs2(i)
280 vxm1(jlt_new) = vxm1(i)
281 vym1(jlt_new) = vym1(i)
282 vzm1(jlt_new) = vzm1(i)
283 vxm2(jlt_new) = vxm2(i)
284 vym2(jlt_new) = vym2(i)
285 vzm2(jlt_new) = vzm2(i)
286 ms1(jlt_new) = ms1(i)
287 ms2(jlt_new) = ms2(i)
288 mm1(jlt_new) = mm1(i)
289 mm2(jlt_new) = mm2(i)
290 nsms(jlt_new)= nsms(i)
291 index(jlt_new) = index(i)
292 ENDIF
293 ENDDO
294 ELSE
295 DO i=1,jlt
296 IF(pene2(i)/=zero.AND.stif(i)/=zero)THEN
297 jlt_new = jlt_new + 1
298 cand_s(jlt_new) = cand_s(i)
299 cand_m(jlt_new) = cand_m(i)
300 n1(jlt_new) = n1(i)
301 n2(jlt_new) = n2(i)
302 m1(jlt_new) = m1(i)
303 m2(jlt_new) = m2(i)
304 h1s(jlt_new) = h1s(i)
305 h2s(jlt_new) = h2s(i)
306 h1m(jlt_new) = h1m(i)
307 h2m(jlt_new) = h2m(i)
308 nx(jlt_new) = nx(i)
309 ny(jlt_new) = ny(i)
310 nz(jlt_new) = nz(i)
311 stif(jlt_new)= stif(i)
312 gapv(jlt_new)= gapv(i)
313 vxs1(jlt_new) = vxs1(i)
314 vys1(jlt_new) = vys1(i)
315 vzs1(jlt_new) = vzs1(i)
316 vxs2(jlt_new) = vxs2(i)
317 vys2(jlt_new) = vys2(i)
318 vzs2(jlt_new) = vzs2(i)
319 vxm1(jlt_new) = vxm1(i)
320 vym1(jlt_new) = vym1(i)
321 vzm1(jlt_new) = vzm1(i)
322 vxm2(jlt_new) = vxm2(i)
323 vym2(jlt_new) = vym2(i)
324 vzm2(jlt_new) = vzm2(i)
325 ms1(jlt_new) = ms1(i)
326 ms2(jlt_new) = ms2(i)
327 mm1(jlt_new) = mm1(i)
328 mm2(jlt_new) = mm2(i)
329 nsms(jlt_new)= nsms(i)
330 ipartfricsi(jlt_new)=ipartfricsi(i)
331 ipartfricmi(jlt_new)=ipartfricmi(i)
332 index(jlt_new) = index(i)
333 ENDIF
334 ENDDO
335 ENDIF
336 END IF
337
338 RETURN