56
57
58
59
60
61
62
63
65 USE multi_fvm_mod
67
68
69
70#include "implicit_f.inc"
71#include "comlock.inc"
72#include "com04_c.inc"
73#include "tabsiz_c.inc"
74
75
76
77#include "mvsiz_p.inc"
78
79
80
81 INTEGER JLT, JLT_NEW, CAND_N(*),CN_LOC(MVSIZ),
82 . CAND_E(*),CE_LOC(MVSIZ), NSVG(MVSIZ), KINI(*),ICURV
83 INTEGER IX1(MVSIZ), IX2(MVSIZ), IX3(MVSIZ), IX4(MVSIZ),
84 . INDEX(MVSIZ)
85 INTEGER, INTENT(IN) :: ITAB(NUMNOD),IGAP
86 INTEGER, INTENT(in) :: S_XCELL_REMOTE
88 . nx1(mvsiz), nx2(mvsiz), nx3(mvsiz), nx4(mvsiz),
89 . ny1(mvsiz), ny2(mvsiz), ny3(mvsiz), ny4(mvsiz),
90 . nz1(mvsiz), nz2(mvsiz), nz3(mvsiz), nz4(mvsiz),
91 . lb1(mvsiz), lb2(mvsiz), lb3(mvsiz), lb4(mvsiz),
92 . lc1(mvsiz), lc2(mvsiz), lc3(mvsiz), lc4(mvsiz),
93 . kb1(mvsiz), kb2(mvsiz), kb3(mvsiz), kb4(mvsiz),
94 . kc1(mvsiz), kc2(mvsiz), kc3(mvsiz), kc4(mvsiz),
95 . x1(mvsiz), x2(mvsiz), x3(mvsiz), x4(mvsiz),
96 . y1(mvsiz), y2(mvsiz), y3(mvsiz), y4(mvsiz),
97 . z1(mvsiz), z2(mvsiz), z3(mvsiz), z4(mvsiz),
98 . xi(mvsiz), yi(mvsiz), zi(mvsiz), stif(mvsiz),
99 . p1(mvsiz), p2(mvsiz), p3(mvsiz), p4(mvsiz),
100 . gapv(mvsiz), cand_p(*),
101 . vxi(mvsiz), vyi(mvsiz), vzi(mvsiz), msi(mvsiz),
102 . xcell(3,sxcell)
103 TYPE(MULTI_FVM_STRUCT), INTENT(IN) :: MULTI_FVM
104 TYPE(t_connectivity), INTENT(IN) :: ALE_NE_CONNECT
105 my_real,
DIMENSION(S_XCELL_REMOTE),
INTENT(in) :: xcell_remote
106
107
108
109 INTEGER I, IG, J, IEL, IAD1, IAD2, ELEM_ID
111 . x0(mvsiz), y0(mvsiz), z0(mvsiz),
112 . al1(mvsiz), al2(mvsiz), al3(mvsiz), al4(mvsiz),
113 . x01(mvsiz), x02(mvsiz), x03(mvsiz), x04(mvsiz),
114 . y01(mvsiz), y02(mvsiz), y03(mvsiz), y04(mvsiz),
115 . z01(mvsiz), z02(mvsiz), z03(mvsiz), z04(mvsiz),
116 . xi1(mvsiz), xi2(mvsiz), xi3(mvsiz), xi4(mvsiz),
117 . yi1(mvsiz), yi2(mvsiz), yi3(mvsiz), yi4(mvsiz),
118 . zi1(mvsiz), zi2(mvsiz), zi3(mvsiz), zi4(mvsiz),
119 . pene2(mvsiz),
120 . hlb1(mvsiz), hlc1(mvsiz), hlb2(mvsiz),hlc2(mvsiz),
121 . hlb3(mvsiz), hlc3(mvsiz), hlb4(mvsiz),hlc4(mvsiz)
123 . s2,d1,d2,d3,d4,dl,
124 . x12,x23,x34,x41,xi0,sx1,sx2,sx3,sx4,sx0,
125 . y12,y23,y34,y41,yi0,sy1,sy2,sy3,sy4,sy0,
126 . z12,z23,z34,z41,zi0,sz1,sz2,sz3,sz4,sz0,
127 . gap2(mvsiz), ds2,t1,t2,t3,
128 . al1num,al2num,al3num,al4num,al1den,al2den,al3den,al4den,
129 . x23d,y23d,z23d,x34d,y34d,z34d,x41d,y41d,z41d,
130 . x12d,y12d,z12d,gap2d,xi0d,yi0d,zi0d,s2d, la, hla, aaa,
131 . xi0v(mvsiz), yi0v(mvsiz), zi0v(mvsiz)
133 . nnx1(mvsiz), nnx2(mvsiz), nnx3(mvsiz), nnx4(mvsiz),
134 . nny1(mvsiz), nny2(mvsiz), nny3(mvsiz), nny4(mvsiz),
135 . nnz1(mvsiz), nnz2(mvsiz), nnz3(mvsiz), nnz4(mvsiz)
136 LOGICAL IS_GAP_CONSTANT, IS_GAP_VARIABLE
137 INTEGER :: NODE_ID
138
139
140
141
142 is_gap_constant = .true.
143 is_gap_variable = .false.
144 IF(igap == 1)THEN
145 is_gap_constant = .false.
146 is_gap_variable = .true.
147 ENDIF
148
149
150 DO i=1,jlt
151 IF(ix3(i) /= ix4(i))THEN
152
153 x0(i) = fourth*(x1(i)+x2(i)+x3(i)+x4(i))
154 y0(i) = fourth*(y1(i)+y2(i)+y3(i)+y4(i))
155 z0(i) = fourth*(z1(i)+z2(i)+z3(i)+z4(i))
156 ELSE
157
158 x0(i) = x3(i)
159 y0(i) = y3(i)
160 z0(i) = z3(i)
161 ENDIF
162 ENDDO
163
164
165 DO i=1,jlt
166
167 x01(i) = x1(i) - x0(i)
168 y01(i) = y1(i) - y0(i)
169 z01(i) = z1(i) - z0(i)
170
171 x02(i) = x2(i) - x0(i)
172 y02(i) = y2(i) - y0(i)
173 z02(i) = z2(i) - z0(i)
174
175 x03(i) = x3(i) - x0(i)
176 y03(i) = y3(i) - y0(i)
177 z03(i) = z3(i) - z0(i)
178
179 x04(i) = x4(i) - x0(i)
180 y04(i) = y4(i) - y0(i)
181 z04(i) = z4(i) - z0(i)
182
183 xi0v(i) = x0(i) - xi(i)
184 yi0v(i) = y0(i) - yi(i)
185 zi0v(i) = z0(i) - zi(i)
186
187 xi1(i) = x1(i) - xi(i)
188 yi1(i) = y1(i) - yi(i)
189 zi1(i) = z1(i) - zi(i)
190
191 xi2(i) = x2(i) - xi(i)
192 yi2(i) = y2(i) - yi(i)
193 zi2(i) = z2(i) - zi(i)
194
195 xi3(i) = x3(i) - xi(i)
196 yi3(i) = y3(i) - yi(i)
197 zi3(i) = z3(i) - zi(i)
198
199 xi4(i) = x4(i) - xi(i)
200 yi4(i) = y4(i) - yi(i)
201 zi4(i) = z4(i) - zi(i)
202 ENDDO
203
204 DO i=1,jlt
205
206 sx1 = yi0v(i)*zi1(i) - zi0v(i)*yi1(i)
207 sy1 = zi0v(i)*xi1(i) - xi0v(i)*zi1(i)
208 sz1 = xi0v(i)*yi1(i) - yi0v(i)*xi1(i)
209
210 sx2 = yi0v(i)*zi2(i) - zi0v(i)*yi2(i)
211 sy2 = zi0v(i)*xi2(i) - xi0v(i)*zi2(i)
212 sz2 = xi0v(i)*yi2(i) - yi0v(i)*xi2(i)
213
214 sx3 = yi0v(i)*zi3(i) - zi0v(i)*yi3(i)
215 sy3 = zi0v(i)*xi3(i) - xi0v(i)*zi3(i)
216 sz3 = xi0v(i)*yi3(i) - yi0v(i)*xi3(i)
217
218 sx4 = yi0v(i)*zi4(i) - zi0v(i)*yi4(i)
219 sy4 = zi0v(i)*xi4(i) - xi0v(i)*zi4(i)
220 sz4 = xi0v(i)*yi4(i) - yi0v(i)*xi4(i)
221
222
223
224 sx0 = y01(i)*z02(i) - z01(i)*y02(i)
225 sy0 = z01(i)*x02(i) - x01(i)*z02(i)
226 sz0 = x01(i)*y02(i) - y01(i)*x02(i)
227 s2 = one/
max(em30,sx0*sx0 + sy0*sy0 + sz0*sz0)
228 nnx1(i) = sx0
229 nny1(i) = sy0
230 nnz1(i) = sz0
231 lb1(i) = -(sx0*sx2 + sy0*sy2 + sz0*sz2) * s2
232 lc1(i) = (sx0*sx1 + sy0*sy1 + sz0*sz1) * s2
233
234
235
236 sx0 = y02(i)*z03(i) - z02(i)*y03(i)
237 sy0 = z02(i)*x03(i) - x02(i)*z03(i)
238 sz0 = x02(i)*y03(i) - y02(i)*x03(i)
239 s2 = 1./
max(em30,sx0*sx0 + sy0*sy0 + sz0*sz0)
240 nnx2(i) = sx0
241 nny2(i) = sy0
242 nnz2(i) = sz0
243 lb2(i) = -(sx0*sx3 + sy0*sy3 + sz0*sz3) * s2
244 lc2(i) = (sx0*sx2 + sy0*sy2 + sz0*sz2) * s2
245
246
247
248 sx0 = y03(i)*z04(i) - z03(i)*y04(i)
249 sy0 = z03(i)*x04(i) - x03(i)*z04(i)
250 sz0 = x03(i)*y04(i) - y03(i)*x04(i)
251 s2 = one/
max(em30,sx0*sx0 + sy0*sy0 + sz0*sz0)
252 nnx3(i) = sx0
253 nny3(i) = sy0
254 nnz3(i) = sz0
255 lb3(i) = -(sx0*sx4 + sy0*sy4 + sz0*sz4) * s2
256 lc3(i) = (sx0*sx3 + sy0*sy3 + sz0*sz3) * s2
257
258
259
260 sx0 = y04(i)*z01(i) - z04(i)*y01(i)
261 sy0 = z04(i)*x01(i) - x04(i)*z01(i)
262 sz0 = x04(i)*y01(i) - y04(i)*x01(i)
263 s2 = one/
max(em30,sx0*sx0 + sy0*sy0 + sz0*sz0)
264 nnx4(i) = sx0
265 nny4(i) = sy0
266 nnz4(i) = sz0
267 lb4(i) = -(sx0*sx1 + sy0*sy1 + sz0*sz1) * s2
268 lc4(i) = (sx0*sx4 + sy0*sy4 + sz0*sz4) * s2
269 ENDDO
270 IF(igap == 1)THEN
271 IF (.NOT.multi_fvm%IS_USED) THEN
272 DO i=1,jlt
273 node_id = nsvg(i)
274
275
276 IF(node_id>0) THEN
277 iad1 = ale_ne_connect%IAD_CONNECT(node_id)
278 iad2 = ale_ne_connect%IAD_CONNECT(node_id + 1) - 1
279 dl = zero
280 DO iel = iad1, iad2
281 elem_id = ale_ne_connect%CONNECTED(iel)
282 dl=
max(dl, xcell(1,elem_id))
283 ENDDO
284 ELSE
285 node_id = abs(node_id)
286 dl = xcell_remote(node_id)
287 ENDIF
288
289 gapv(i) = three_half*dl
290 enddo
291 ELSEIF(multi_fvm%IS_USED) THEN
292 DO i=1,jlt
293 IF(nsvg(i)>0) THEN
294 elem_id = nsvg(i) - numnod
295 dl=xcell(1,elem_id)
296 ELSE
297 node_id = abs(nsvg(i))
298 dl = xcell_remote(node_id)
299 ENDIF
300
301 gapv(i) = three_half*dl
302 enddo
303 ENDIF
304 ELSE
305 gapv(1:jlt)=gapv(1)
306 ENDIF
307
308 IF(is_gap_constant)THEN
309
310 aaa =gapv(1)*gapv(1)
311 DO i=1,jlt
312 gap2(i)=aaa
313 ENDDO
314 ELSE
315 DO i=1,jlt
316 gap2(i)=gapv(i)*gapv(i)
317 ENDDO
318 ENDIF
319
320 DO i=1,jlt
321 aaa = one/
max(em30,x01(i)*x01(i)+y01(i)*y01(i)+z01(i)*z01(i))
322 hlc1(i)= lc1(i)*abs(lc1(i))*aaa
323 hlb4(i)= lb4(i)*abs(lb4(i))*aaa
324 al1(i) = -(xi0v(i)*x01(i)+yi0v(i)*y01(i)+zi0v(i)*z01(i))*aaa
325 al1(i) =
max(zero,
min(one,al1(i)))
326 aaa = one/
max(em30,x02(i)*x02(i)+y02(i)*y02(i)+z02(i)*z02(i))
327 hlc2(i)= lc2(i)*abs(lc2(i))*aaa
328 hlb1(i)= lb1(i)*abs(lb1(i))*aaa
329 al2(i) = -(xi0v(i)*x02(i)+yi0v(i)*y02(i)+zi0v(i)*z02(i))*aaa
330 al2(i) =
max(zero,
min(one,al2(i)))
331 aaa = one/
max(em30,x03(i)*x03(i)+y03(i)*y03(i)+z03(i)*z03(i))
332 hlc3(i)= lc3(i)*abs(lc3(i))*aaa
333 hlb2(i)= lb2(i)*abs(lb2(i))*aaa
334 al3(i) = -(xi0v(i)*x03(i)+yi0v(i)*y03(i)+zi0v(i)*z03(i))*aaa
335 al3(i) =
max(zero,
min(one,al3(i)))
336 aaa = one/
max(em30,x04(i)*x04(i)+y04(i)*y04(i)+z04(i)*z04(i))
337 hlc4(i)= lc4(i)*abs(lc4(i))*aaa
338 hlb3(i)= lb3(i)*abs(lb3(i))*aaa
339 al4(i) = -(xi0v(i)*x04(i)+yi0v(i)*y04(i)+zi0v(i)*z04(i))*aaa
340 al4(i) =
max(zero,
min(one,al4(i)))
341 ENDDO
342
343
344
345
346
347 DO i=1,jlt
348 x12 = x2(i) - x1(i)
349 y12 = y2(i) - y1(i)
350 z12 = z2(i) - z1(i)
351 la = one - lb1(i) - lc1(i)
352 aaa = one /
max(em20,x12*x12+y12*y12+z12*z12)
353 hla= la*abs(la) * aaa
354 IF(la < zero.AND.hla <= hlb1(i).AND.hla <= hlc1(i))THEN
355 lb1(i) = (xi2(i)*x12+yi2(i)*y12+zi2(i)*z12) * aaa
356 lb1(i) =
max(zero,
min(one,lb1(i)))
357 lc1(i) = one - lb1(i)
358 ELSEIF(lb1(i) < zero.AND.hlb1(i) <= hlc1(i).AND.hlb1(i) <= hla)THEN
359 lb1(i) = zero
360 lc1(i) = al2(i)
361 ELSEIF(lc1(i) < zero.AND.hlc1(i) <= hla.AND.hlc1(i) <= hlb1(i))THEN
362 lc1(i) = zero
363 lb1(i) = al1(i)
364 ENDIF
365 ENDDO
366
367
368
369
370
371 DO i=1,jlt
372 x23 = x3(i) - x2(i)
373 y23 = y3(i) - y2(i)
374 z23 = z3(i) - z2(i)
375 la = one - lb2(i) - lc2(i)
376 aaa = one /
max(em20,x23*x23+y23*y23+z23*z23)
377 hla= la*abs(la) * aaa
378 IF(la < zero.AND.hla <= hlb2(i).AND.hla <= hlc2(i))THEN
379 lb2(i) = (xi3(i)*x23+yi3(i)*y23+zi3(i)*z23)*aaa
380 lb2(i) =
max(zero,
min(one,lb2(i)))
381 lc2(i) = one - lb2(i)
382 ELSEIF(lb2(i) < zero.AND.hlb2(i) <= hlc2(i).AND.hlb2(i) <= hla)THEN
383 lb2(i) = zero
384 lc2(i) = al3(i)
385 ELSEIF(lc2(i) < zero.AND.hlc2(i) <= hla.AND.hlc2(i) <= hlb2(i))THEN
386 lc2(i) = zero
387 lb2(i) = al2(i)
388 ENDIF
389 ENDDO
390
391
392
393
394
395 DO i=1,jlt
396 x34 = x4(i) - x3(i)
397 y34 = y4(i) - y3(i)
398 z34 = z4(i) - z3(i)
399 la = one - lb3(i) - lc3(i)
400 aaa = one /
max(em20,x34*x34+y34*y34+z34*z34)
401 hla= la*abs(la) * aaa
402 IF(la < zero.AND.hla <= hlb3(i).AND.hla <= hlc3(i))THEN
403 lb3(i) = (xi4(i)*x34+yi4(i)*y34+zi4(i)*z34)*aaa
404 lb3(i) =
max(zero,
min(one,lb3(i)))
405 lc3(i) = one - lb3(i)
406 ELSEIF(lb3(i) < zero.AND.hlb3(i) <= hlc3(i).AND.hlb3(i) <= hla)THEN
407 lb3(i) = zero
408 lc3(i) = al4(i)
409 ELSEIF(lc3(i) < zero.AND.hlc3(i) <= hla.AND.hlc3(i) <= hlb3(i))THEN
410 lc3(i) = zero
411 lb3(i) = al3(i)
412 ENDIF
413 ENDDO
414
415
416
417
418
419 DO i=1,jlt
420 x41 = x1(i) - x4(i)
421 y41 = y1(i) - y4(i)
422 z41 = z1(i) - z4(i)
423 la = one - lb4(i) - lc4(i)
424 aaa = one /
max(em20,x41*x41+y41*y41+z41*z41)
425 hla= la*abs(la) * aaa
426 IF(la < zero.AND.hla <= hlb4(i).AND.hla <= hlc4(i))THEN
427 lb4(i) = (xi1(i)*x41+yi1(i)*y41+zi1(i)*z41)*aaa
428 lb4(i) =
max(zero,
min(one,lb4(i)))
429 lc4(i) = one - lb4(i)
430 ELSEIF(lb4(i) < zero.AND.hlb4(i) <= hlc4(i).AND.hlb4(i) <= hla)THEN
431 lb4(i) = zero
432 lc4(i) = al1(i)
433 ELSEIF(lc4(i) < zero.AND.hlc4(i) <= hla.AND.hlc4(i) <= hlb4(i))THEN
434 lc4(i) = zero
435 lb4(i) = al4(i)
436 ENDIF
437 ENDDO
438
439
440 DO i=1,jlt
441
442 nx1(i) = xi(i)-(x0(i) + lb1(i)*x01(i) + lc1(i)*x02(i))
443 ny1(i) = yi(i)-(y0(i) + lb1(i)*y01(i) + lc1(i)*y02(i))
444 nz1(i) = zi(i)-(z0(i) + lb1(i)*z01(i) + lc1(i)*z02(i))
445 p1(i) = nx1(i)*nx1(i) + ny1(i)*ny1(i) +nz1(i)*nz1(i)
446
447 nx2(i) = xi(i)-(x0(i) + lb2(i)*x02(i) + lc2(i)*x03(i))
448 ny2(i) = yi(i)-(y0(i) + lb2(i)*y02(i) + lc2(i)*y03(i))
449 nz2(i) = zi(i)-(z0(i) + lb2(i)*z02(i) + lc2(i)*z03(i))
450 p2(i) = nx2(i)*nx2(i) + ny2(i)*ny2(i) +nz2(i)*nz2(i)
451
452 nx3(i) = xi(i)-(x0(i) + lb3(i)*x03(i) + lc3(i)*x04(i))
453 ny3(i) = yi(i)-(y0(i) + lb3(i)*y03(i) + lc3(i)*y04(i))
454 nz3(i) = zi(i)-(z0(i) + lb3(i)*z03(i) + lc3(i)*z04(i))
455 p3(i) = nx3(i)*nx3(i) + ny3(i)*ny3(i) +nz3(i)*nz3(i)
456
457 nx4(i) = xi(i)-(x0(i) + lb4(i)*x04(i) + lc4(i)*x01(i))
458 ny4(i) = yi(i)-(y0(i) + lb4(i)*y04(i) + lc4(i)*y01(i))
459 nz4(i) = zi(i)-(z0(i) + lb4(i)*z04(i) + lc4(i)*z01(i))
460 p4(i) = nx4(i)*nx4(i) + ny4(i)*ny4(i) +nz4(i)*nz4(i)
461
462 d1 =
max(zero, gap2(i) - p1(i))
463 d2 =
max(zero, gap2(i) - p2(i))
464 d3 =
max(zero, gap2(i) - p3(i))
465 d4 =
max(zero, gap2(i) - p4(i))
466 pene2(i) =
max(d1,d2,d3,d4)
467
468 ENDDO
469
470
471
472 DO i=1,jlt
473 IF(pene2(i) == zero.OR.stif(i) == zero)THEN
474 cand_p(index(i))=0
475 ENDIF
476 ENDDO
477 DO i=1,jlt
478 IF(pene2(i) /= zero.AND.stif(i) /= zero)THEN
479 jlt_new = jlt_new + 1
480 cn_loc(jlt_new) = cand_n(i)
481 ce_loc(jlt_new) = cand_e(i)
482 ix1(jlt_new) = ix1(i)
483 ix2(jlt_new) = ix2(i)
484 ix3(jlt_new) = ix3(i)
485 ix4(jlt_new) = ix4(i)
486 nsvg(jlt_new) = nsvg(i)
487 nx1(jlt_new) = nnx1(i)
488 nx2(jlt_new) = nnx2(i)
489 nx3(jlt_new) = nnx3(i)
490 nx4(jlt_new) = nnx4(i)
491 ny1(jlt_new) = nny1(i)
492 ny2(jlt_new) = nny2(i)
493 ny3(jlt_new) = nny3(i)
494 ny4(jlt_new) = nny4(i)
495 nz1(jlt_new) = nnz1(i)
496 nz2(jlt_new) = nnz2(i)
497 nz3(jlt_new) = nnz3(i)
498 nz4(jlt_new) = nnz4(i)
499 p1(jlt_new) = p1(i)
500 p2(jlt_new) = p2(i)
501 p3(jlt_new) = p3(i)
502 p4(jlt_new) = p4(i)
503 lb1(jlt_new) = lb1(i)
504 lb2(jlt_new) = lb2(i)
505 lb3(jlt_new) = lb3(i)
506 lb4(jlt_new) = lb4(i)
507 lc1(jlt_new) = lc1(i)
508 lc2(jlt_new) = lc2(i)
509 lc3(jlt_new) = lc3(i)
510 lc4(jlt_new) = lc4(i)
511 stif(jlt_new) = stif(i)
512 gapv(jlt_new) = gapv(i)
513 index(jlt_new) = index(i)
514 kini(jlt_new) = kini(i)
515 vxi(jlt_new) = vxi(i)
516 vyi(jlt_new) = vyi(i)
517 vzi(jlt_new) = vzi(i)
518 msi(jlt_new) = msi(i)
519 ENDIF
520 ENDDO
521
522 RETURN