261
262
263
264#include "implicit_f.inc"
265
266
267
268 INTEGER IRECT(4,*),NSV(*),IRTL(*),IPARI(*),IRUPT(*)
269 my_real x(3,*),in(*),stifn(*),stfn(*),stfr(*),crst(2,*),visc
270
271
272
273 INTEGER NIR,I,J,II,JJ,L,W,NN,KK,
274 . IX1, IX2, IX3, IX4,NSVG,NSN
276 . s,t,sp,sm,tp,tm,e1x,e1y,e1z,e2x,e2y,e2z,e3x,e3y,e3z,
277 . xsm,ysm,zsm,xm,ym,zm,x1,x2,x3,x4,y1,y2,y3,y4,z1,z2,z3,z4,x0,y0,z0,xs,ys,zs,stf_mom,
278 . stf,str,stbrk
280 . h(4),h2(4),rx(4),ry(4),rz(4),rm(3),rs(3),stif,vis
282 . len2,fac_triang,irot,skew(9),tt,derx,dery,derz,det,b1,b2,b3,c1,c2,c3,bid,bid3(3)
283
284 nsn = ipari(5)
285
286 bid = zero
287 bid3(1:3)=zero
288 tt = zero
289
290 DO ii=1,nsn
291 IF (irupt(ii) == 0) cycle
292 i = nsv(ii)
293 l = irtl(ii)
294
295 ix1 = irect(1,l)
296 ix2 = irect(2,l)
297 ix3 = irect(3,l)
298 ix4 = irect(4,l)
299
300 IF (i > 0) THEN
301 s = crst(1,ii)
302 t = crst(2,ii)
303 l = irtl(ii)
304
305 ix1 = irect(1,l)
306 ix2 = irect(2,l)
307 ix3 = irect(3,l)
308 ix4 = irect(4,l)
309
310 nir= 4
311 sp = one + s
312 sm = one - s
313 tp = fourth*(one + t)
314 tm = fourth*(one - t)
315
316 h(1)=fourth
317 h(2)=fourth
318 h(3)=fourth
319 h(4)=fourth
320
321 h2(1)=tm*sm
322 h2(2)=tm*sp
323 h2(3)=tp*sp
324 h2(4)=tp*sm
325
326 IF (ix3 == ix4) THEN
327 nir = 3
328 h(1)=third
329 h(2)=third
330 h(3)=third
331 h(4) = zero
332 ENDIF
333
334
335
336 x1 = x(1,ix1)
337 y1 = x(2,ix1)
338 z1 = x(3,ix1)
339 x2 = x(1,ix2)
340 y2 = x(2,ix2)
341 z2 = x(3,ix2)
342 x3 = x(1,ix3)
343 y3 = x(2,ix3)
344 z3 = x(3,ix3)
345 x4 = x(1,ix4)
346 y4 = x(2,ix4)
347 z4 = x(3,ix4)
348 xs = x(1,i)
349 ys = x(2,i)
350 zs = x(3,i)
351
352
353 CALL i2rep(x1 ,x2 ,x3 ,x4 ,
354 . y1 ,y2 ,y3 ,y4 ,
355 . z1 ,z2 ,z3 ,z4 ,
356 . e1x ,e1y ,e1z ,
357 . e2x ,e2y ,e2z ,
358 . e3x ,e3y ,e3z ,nir )
359
360 IF (nir == 4) THEN
361 fac_triang = one
362
363 xm = x1*h2(1) + x2*h2(2) + x3*h2(3) + x4*h2(4)
364 ym = y1*h2(1) + y2*h2(2) + y3*h2(3) + y4*h2(4)
365 zm = z1*h2(1) + z2*h2(2) + z3*h2(3) + z4*h2(4)
366 x0 = (x1 + x2 + x3 + x4)/nir
367 y0 = (y1 + y2 + y3 + y4)/nir
368 z0 = (z1 + z2 + z3 + z4)/nir
369
370 xm = xm - x0
371 ym = ym - y0
372 zm = zm - z0
373 xs = xs - x0
374 ys = ys - y0
375 zs = zs - z0
376 xsm = xs - xm
377 ysm = ys - ym
378 zsm = zs - zm
379
380 ELSE
381 x0 = (x1 + x2 + x3)/nir
382 y0 = (y1 + y2 + y3)/nir
383 z0 = (z1 + z2 + z3)/nir
384
385 xm = x1*h(1) + x2*h(2) + x3*h(3)
386 ym = y1*h(1) + y2*h(2) + y3*h(3)
387 zm = z1*h(1) + z2*h(2) + z3*h(3)
388
389 xm = xm - x0
390 ym = ym - y0
391 zm = zm - z0
392 xs = xs - x0
393 ys = ys - y0
394 zs = zs - z0
395 xsm = xs - xm
396 ysm = ys - ym
397 zsm = zs - zm
398 ENDIF
399
400 x1 = x1 - x0
401 y1 = y1 - y0
402 z1 = z1 - z0
403 x2 = x2 - x0
404 y2 = y2 - y0
405 z2 = z2 - z0
406 x3 = x3 - x0
407 y3 = y3 - y0
408 z3 = z3 - z0
409 x4 = x4 - x0
410 y4 = y4 - y0
411 z4 = z4 - z0
412
413
414
415 rs(1) = xs*e1x + ys*e1y + zs*e1z
416 rs(2) = xs*e2x + ys*e2y + zs*e2z
417 rs(3) = xs*e3x + ys*e3y + zs*e3z
418 rm(1) = xm*e1x + ym*e1y + zm*e1z
419 rm(2) = xm*e2x + ym*e2y + zm*e2z
420 rm(3) = xm*e3x + ym*e3y + zm*e3z
421
422 rx(1) = e1x*x1 + e1y*y1 + e1z*z1
423 ry(1) = e2x*x1 + e2y*y1 + e2z*z1
424 rz(1) = e3x*x1 + e3y*y1 + e3z*z1
425 rx(2) = e1x*x2 + e1y*y2 + e1z*z2
426 ry(2) = e2x*x2 + e2y*y2 + e2z*z2
427 rz(2) = e3x*x2 + e3y*y2 + e3z*z2
428 rx(3) = e1x*x3 + e1y*y3 + e1z*z3
429 ry(3) = e2x*x3 + e2y*y3 + e2z*z3
430 rz(3) = e3x*x3 + e3y*y3 + e3z*z3
431 rx(4) = e1x*x4 + e1y*y4 + e1z*z4
432 ry(4) = e2x*x4 + e2y*y4 + e2z*z4
433 rz(4) = e3x*x4 + e3y*y4 + e3z*z4
434
435 IF (nir==3) THEN
436 rx(4)=zero
437 ry(4)=zero
438 rz(4)=zero
439 END IF
440
441
443 . skew ,tt ,bid ,stbrk,
444 . rs ,rm ,bid3 ,bid3 ,bid3 ,
445 . rx ,ry ,rz ,bid3 ,bid3 ,
446 . bid3 ,bid3 ,bid3 ,bid3 ,det ,
447 . b1 ,b2 ,b3 ,c1 ,c2 ,
448 . c3 ,in(i))
449
450
451
452 stf = stfn(ii)*(visc + sqrt(visc**2 + (one+stbrk)))**2
453
454 len2 = xsm**2+ysm**2+zsm**2
455 str = (stfr(ii)+stfn(ii)*len2)*(visc + sqrt(visc**2 + one))**2
456
457 derx = (b1+c3+c2)
458 dery = (b2+c1+c3)
459 derz = (b3+c2+c1)
460
461 stf_mom = det*
max(derx,dery,derz)*(str+stf*(xm*xm+ym*ym+zm*zm))
462
463
464 stifn(ix1) = stifn(ix1)+abs(stf*h(1))+stf_mom
465 stifn(ix2) = stifn(ix2)+abs(stf*h(2))+stf_mom
466 stifn(ix3) = stifn(ix3)+abs(stf*h(3))+stf_mom
467 stifn(ix4) = stifn(ix4)+abs(stf*h(4))+stf_mom
468
469 END IF
470 ENDDO
471
472
473 RETURN
subroutine i2pen_rot28(skew, tt, dt1, stif, rs, rm, vx, vy, vz, rx, ry, rz, va, vb, vc, vd, vrm, vrs, det, b1, b2, b3, c1, c2, c3, in_secnd)
subroutine i2rep(x1, x2, x3, x4, y1, y2, y3, y4, z1, z2, z3, z4, e1x, e1y, e1z, e2x, e2y, e2z, e3x, e3y, e3z, nir)