342
343
344
345#include "implicit_f.inc"
346
347
348
349 INTEGER IRECT(4,*),NSV(*),IRTL(*),IPARI(*)
350 my_real x(3,*),stifn(*),stfn(*),crst(2,*),visc
351
352
353
354 INTEGER NIR,I,J,II,JJ,L,W,NN,KK,LLT,
355 . IX1, IX2, IX3, IX4,NSVG,NSN
357 . s,t,sp,sm,tp,tm,e1x,e1y,e1z,e2x,e2y,e2z,e3x,e3y,e3z,
358 . xsm,ysm,zsm,xm,ym,zm,x1,x2,x3,x4,y1,y2,y3,y4,z1,z2,z3,z4,x0,y0,z0,xs,ys,zs,stifm,
359 . stf,str,stbrk
361 . h(4),rx(4),ry(4),rz(4),rm(3),rs(3),stif, vis
363 . len2,fac_triang,irot,skew(9),tt,bid,bid3(3),bid4(4,3)
364
365 nsn = ipari(5)
366
367 bid = zero
368 bid3(1:3) = zero
369 bid4(1:4,1:3) = zero
370 tt = zero
371
372 DO ii=1,nsn
373
374 i = nsv(ii)
375 l = irtl(ii)
376
377 ix1 = irect(1,l)
378 ix2 = irect(2,l)
379 ix3 = irect(3,l)
380 ix4 = irect(4,l)
381
382 IF (i > 0) THEN
383 s = crst(1,ii)
384 t = crst(2,ii)
385 l = irtl(ii)
386
387 ix1 = irect(1,l)
388 ix2 = irect(2,l)
389 ix3 = irect(3,l)
390 ix4 = irect(4,l)
391
392 nir= 4
393 sp = one + s
394 sm = one - s
395 tp = fourth*(one + t)
396 tm = fourth*(one - t)
397
398 h(1)=tm*sm
399 h(2)=tm*sp
400 h(3)=tp*sp
401 h(4)=tp*sm
402
403 IF (ix3 == ix4) THEN
404 nir = 3
405 h(3) = h(3) + h(4)
406 h(4) = zero
407 ENDIF
408
409
410
411 x1 = x(1,ix1)
412 y1 = x(2,ix1)
413 z1 = x(3,ix1)
414 x2 = x(1,ix2)
415 y2 = x(2,ix2)
416 z2 = x(3,ix2)
417 x3 = x(1,ix3)
418 y3 = x(2,ix3)
419 z3 = x(3,ix3)
420 x4 = x(1,ix4)
421 y4 = x(2,ix4)
422 z4 = x(3,ix4)
423 xs = x(1,i)
424 ys = x(2,i)
425 zs = x(3,i)
426
427
428 CALL i2rep(x1 ,x2 ,x3 ,x4 ,
429 . y1 ,y2 ,y3 ,y4 ,
430 . z1 ,z2 ,z3 ,z4 ,
431 . e1x ,e1y ,e1z ,
432 . e2x ,e2y ,e2z ,
433 . e3x ,e3y ,e3z ,nir )
434
435 IF (nir == 4) THEN
436 fac_triang = one
437
438 xm = x1*h(1) + x2*h(2) + x3*h(3) + x4*h(4)
439 ym = y1*h(1) + y2*h(2) + y3*h(3) + y4*h(4)
440 zm = z1*h(1) + z2*h(2) + z3*h(3) + z4*h(4)
441 x0 = (x1 + x2 + x3 + x4)/nir
442 y0 = (y1 + y2 + y3 + y4)/nir
443 z0 = (z1 + z2 + z3 + z4)/nir
444
445 xm = xm - x0
446 ym = ym - y0
447 zm = zm - z0
448 xs = xs - x0
449 ys = ys - y0
450 zs = zs - z0
451 xsm = xs - xm
452 ysm = ys - ym
453 zsm = zs - zm
454
455 ELSE
456 fac_triang = zero
457
458 x0 = (x1 + x2 + x3)/nir
459 y0 = (y1 + y2 + y3)/nir
460 z0 = (z1 + z2 + z3)/nir
461
462 xm = x1*h(1) + x2*h(2) + x3*h(3)
463 ym = y1*h(1) + y2*h(2) + y3*h(3)
464 zm = z1*h(1) + z2*h(2) + z3*h(3)
465
466 xm = xm - x0
467 ym = ym - y0
468 zm = zm - z0
469 xs = xs - x0
470 ys = ys - y0
471 zs = zs - z0
472 xsm = xs - xm
473 ysm = ys - ym
474 zsm = zs - zm
475 ENDIF
476
477 x1 = x1 - x0
478 y1 = y1 - y0
479 z1 = z1 - z0
480 x2 = x2 - x0
481 y2 = y2 - y0
482 z2 = z2 - z0
483 x3 = x3 - x0
484 y3 = y3 - y0
485 z3 = z3 - z0
486 x4 = x4 - x0
487 y4 = y4 - y0
488 z4 = z4 - z0
489
490
491
492 rs(1) = xs*e1x + ys*e1y + zs*e1z
493 rs(2) = xs*e2x + ys*e2y + zs*e2z
494 rs(3) = xs*e3x + ys*e3y + zs*e3z
495 rm(1) = xm*e1x + ym*e1y + zm*e1z
496 rm(2) = xm*e2x + ym*e2y + zm*e2z
497 rm(3) = xm*e3x + ym*e3y + zm*e3z
498
499 rx(1) = e1x*x1 + e1y*y1 + e1z*z1
500 ry(1) = e2x*x1 + e2y*y1 + e2z*z1
501 rz(1) = e3x*x1 + e3y*y1 + e3z*z1
502 rx(2) = e1x*x2 + e1y*y2 + e1z*z2
503 ry(2) = e2x*x2 + e2y*y2 + e2z*z2
504 rz(2) = e3x*x2 + e3y*y2 + e3z*z2
505 rx(3) = e1x*x3 + e1y*y3 + e1z*z3
506 ry(3) = e2x*x3 + e2y*y3 + e2z*z3
507 rz(3) = e3x*x3 + e3y*y3 + e3z*z3
508 rx(4) = e1x*x4 + e1y*y4 + e1z*z4
509 ry(4) = e2x*x4 + e2y*y4 + e2z*z4
510 rz(4) = e3x*x4 + e3y*y4 + e3z*z4
511
512 IF (nir==3) THEN
513 rx(4)=zero
514 ry(4)=zero
515 rz(4)=zero
516 END IF
517
518
520 . rs ,rm ,bid ,bid ,bid ,
521 . rx ,ry ,rz ,bid3 ,bid3 ,
522 . bid3 ,bid3)
523
524
525 stf = stfn(ii)*(visc + sqrt(visc**2 + (one+stbrk)))**2
526 stifm=zero
527
528
529
530
531
532 str = zero
533
534 CALL i2loceq( nir ,rs ,rx ,ry ,rz ,
535 . bid4(1:4,1) ,bid4(1:4,2) ,bid4(1:4,3) ,h(1) ,stifm)
536
537
538
539 stifn(ix1) = stifn(ix1)+abs(stf*h(1))+stifm*stf
540 stifn(ix2) = stifn(ix2)+abs(stf*h(2))+stifm*stf
541 stifn(ix3) = stifn(ix3)+abs(stf*h(3))+stifm*stf
542 stifn(ix4) = stifn(ix4)+abs(stf*h(4))+stifm*stf*fac_triang
543
544 END IF
545
546 ENDDO
547
548
549 RETURN
subroutine i2loceq(nir, rs, rx, ry, rz, fmx, fmy, fmz, h, stifm)
subroutine i2pen_rot(skew, tt, dt1, stif, rs, rm, v1, v2, v3, rx, ry, rz, va, vb, vc, vd)
subroutine i2rep(x1, x2, x3, x4, y1, y2, y3, y4, z1, z2, z3, z4, e1x, e1y, e1z, e2x, e2y, e2z, e3x, e3y, e3z, nir)