OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
i24for3.F File Reference
#include "implicit_f.inc"
#include "comlock.inc"
#include "mvsiz_p.inc"
#include "com01_c.inc"
#include "com04_c.inc"
#include "com06_c.inc"
#include "com08_c.inc"
#include "scr05_c.inc"
#include "scr07_c.inc"
#include "scr11_c.inc"
#include "scr14_c.inc"
#include "scr16_c.inc"
#include "scr18_c.inc"
#include "sms_c.inc"
#include "parit_c.inc"
#include "param_c.inc"
#include "impl1_c.inc"
#include "lockon.inc"
#include "lockoff.inc"
#include "task_c.inc"

Go to the source code of this file.

Functions/Subroutines

subroutine i24for3 (jlt, a, v, ibcc, icodt, fsav, gap, fric, ms, visc, viscf, noint, stfn, itab, cn_loc, stiglo, stifn, stif, fskyi, isky, n1, n2, n3, h1, h2, h3, h4, fcont, pene, ix1, ix2, ix3, ix4, nsvg, ivis2, neltst, ityptst, dt2t, subtria, gapv, inacti, index, niskyfi, kinet, newfront, isecin, nstrf, secfcum, x, irect, ce_loc, mfrot, ifq, frot_p, secnd_fr, alpha0, ibag, icontact, irtlm, viscn, vxi, vyi, vzi, msi, kini, nin, nisub, lisub, addsubs, addsubm, lisubs, lisubm, fsavsub, cand_n, ilagm, icurv, fncont, ftcont, nsn, x1, x2, x3, x4, y1, y2, y3, y4, z1, z2, z3, z4, xi, yi, zi, iadm, rcurvi, rcontact, acontact, pcontact, anglmi, padm, intth, phi, fthe, ftheskyi, temp, tempi, rstif, iform, mskyi_sms, iskyi_sms, nsms, cand_n_n, pene_old, stif_old, mbinflg, ilev, igsti, kmin, intply, iply, inod_pxfem, nm1, nm2, nm3, nrebou, irtse, nsne, is2se, is2pt, msegtyp, jtask, isensint, fsavparit, nft, h3d_data, fricc, viscffric, fric_coefs, t2main_sms, intnitsche, forneqsi, iorthfric, fric_coefs2, fricc2, viscffric2, nforth, nfisot, indexorth, indexisot, dir1, dir2, t2fac_sms, f_pfit, tagncont, kloadpinter, loadpinter, loadp_hyd_inter, typsub, inflg_subs, inflg_subm, ninloadp, dgaploadint, s_loadpinter, dist, ixx, interefric, intcarea, parameters, penref, kmax, s_addsubm, s_lisubm, s_typsub, nisubmax, i_stok, nrtm, nrtse, nsnr)
subroutine i24ass2 (jlt, ix1, ix2, ix3, ix4, nsvg, h1, h2, h3, h4, stif, fx1, fy1, fz1, fx2, fy2, fz2, fx3, fy3, fz3, fx4, fy4, fz4, fxi, fyi, fzi, fskyi, isky, niskyfi, nin, noint, intth, phi, ftheskyi, phi1, phi2, phi3, phi4, itab, intply, iply, inod, irtse, nsne, is2se, is2pt, tagip)
subroutine i24ass0 (jlt, ix1, ix2, ix3, ix4, nsvg, h1, h2, h3, h4, stif, fx1, fy1, fz1, fx2, fy2, fz2, fx3, fy3, fz3, fx4, fy4, fz4, fxi, fyi, fzi, a, stifn, nin, intth, phi, fthe, phi1, phi2, phi3, phi4, intply, iply, inod, irtse, nsne, is2se, is2pt, jtask)
subroutine i24sms2 (jlt, ix1, ix2, ix3, ix4, nsvg, h1, h2, h3, h4, stif, nin, noint, mskyi_sms, iskyi_sms, nsms, kt, c, cf, dtmini, dti, irtse, nsne, is2se, is2pt, t2main_sms, t2fac_sms)

Function/Subroutine Documentation

◆ i24ass0()

subroutine i24ass0 ( integer jlt,
integer, dimension(mvsiz) ix1,
integer, dimension(mvsiz) ix2,
integer, dimension(mvsiz) ix3,
integer, dimension(mvsiz) ix4,
integer, dimension(mvsiz) nsvg,
h1,
h2,
h3,
h4,
stif,
fx1,
fy1,
fz1,
fx2,
fy2,
fz2,
fx3,
fy3,
fz3,
fx4,
fy4,
fz4,
fxi,
fyi,
fzi,
a,
stifn,
integer nin,
integer intth,
phi,
fthe,
phi1,
phi2,
phi3,
phi4,
integer intply,
integer, dimension(4,*) iply,
integer, dimension(*) inod,
integer, dimension(5,*) irtse,
integer nsne,
integer, dimension(2,*) is2se,
integer, dimension(*) is2pt,
integer jtask )

Definition at line 4426 of file i24for3.F.

4434C-----------------------------------------------
4435C M o d u l e s
4436C-----------------------------------------------
4437 USE tri7box
4438 USE plyxfem_mod
4439C-----------------------------------------------
4440C I m p l i c i t T y p e s
4441C-----------------------------------------------
4442#include "implicit_f.inc"
4443C-----------------------------------------------
4444C G l o b a l P a r a m e t e r s
4445C-----------------------------------------------
4446#include "mvsiz_p.inc"
4447C-----------------------------------------------
4448C C o m m o n B l o c k s
4449C-----------------------------------------------
4450#include "com04_c.inc"
4451C-----------------------------------------------
4452C D u m m y A r g u m e n t s
4453C-----------------------------------------------
4454 INTEGER JLT, NIN,INTTH,
4455 . IX1(MVSIZ),IX2(MVSIZ),IX3(MVSIZ),IX4(MVSIZ),NSVG(MVSIZ),
4456 . INTPLY,IPLY(4,*),INOD(*),IRTSE(5,*),NSNE,IS2SE(2,*),
4457 . IS2PT(*) ,JTASK
4458 my_real
4459 . h1(mvsiz),h2(mvsiz),h3(mvsiz),h4(mvsiz),stif(mvsiz),
4460 . fx1(mvsiz),fy1(mvsiz),fz1(mvsiz),
4461 . fx2(mvsiz),fy2(mvsiz),fz2(mvsiz),
4462 . fx3(mvsiz),fy3(mvsiz),fz3(mvsiz),
4463 . fx4(mvsiz),fy4(mvsiz),fz4(mvsiz),
4464 . fxi(mvsiz),fyi(mvsiz),fzi(mvsiz),
4465 . a(3,*), stifn(*),phi(*), fthe(*),
4466 . phi1(*), phi2(*), phi3(*), phi4(*)
4467C-----------------------------------------------
4468C L o c a l V a r i a b l e s
4469C-----------------------------------------------
4470 my_real
4471 . hh(mvsiz),fici(5),fics(4,4),ficsth(4,5)
4472 INTEGER I, J1, IG,ILY,NN,JG,IXSS(4),NFIC,J,ISHIFT,NODFI
4473
4474C
4475 IF(intth == 0) THEN
4476 DO i=1,jlt
4477 hh(i)=h1(i)+h2(i)+h3(i)+h4(i)
4478 IF(hh(i)==zero)cycle
4479 j1=ix1(i)
4480 a(1,j1)=a(1,j1)+fx1(i)
4481 a(2,j1)=a(2,j1)+fy1(i)
4482 a(3,j1)=a(3,j1)+fz1(i)
4483 stifn(j1) = stifn(j1) + stif(i)*abs(h1(i))
4484C
4485 j1=ix2(i)
4486 a(1,j1)=a(1,j1)+fx2(i)
4487 a(2,j1)=a(2,j1)+fy2(i)
4488 a(3,j1)=a(3,j1)+fz2(i)
4489 stifn(j1) = stifn(j1) + stif(i)*abs(h2(i))
4490C
4491 j1=ix3(i)
4492 a(1,j1)=a(1,j1)+fx3(i)
4493 a(2,j1)=a(2,j1)+fy3(i)
4494 a(3,j1)=a(3,j1)+fz3(i)
4495 stifn(j1) = stifn(j1) + stif(i)*abs(h3(i))
4496C
4497 j1=ix4(i)
4498 a(1,j1)=a(1,j1)+fx4(i)
4499 a(2,j1)=a(2,j1)+fy4(i)
4500 a(3,j1)=a(3,j1)+fz4(i)
4501 stifn(j1) = stifn(j1) + stif(i)*abs(h4(i))
4502 ENDDO
4503 ELSE
4504 DO i=1,jlt
4505 hh(i)=h1(i)+h2(i)+h3(i)+h4(i)
4506 IF(hh(i)==zero)cycle
4507 j1=ix1(i)
4508 a(1,j1)=a(1,j1)+fx1(i)
4509 a(2,j1)=a(2,j1)+fy1(i)
4510 a(3,j1)=a(3,j1)+fz1(i)
4511 stifn(j1) = stifn(j1) + stif(i)*abs(h1(i))
4512 fthe(j1) = fthe(j1) + phi1(i)
4513C
4514 j1=ix2(i)
4515 a(1,j1)=a(1,j1)+fx2(i)
4516 a(2,j1)=a(2,j1)+fy2(i)
4517 a(3,j1)=a(3,j1)+fz2(i)
4518 stifn(j1) = stifn(j1) + stif(i)*abs(h2(i))
4519 fthe(j1) = fthe(j1) + phi2(i)
4520C
4521 j1=ix3(i)
4522 a(1,j1)=a(1,j1)+fx3(i)
4523 a(2,j1)=a(2,j1)+fy3(i)
4524 a(3,j1)=a(3,j1)+fz3(i)
4525 stifn(j1) = stifn(j1) + stif(i)*abs(h3(i))
4526 fthe(j1) = fthe(j1) + phi3(i)
4527C
4528 j1=ix4(i)
4529 a(1,j1)=a(1,j1)+fx4(i)
4530 a(2,j1)=a(2,j1)+fy4(i)
4531 a(3,j1)=a(3,j1)+fz4(i)
4532 stifn(j1) = stifn(j1) + stif(i)*abs(h4(i))
4533 fthe(j1) = fthe(j1) + phi4(i)
4534 ENDDO
4535 ENDIF
4536 IF(intply > 0) THEN
4537 DO i=1,jlt
4538 hh(i)=h1(i)+h2(i)+h3(i)+h4(i)
4539 IF(hh(i)==zero)cycle
4540 j1=ix1(i)
4541 nn = inod(j1)
4542 ily = iply(1,i)
4543 IF(inod(j1) > 0 .AND. ily > 0 ) THEN
4544 ply(ily)%A(1,nn)=ply(ily)%A(1,nn) + fx1(i)
4545 ply(ily)%A(2,nn)=ply(ily)%A(2,nn) + fy1(i)
4546 ply(ily)%A(3,nn)=ply(ily)%A(3,nn) + fz1(i)
4547 ply(ily)%A(4,nn)=ply(ily)%A(4,nn) + stif(i)*abs(h1(i))
4548 ENDIF
4549C
4550 j1=ix2(i)
4551 nn = inod(j1)
4552 ily = iply(2,i)
4553 IF(inod(j1) > 0 .AND. ily > 0) THEN
4554 ply(ily)%A(1,nn)=ply(ily)%A(1,nn) +fx2(i)
4555 ply(ily)%A(2,nn)=ply(ily)%A(2,nn) +fy2(i)
4556 ply(ily)%A(3,nn)=ply(ily)%A(3,nn) +fz2(i)
4557 ply(ily)%A(4,nn)=ply(ily)%A(4,nn) + stif(i)*abs(h2(i))
4558 ENDIF
4559C
4560 j1=ix3(i)
4561 nn = inod(j1)
4562 ily = iply(3,i)
4563 IF(inod(j1) > 0 .AND. ily > 0) THEN
4564 ply(ily)%A(1,nn)=ply(ily)%A(1,nn) +fx3(i)
4565 ply(ily)%A(2,nn)=ply(ily)%A(2,nn) +fy3(i)
4566 ply(ily)%A(3,nn)=ply(ily)%A(3,nn) +fz3(i)
4567 ply(ily)%A(4,nn)=ply(ily)%A(4,nn) + stif(i)*abs(h3(i))
4568 ENDIF
4569C
4570 j1=ix4(i)
4571 nn = inod(j1)
4572 ily = iply(4,i)
4573 IF(inod(j1) > 0 .AND. ily > 0) THEN
4574 ply(ily)%A(1,nn)=ply(ily)%A(1,nn) +fx4(i)
4575 ply(ily)%A(2,nn)=ply(ily)%A(2,nn) +fy4(i)
4576 ply(ily)%A(3,nn)=ply(ily)%A(3,nn) +fz4(i)
4577 ply(ily)%A(4,nn)=ply(ily)%A(4,nn) + stif(i)*abs(h4(i))
4578 ENDIF
4579 ENDDO
4580 ENDIF
4581C
4582 nodfi = nlskyfi(nin)
4583 ishift = nodfi*(jtask-1)
4584 IF(intth == 0 ) THEN
4585 DO i=1,jlt
4586 IF(hh(i)==zero)cycle
4587 ig=nsvg(i)
4588 IF(ig>0)THEN
4589 IF (ig >numnod) THEN
4590 jg = ig - numnod
4591 nfic=1
4592 fici(nfic)=fxi(i)
4593 nfic = nfic + 1
4594 fici(nfic)=fyi(i)
4595 nfic = nfic + 1
4596 fici(nfic)=fzi(i)
4597 nfic = nfic + 1
4598 fici(nfic)=stif(i)
4599 CALL i24forc_fic( 3,irtse ,nsne ,is2se ,is2pt ,
4600 + jg ,nfic ,fici ,fics ,ixss )
4601 DO j = 1,4
4602 j1 = ixss(j)
4603 a(1,j1)=a(1,j1)-fics(j,1)
4604 a(2,j1)=a(2,j1)-fics(j,2)
4605 a(3,j1)=a(3,j1)-fics(j,3)
4606 stifn(j1) = stifn(j1) + fics(j,4)
4607 END DO
4608 ELSE
4609 a(1,ig)=a(1,ig)-fxi(i)
4610 a(2,ig)=a(2,ig)-fyi(i)
4611 a(3,ig)=a(3,ig)-fzi(i)
4612 stifn(ig) = stifn(ig) + stif(i)
4613 END IF !(IG >NUMNOD) THEN
4614 ELSE
4615 ig = -ig
4616 IF(isedge_fi(nin)%P(ig)==1)THEN
4617 jg = ig
4618 nfic=1
4619 fici(nfic)=fxi(i)
4620 nfic = nfic + 1
4621 fici(nfic)=fyi(i)
4622 nfic = nfic + 1
4623 fici(nfic)=fzi(i)
4624 nfic = nfic + 1
4625 fici(nfic)=stif(i)
4626 CALL i24forc_fic( 3,irtse_fi(nin)%P(1,1) ,nsne ,is2se_fi(nin)%P(1,1) ,is2pt_fi(nin)%P(1) ,
4627 + jg ,nfic ,fici ,fics ,ixss )
4628 DO j = 1,4
4629 j1 = ixss(j)
4630 afi(nin)%P(1,j1+ishift)=afi(nin)%P(1,j1+ishift)-fics(j,1)
4631 afi(nin)%P(2,j1+ishift)=afi(nin)%P(2,j1+ishift)-fics(j,2)
4632 afi(nin)%P(3,j1+ishift)=afi(nin)%P(3,j1+ishift)-fics(j,3)
4633 stnfi(nin)%P(j1+ishift)=stnfi(nin)%P(j1+ishift) + fics(j,4)
4634 END DO
4635
4636 ELSE
4637 afi(nin)%P(1,ig+ishift)=afi(nin)%P(1,ig+ishift)-fxi(i)
4638 afi(nin)%P(2,ig+ishift)=afi(nin)%P(2,ig+ishift)-fyi(i)
4639 afi(nin)%P(3,ig+ishift)=afi(nin)%P(3,ig+ishift)-fzi(i)
4640 stnfi(nin)%P(ig+ishift)=stnfi(nin)%P(ig+ishift) + stif(i)
4641 ENDIF
4642 ENDIF
4643 ENDDO
4644C
4645 ELSE
4646 DO i=1,jlt
4647 IF(hh(i)==zero)cycle
4648 ig=nsvg(i)
4649 IF(ig>0)THEN
4650 IF (ig >numnod) THEN
4651 jg = ig - numnod
4652 nfic=1
4653 fici(nfic)=fxi(i)
4654 nfic = nfic + 1
4655 fici(nfic)=fyi(i)
4656 nfic = nfic + 1
4657 fici(nfic)=fzi(i)
4658 nfic = nfic + 1
4659 fici(nfic)=stif(i)
4660 nfic = nfic + 1
4661 fici(nfic)=phi(i)
4662 CALL i24forc_fic(3 ,irtse ,nsne ,is2se ,is2pt ,
4663 + jg ,nfic ,fici ,ficsth ,ixss )
4664 DO j = 1,4
4665 j1 = ixss(j)
4666 a(1,j1)=a(1,j1)-ficsth(j,1)
4667 a(2,j1)=a(2,j1)-ficsth(j,2)
4668 a(3,j1)=a(3,j1)-ficsth(j,3)
4669 stifn(j1) = stifn(j1) + ficsth(j,4)
4670 fthe(j1)=fthe(j1) + ficsth(j,5)
4671 END DO
4672 ELSE
4673 a(1,ig)=a(1,ig)-fxi(i)
4674 a(2,ig)=a(2,ig)-fyi(i)
4675 a(3,ig)=a(3,ig)-fzi(i)
4676 stifn(ig) = stifn(ig) + stif(i)
4677 fthe(ig)=fthe(ig) + phi(i)
4678 END IF !(IG >NUMNOD) THEN
4679 ELSE
4680 ig = -ig
4681 IF(isedge_fi(nin)%P(ig)==1)THEN
4682 jg = ig
4683 nfic=1
4684 fici(nfic)=fxi(i)
4685 nfic = nfic + 1
4686 fici(nfic)=fyi(i)
4687 nfic = nfic + 1
4688 fici(nfic)=fzi(i)
4689 nfic = nfic + 1
4690 fici(nfic)=stif(i)
4691 nfic = nfic + 1
4692 fici(nfic)=phi(i)
4693 CALL i24forc_fic(3 ,irtse_fi(nin)%P(1,1) ,nsne ,is2se_fi(nin)%P(1,1) ,is2pt_fi(nin)%P(1) ,
4694 + jg ,nfic ,fici ,ficsth ,ixss )
4695 DO j = 1,4
4696 j1 = ixss(j)
4697 afi(nin)%P(1,j1+ishift)=afi(nin)%P(1,j1+ishift)-ficsth(j,1)
4698 afi(nin)%P(2,j1+ishift)=afi(nin)%P(2,j1+ishift)-ficsth(j,2)
4699 afi(nin)%P(3,j1+ishift)=afi(nin)%P(3,j1+ishift)-ficsth(j,3)
4700
4701 stnfi(nin)%P(j1+ishift)=stnfi(nin)%P(j1+ishift) + ficsth(j,4)
4702 fthefi(nin)%P(j1+ishift)= fthefi(nin)%P(j1+ishift) + ficsth(j,5)
4703 END DO
4704
4705 ELSE
4706 afi(nin)%P(1,ig+ishift)=afi(nin)%P(1,ig+ishift)-fxi(i)
4707 afi(nin)%P(2,ig+ishift)=afi(nin)%P(2,ig+ishift)-fyi(i)
4708 afi(nin)%P(3,ig+ishift)=afi(nin)%P(3,ig+ishift)-fzi(i)
4709 stnfi(nin)%P(ig+ishift)=stnfi(nin)%P(ig+ishift) + stif(i)
4710 fthefi(nin)%P(ig+ishift)= fthefi(nin)%P(ig+ishift) + phi(i)
4711 ENDIF
4712 ENDIF
4713 ENDDO
4714 ENDIF
4715
4716cc IF(INTPLY > 0) THEN
4717cc DO I=1,JLT
4718cc IF(HH(I)==ZERO)CYCLE
4719cc IG=NSVG(I)
4720cc NN = INOD(IG)
4721cc IF(IG>0 .AND. NN > 0)THEN
4722cc ILY = IPLY(5,I)
4723cc PLY(ILY)%A(1,NN) = PLY(ILY)%A(1,NN) - FXI(I)
4724cc PLY(ILY)%A(2,NN) = PLY(ILY)%A(2,NN) - FYI(I)
4725cc PLY(ILY)%A(3,NN) = PLY(ILY)%A(3,NN) - FZI(I)
4726cc PLY(ILY)%A(4,NN) = PLY(ILY)%A(4,NN) + STIF(I)
4727cc ELSE
4728C avoir l'aspect spmd
4729cc IG = -IG
4730cc AFI(NIN)%P(1,IG)=AFI(NIN)%P(1,IG)-FXI(I)
4731cc AFI(NIN)%P(2,IG)=AFI(NIN)%P(2,IG)-FYI(I)
4732cc AFI(NIN)%P(3,IG)=AFI(NIN)%P(3,IG)-FZI(I)
4733cc ENDIF
4734cc ENDDO
4735cc ENDIF
4736C
4737 RETURN
my_real
#define my_real
Definition cppsort.cpp:32
i24forc_fic
subroutine i24forc_fic(npt, irtse, nsne, is2se, is2pt, ns, nfic, fici, fics, ixss)
Definition i24for3e.F:39
plyxfem_mod
Definition plyxfem_mod.F:31
plyxfem_mod::ply
type(ply_data), dimension(:), allocatable ply
Definition plyxfem_mod.F:91
tri7box::is2pt_fi
type(int_pointer), dimension(:), allocatable is2pt_fi
Definition tri7box.F:537
tri7box::is2se_fi
type(int_pointer2), dimension(:), allocatable is2se_fi
Definition tri7box.F:536
tri7box::stnfi
type(real_pointer), dimension(:), allocatable stnfi
Definition tri7box.F:449
tri7box::irtse_fi
type(int_pointer2), dimension(:), allocatable irtse_fi
Definition tri7box.F:535
tri7box::afi
type(real_pointer2), dimension(:), allocatable afi
Definition tri7box.F:459
tri7box::isedge_fi
type(int_pointer), dimension(:), allocatable isedge_fi
Definition tri7box.F:540
tri7box::nlskyfi
integer, dimension(:), allocatable nlskyfi
Definition tri7box.F:512
tri7box::fthefi
type(real_pointer), dimension(:), allocatable fthefi
Definition tri7box.F:449

◆ i24ass2()

subroutine i24ass2 ( integer jlt,
integer, dimension(mvsiz) ix1,
integer, dimension(mvsiz) ix2,
integer, dimension(mvsiz) ix3,
integer, dimension(mvsiz) ix4,
integer, dimension(mvsiz) nsvg,
h1,
h2,
h3,
h4,
stif,
fx1,
fy1,
fz1,
fx2,
fy2,
fz2,
fx3,
fy3,
fz3,
fx4,
fy4,
fz4,
fxi,
fyi,
fzi,
fskyi,
integer, dimension(*) isky,
integer niskyfi,
integer nin,
integer noint,
integer intth,
phi,
ftheskyi,
phi1,
phi2,
phi3,
phi4,
integer, dimension(*) itab,
integer intply,
integer, dimension(4,*) iply,
integer, dimension(*) inod,
integer, dimension(5,*) irtse,
integer nsne,
integer, dimension(2,*) is2se,
integer, dimension(*) is2pt,
integer, dimension(mvsiz), intent(inout) tagip )

Definition at line 3792 of file i24for3.F.

3800C-----------------------------------------------
3801C M o d u l e s
3802C-----------------------------------------------
3803 USE tri7box
3804 USE message_mod
3805 USE plyxfem_mod
3806 USE debug_mod
3807C-----------------------------------------------
3808C I m p l i c i t T y p e s
3809C-----------------------------------------------
3810#include "implicit_f.inc"
3811#include "comlock.inc"
3812C-----------------------------------------------
3813C G l o b a l P a r a m e t e r s
3814C-----------------------------------------------
3815#include "mvsiz_p.inc"
3816C-----------------------------------------------
3817C C o m m o n B l o c k s
3818C-----------------------------------------------
3819#include "parit_c.inc"
3820#include "com04_c.inc"
3821C-----------------------------------------------
3822C D u m m y A r g u m e n t s
3823C-----------------------------------------------
3824 INTEGER JLT,NISKYFI,NIN,NOINT,INTTH,
3825 . ISKY(*),ITAB(*),
3826 . IX1(MVSIZ),IX2(MVSIZ),IX3(MVSIZ),IX4(MVSIZ),NSVG(MVSIZ),
3827 . INTPLY,IPLY(4,*),INOD(*),
3828 . IRTSE(5,*),NSNE ,IS2SE(2,*),IS2PT(*)
3829 INTEGER , INTENT(INOUT) :: TAGIP(MVSIZ)
3830 my_real
3831 . h1(mvsiz),h2(mvsiz),h3(mvsiz),h4(mvsiz),stif(mvsiz),
3832 . fx1(mvsiz),fy1(mvsiz),fz1(mvsiz),
3833 . fx2(mvsiz),fy2(mvsiz),fz2(mvsiz),
3834 . fx3(mvsiz),fy3(mvsiz),fz3(mvsiz),
3835 . fx4(mvsiz),fy4(mvsiz),fz4(mvsiz),
3836 . fxi(mvsiz),fyi(mvsiz),fzi(mvsiz),
3837 . fskyi(lskyi,nfskyi),ftheskyi(lskyi),phi(mvsiz),
3838 . phi1(*),phi2(*) ,phi3(*) ,phi4(*)
3839C-----------------------------------------------
3840C L o c a l V a r i a b l e s
3841C-----------------------------------------------
3842 my_real
3843 . hh(mvsiz),fici(5),fics(4,4),ficsth(4,5)
3844 INTEGER I, J1, IG, NISKYL1, NISKYL,IGP,IGM,IDR,NISKYFIL,J
3845 INTEGER ITG,NN,ILY,JG,IXSS(4),NFIC
3846 INTEGER se
3847C
3848 niskyl1 = 0
3849 tagip(1:jlt) = 0
3850 DO i = 1, jlt
3851 hh(i)=h1(i)+h2(i)+h3(i)+h4(i)
3852 IF(hh(i)==zero.AND.tagip(i)==0)cycle
3853 IF (h1(i)/=zero.OR.tagip(i)==1) niskyl1 = niskyl1 + 1
3854 IF (h2(i)/=zero.OR.tagip(i)==1) niskyl1 = niskyl1 + 1
3855 IF (h3(i)/=zero.OR.tagip(i)==1) niskyl1 = niskyl1 + 1
3856 IF (h4(i)/=zero.OR.tagip(i)==1) niskyl1 = niskyl1 + 1
3857 ENDDO
3858C
3859C Precalcul impact locaux / remote
3860C
3861 igp = 0
3862 igm = 0
3863 DO i=1,jlt
3864 IF(hh(i)==zero.AND.tagip(i)==0)cycle
3865 ig =nsvg(i)
3866 IF(ig>0) THEN
3867 igp = igp+1
3868 IF (ig > numnod) igp = igp+3
3869 ELSE
3870 igm = igm+1
3871 IF(isedge_fi(nin)%P(-ig)==1)igm = igm+3
3872 ENDIF
3873 ENDDO
3874C
3875#include "lockon.inc"
3876 niskyl = nisky
3877 nisky = nisky + niskyl1 + igp
3878 niskyfil = niskyfi
3879 niskyfi = niskyfi + igm
3880#include "lockoff.inc"
3881C
3882 IF (niskyl+niskyl1+igp > lskyi) THEN
3883 CALL ancmsg(msgid=26,anmode=aninfo)
3884 CALL arret(2)
3885 ENDIF
3886 IF (niskyfil+igm > nlskyfi(nin)) THEN
3887 CALL ancmsg(msgid=26,anmode=aninfo)
3888 CALL arret(2)
3889 ENDIF
3890C
3891 IF(intth == 0 ) THEN
3892 IF(intply == 0 )THEN
3893 DO i=1,jlt
3894 IF (h1(i)/=zero.OR.tagip(i)==1) THEN
3895 niskyl = niskyl + 1
3896 fskyi(niskyl,1)=fx1(i)
3897 fskyi(niskyl,2)=fy1(i)
3898 fskyi(niskyl,3)=fz1(i)
3899 fskyi(niskyl,4)=stif(i)*abs(h1(i))
3900 isky(niskyl) = ix1(i)
3901 ENDIF
3902 ENDDO
3903 DO i=1,jlt
3904 IF (h2(i)/=zero.OR.tagip(i)==1) THEN
3905 niskyl = niskyl + 1
3906 fskyi(niskyl,1)=fx2(i)
3907 fskyi(niskyl,2)=fy2(i)
3908 fskyi(niskyl,3)=fz2(i)
3909 fskyi(niskyl,4)=stif(i)*abs(h2(i))
3910 isky(niskyl) = ix2(i)
3911 ENDIF
3912 ENDDO
3913 DO i=1,jlt
3914 IF (h3(i)/=zero.OR.tagip(i)==1) THEN
3915 niskyl = niskyl + 1
3916 fskyi(niskyl,1)=fx3(i)
3917 fskyi(niskyl,2)=fy3(i)
3918 fskyi(niskyl,3)=fz3(i)
3919 fskyi(niskyl,4)=stif(i)*abs(h3(i))
3920 isky(niskyl) = ix3(i)
3921 ENDIF
3922 ENDDO
3923 DO i=1,jlt
3924 IF (h4(i)/=zero.OR.tagip(i)==1) THEN
3925 niskyl = niskyl + 1
3926 fskyi(niskyl,1)=fx4(i)
3927 fskyi(niskyl,2)=fy4(i)
3928 fskyi(niskyl,3)=fz4(i)
3929 fskyi(niskyl,4)=stif(i)*abs(h4(i))
3930 isky(niskyl) = ix4(i)
3931 ENDIF
3932 ENDDO
3933C
3934 DO i=1,jlt
3935 IF(hh(i)==zero.AND.tagip(i)==0)cycle
3936 ig =nsvg(i)
3937 IF(ig>0) THEN
3938 IF (ig >numnod) THEN
3939 jg = ig - numnod
3940 nfic=1
3941 fici(nfic)=fxi(i)
3942 nfic = nfic + 1
3943 fici(nfic)=fyi(i)
3944 nfic = nfic + 1
3945 fici(nfic)=fzi(i)
3946 nfic = nfic + 1
3947 fici(nfic)=stif(i)
3948 CALL i24forc_fic(3 ,irtse ,nsne ,is2se ,is2pt ,
3949 + jg ,nfic ,fici ,fics ,ixss )
3950 DO j = 1,4
3951 j1 = ixss(j)
3952 niskyl = niskyl + 1
3953 fskyi(niskyl,1)=-fics(j,1)
3954 fskyi(niskyl,2)=-fics(j,2)
3955 fskyi(niskyl,3)=-fics(j,3)
3956 fskyi(niskyl,4)= fics(j,4)
3957 isky(niskyl) = j1
3958 END DO
3959 ELSE
3960 niskyl = niskyl + 1
3961 fskyi(niskyl,1)=-fxi(i)
3962 fskyi(niskyl,2)=-fyi(i)
3963 fskyi(niskyl,3)=-fzi(i)
3964 fskyi(niskyl,4)= stif(i)
3965 isky(niskyl) = ig
3966 END IF !(IG >NUMNOD) THEN
3967 ELSE
3968 ig = -ig
3969 IF(isedge_fi(nin)%P(ig)==1)THEN
3970 nfic=1
3971 fici(nfic)=fxi(i)
3972 nfic = nfic + 1
3973 fici(nfic)=fyi(i)
3974 nfic = nfic + 1
3975 fici(nfic)=fzi(i)
3976 nfic = nfic + 1
3977 fici(nfic)=stif(i)
3978 CALL i24forc_fic(3 ,irtse_fi(nin)%P(1,1) ,nsne ,is2se_fi(nin)%P(1,1) , is2pt_fi(nin)%P(1) ,
3979 + ig ,nfic ,fici ,fics ,ixss )
3980
3981 DO j = 1,4
3982 j1 = ixss(j)
3983 niskyfil = niskyfil + 1
3984 fskyfi(nin)%P(1,niskyfil)=-fics(j,1)
3985 fskyfi(nin)%P(2,niskyfil)=-fics(j,2)
3986 fskyfi(nin)%P(3,niskyfil)=-fics(j,3)
3987 fskyfi(nin)%P(4,niskyfil)= fics(j,4)
3988 iskyfi(nin)%P(niskyfil) = j1
3989 END DO
3990
3991 ELSE
3992 niskyfil = niskyfil + 1
3993 fskyfi(nin)%P(1,niskyfil)=-fxi(i)
3994 fskyfi(nin)%P(2,niskyfil)=-fyi(i)
3995 fskyfi(nin)%P(3,niskyfil)=-fzi(i)
3996 fskyfi(nin)%P(4,niskyfil)= stif(i)
3997 iskyfi(nin)%P(niskyfil) = ig
3998 ENDIF
3999 ENDIF
4000 ENDDO
4001 ELSE
4002C with plyxfem formulation s
4003 DO i=1,jlt
4004 IF (h1(i)/=zero.OR.tagip(i)==1) THEN
4005 niskyl = niskyl + 1
4006 fskyi(niskyl,1)=fx1(i)
4007 fskyi(niskyl,2)=fy1(i)
4008 fskyi(niskyl,3)=fz1(i)
4009 fskyi(niskyl,4)=stif(i)*abs(h1(i))
4010 isky(niskyl) = ix1(i)
4011C
4012 plyskyi%FSKYI(niskyl,1)=fx1(i)
4013 plyskyi%FSKYI(niskyl,2)=fy1(i)
4014 plyskyi%FSKYI(niskyl,3)=fz1(i)
4015 plyskyi%FSKYI(niskyl,4)=stif(i)*abs(h1(i))
4016 plyskyi%FSKYI(niskyl,5)=iply(1,i)
4017 ENDIF
4018 ENDDO
4019 DO i=1,jlt
4020 IF (h2(i)/=zero.OR.tagip(i)==1) THEN
4021 niskyl = niskyl + 1
4022 fskyi(niskyl,1)=fx2(i)
4023 fskyi(niskyl,2)=fy2(i)
4024 fskyi(niskyl,3)=fz2(i)
4025 fskyi(niskyl,4)=stif(i)*abs(h2(i))
4026 isky(niskyl) = ix2(i)
4027C
4028 plyskyi%FSKYI(niskyl,1)=fx2(i)
4029 plyskyi%FSKYI(niskyl,2)=fy2(i)
4030 plyskyi%FSKYI(niskyl,3)=fz2(i)
4031 plyskyi%FSKYI(niskyl,4)=stif(i)*abs(h2(i))
4032 plyskyi%FSKYI(niskyl,5)=iply(2,i)
4033 ENDIF
4034 ENDDO
4035 DO i=1,jlt
4036 IF (h3(i)/=zero.OR.tagip(i)==1) THEN
4037 niskyl = niskyl + 1
4038 fskyi(niskyl,1)=fx3(i)
4039 fskyi(niskyl,2)=fy3(i)
4040 fskyi(niskyl,3)=fz3(i)
4041 fskyi(niskyl,4)=stif(i)*abs(h3(i))
4042 isky(niskyl) = ix3(i)
4043C
4044 plyskyi%FSKYI(niskyl,1)=fx3(i)
4045 plyskyi%FSKYI(niskyl,2)=fy3(i)
4046 plyskyi%FSKYI(niskyl,3)=fz3(i)
4047 plyskyi%FSKYI(niskyl,4)=stif(i)*abs(h3(i))
4048 plyskyi%FSKYI(niskyl,5)=iply(3,i)
4049 ENDIF
4050 ENDDO
4051 DO i=1,jlt
4052 IF (h4(i)/=zero.OR.tagip(i)==1) THEN
4053 niskyl = niskyl + 1
4054 fskyi(niskyl,1)=fx4(i)
4055 fskyi(niskyl,2)=fy4(i)
4056 fskyi(niskyl,3)=fz4(i)
4057 fskyi(niskyl,4)=stif(i)*abs(h4(i))
4058 isky(niskyl) = ix4(i)
4059C
4060 plyskyi%FSKYI(niskyl,1)=fx4(i)
4061 plyskyi%FSKYI(niskyl,2)=fy4(i)
4062 plyskyi%FSKYI(niskyl,3)=fz4(i)
4063 plyskyi%FSKYI(niskyl,4)=stif(i)*abs(h4(i))
4064 plyskyi%FSKYI(niskyl,5)=iply(4,i)
4065 ENDIF
4066 ENDDO
4067C
4068 DO i=1,jlt
4069 IF(hh(i)==zero.AND.tagip(i)==0)cycle
4070 ig =nsvg(i)
4071 IF(ig>0) THEN
4072 IF (ig >numnod) THEN
4073 jg = ig - numnod
4074 nfic=1
4075 fici(nfic)=fxi(i)
4076 nfic = nfic + 1
4077 fici(nfic)=fyi(i)
4078 nfic = nfic + 1
4079 fici(nfic)=fzi(i)
4080 nfic = nfic + 1
4081 fici(nfic)=stif(i)
4082 CALL i24forc_fic(3 ,irtse ,nsne ,is2se ,is2pt ,
4083 + jg ,nfic ,fici ,fics ,ixss )
4084 DO j = 1,4
4085 j1 = ixss(j)
4086 niskyl = niskyl + 1
4087 fskyi(niskyl,1)=-fics(j,1)
4088 fskyi(niskyl,2)=-fics(j,2)
4089 fskyi(niskyl,3)=-fics(j,3)
4090 fskyi(niskyl,4)= fics(j,4)
4091 isky(niskyl) = j1
4092 END DO
4093 ELSE
4094 niskyl = niskyl + 1
4095 fskyi(niskyl,1)=-fxi(i)
4096 fskyi(niskyl,2)=-fyi(i)
4097 fskyi(niskyl,3)=-fzi(i)
4098 fskyi(niskyl,4)= stif(i)
4099 isky(niskyl) = ig
4100 END IF !(IG >NUMNOD) THEN
4101C
4102 ELSE
4103 ig = -ig
4104 IF(isedge_fi(nin)%P(ig)==1)THEN
4105 jg = ig
4106 nfic=1
4107 fici(nfic)=fxi(i)
4108 nfic = nfic + 1
4109 fici(nfic)=fyi(i)
4110 nfic = nfic + 1
4111 fici(nfic)=fzi(i)
4112 nfic = nfic + 1
4113 fici(nfic)=stif(i)
4114 CALL i24forc_fic(3 ,irtse_fi(nin)%P(1,1) ,nsne ,is2se_fi(nin)%P(1,1) ,is2pt_fi(nin)%P(1),
4115 + jg ,nfic ,fici ,fics ,ixss )
4116 DO j = 1,4
4117 j1 = ixss(j)
4118 niskyfil = niskyfil + 1
4119 fskyfi(nin)%P(1,niskyfil)=-fics(j,1)
4120 fskyfi(nin)%P(2,niskyfil)=-fics(j,2)
4121 fskyfi(nin)%P(3,niskyfil)=-fics(j,3)
4122 fskyfi(nin)%P(4,niskyfil)= fics(j,4)
4123 iskyfi(nin)%P(niskyfil) = j1
4124 END DO
4125 ELSE
4126 niskyfil = niskyfil + 1
4127 fskyfi(nin)%P(1,niskyfil)=-fxi(i)
4128 fskyfi(nin)%P(2,niskyfil)=-fyi(i)
4129 fskyfi(nin)%P(3,niskyfil)=-fzi(i)
4130 fskyfi(nin)%P(4,niskyfil)= stif(i)
4131 iskyfi(nin)%P(niskyfil) = ig
4132 ENDIF
4133C must be be done for second/plyxfem node if it must be considered -----
4134 ENDIF
4135 ENDDO
4136 ENDIF
4137C Thermique
4138 ELSE
4139 IF(intply == 0) THEN
4140 DO i=1,jlt
4141 IF (h1(i)/=0.) THEN
4142 niskyl = niskyl + 1
4143 fskyi(niskyl,1)=fx1(i)
4144 fskyi(niskyl,2)=fy1(i)
4145 fskyi(niskyl,3)=fz1(i)
4146 fskyi(niskyl,4)=stif(i)*abs(h1(i))
4147 isky(niskyl) = ix1(i)
4148 ftheskyi(niskyl) = phi1(i)
4149 ENDIF
4150 ENDDO
4151 DO i=1,jlt
4152 IF (h2(i)/=zero) THEN
4153 niskyl = niskyl + 1
4154 fskyi(niskyl,1)=fx2(i)
4155 fskyi(niskyl,2)=fy2(i)
4156 fskyi(niskyl,3)=fz2(i)
4157 fskyi(niskyl,4)=stif(i)*abs(h2(i))
4158 isky(niskyl) = ix2(i)
4159 ftheskyi(niskyl) = phi2(i)
4160 ENDIF
4161 ENDDO
4162 DO i=1,jlt
4163 IF (h3(i)/=zero) THEN
4164 niskyl = niskyl + 1
4165 fskyi(niskyl,1)=fx3(i)
4166 fskyi(niskyl,2)=fy3(i)
4167 fskyi(niskyl,3)=fz3(i)
4168 fskyi(niskyl,4)=stif(i)*abs(h3(i))
4169 isky(niskyl) = ix3(i)
4170 ftheskyi(niskyl) = phi3(i)
4171 ENDIF
4172 ENDDO
4173 DO i=1,jlt
4174 IF (h4(i)/=zero) THEN
4175 niskyl = niskyl + 1
4176 fskyi(niskyl,1)=fx4(i)
4177 fskyi(niskyl,2)=fy4(i)
4178 fskyi(niskyl,3)=fz4(i)
4179 fskyi(niskyl,4)=stif(i)*abs(h4(i))
4180 isky(niskyl) = ix4(i)
4181 ftheskyi(niskyl) = phi4(i)
4182 ENDIF
4183 ENDDO
4184C
4185 DO i=1,jlt
4186 IF(hh(i)==zero)cycle
4187 ig =nsvg(i)
4188 IF(ig>0) THEN
4189 IF (ig >numnod) THEN
4190 jg = ig - numnod
4191 nfic=1
4192 fici(nfic)=fxi(i)
4193 nfic = nfic + 1
4194 fici(nfic)=fyi(i)
4195 nfic = nfic + 1
4196 fici(nfic)=fzi(i)
4197 nfic = nfic + 1
4198 fici(nfic)=stif(i)
4199 nfic = nfic + 1
4200 fici(nfic)=phi(i)
4201 CALL i24forc_fic(3 ,irtse ,nsne ,is2se ,is2pt ,
4202 + jg ,nfic ,fici ,ficsth ,ixss )
4203 DO j = 1,4
4204 j1 = ixss(j)
4205 niskyl = niskyl + 1
4206 fskyi(niskyl,1)=-ficsth(j,1)
4207 fskyi(niskyl,2)=-ficsth(j,2)
4208 fskyi(niskyl,3)=-ficsth(j,3)
4209 fskyi(niskyl,4)= ficsth(j,4)
4210 ftheskyi(niskyl)=ficsth(j,5)
4211 isky(niskyl) = j1
4212 END DO
4213 ELSE
4214 niskyl = niskyl + 1
4215 fskyi(niskyl,1)=-fxi(i)
4216 fskyi(niskyl,2)=-fyi(i)
4217 fskyi(niskyl,3)=-fzi(i)
4218 fskyi(niskyl,4)= stif(i)
4219 isky(niskyl) = ig
4220 ftheskyi(niskyl)=phi(i)
4221 END IF !(IG >NUMNOD) THEN
4222 ELSE
4223 ig = -ig
4224 IF(isedge_fi(nin)%P(ig)==1)THEN
4225 jg = ig
4226 nfic=1
4227 fici(nfic)=fxi(i)
4228 nfic = nfic + 1
4229 fici(nfic)=fyi(i)
4230 nfic = nfic + 1
4231 fici(nfic)=fzi(i)
4232 nfic = nfic + 1
4233 fici(nfic)=stif(i)
4234 nfic = nfic + 1
4235 fici(nfic)=phi(i)
4236 CALL i24forc_fic(3 ,irtse_fi(nin)%P(1,1) ,nsne ,is2se_fi(nin)%P(1,1) ,is2pt_fi(nin)%P(1) ,
4237 + jg ,nfic ,fici ,ficsth ,ixss )
4238
4239 DO j = 1,4
4240 j1 = ixss(j)
4241 niskyfil = niskyfil + 1
4242 fskyfi(nin)%P(1,niskyfil)=-ficsth(j,1)
4243 fskyfi(nin)%P(2,niskyfil)=-ficsth(j,2)
4244 fskyfi(nin)%P(3,niskyfil)=-ficsth(j,3)
4245 fskyfi(nin)%P(4,niskyfil)= ficsth(j,4)
4246 ftheskyfi(nin)%P(niskyfil)=ficsth(j,5)
4247 isky(niskyl) = j1
4248 END DO
4249 ELSE
4250 niskyfil = niskyfil + 1
4251 fskyfi(nin)%P(1,niskyfil)=-fxi(i)
4252 fskyfi(nin)%P(2,niskyfil)=-fyi(i)
4253 fskyfi(nin)%P(3,niskyfil)=-fzi(i)
4254 fskyfi(nin)%P(4,niskyfil)= stif(i)
4255 iskyfi(nin)%P(niskyfil) = ig
4256 ftheskyfi(nin)%P(niskyfil)=phi(i)
4257 ENDIF
4258 ENDIF
4259 ENDDO
4260C with plyxfem formulation
4261 ELSE
4262 DO i=1,jlt
4263 IF (h1(i)/=0.) THEN
4264 niskyl = niskyl + 1
4265 fskyi(niskyl,1)=fx1(i)
4266 fskyi(niskyl,2)=fy1(i)
4267 fskyi(niskyl,3)=fz1(i)
4268 fskyi(niskyl,4)=stif(i)*abs(h1(i))
4269 isky(niskyl) = ix1(i)
4270 ftheskyi(niskyl) = phi1(i)
4271C
4272 plyskyi%FSKYI(niskyl,1)=fx1(i)
4273 plyskyi%FSKYI(niskyl,2)=fy1(i)
4274 plyskyi%FSKYI(niskyl,3)=fz1(i)
4275 plyskyi%FSKYI(niskyl,4)=stif(i)*abs(h1(i))
4276 plyskyi%FSKYI(niskyl,5)=iply(1,i)
4277 ENDIF
4278 ENDDO
4279 DO i=1,jlt
4280 IF (h2(i)/=zero) THEN
4281 niskyl = niskyl + 1
4282 fskyi(niskyl,1)=fx2(i)
4283 fskyi(niskyl,2)=fy2(i)
4284 fskyi(niskyl,3)=fz2(i)
4285 fskyi(niskyl,4)=stif(i)*abs(h2(i))
4286 isky(niskyl) = ix2(i)
4287 ftheskyi(niskyl) = phi2(i)
4288C
4289 plyskyi%FSKYI(niskyl,1)=fx2(i)
4290 plyskyi%FSKYI(niskyl,2)=fy2(i)
4291 plyskyi%FSKYI(niskyl,3)=fz2(i)
4292 plyskyi%FSKYI(niskyl,4)=stif(i)*abs(h2(i))
4293 plyskyi%FSKYI(niskyl,5)=iply(2,i)
4294 ENDIF
4295 ENDDO
4296 DO i=1,jlt
4297 IF (h3(i)/=zero) THEN
4298 niskyl = niskyl + 1
4299 fskyi(niskyl,1)=fx3(i)
4300 fskyi(niskyl,2)=fy3(i)
4301 fskyi(niskyl,3)=fz3(i)
4302 fskyi(niskyl,4)=stif(i)*abs(h3(i))
4303 isky(niskyl) = ix3(i)
4304 ftheskyi(niskyl) = phi3(i)
4305C
4306 plyskyi%FSKYI(niskyl,1)=fx3(i)
4307 plyskyi%FSKYI(niskyl,2)=fy3(i)
4308 plyskyi%FSKYI(niskyl,3)=fz3(i)
4309 plyskyi%FSKYI(niskyl,4)=stif(i)*abs(h3(i))
4310 plyskyi%FSKYI(niskyl,5)=iply(3,i)
4311 ENDIF
4312 ENDDO
4313 DO i=1,jlt
4314 IF (h4(i)/=zero) THEN
4315 niskyl = niskyl + 1
4316 fskyi(niskyl,1)=fx4(i)
4317 fskyi(niskyl,2)=fy4(i)
4318 fskyi(niskyl,3)=fz4(i)
4319 fskyi(niskyl,4)=stif(i)*abs(h4(i))
4320 isky(niskyl) = ix4(i)
4321 ftheskyi(niskyl) = phi4(i)
4322C
4323 plyskyi%FSKYI(niskyl,1)=fx4(i)
4324 plyskyi%FSKYI(niskyl,2)=fy4(i)
4325 plyskyi%FSKYI(niskyl,3)=fz4(i)
4326 plyskyi%FSKYI(niskyl,4)=stif(i)*abs(h4(i))
4327 plyskyi%FSKYI(niskyl,5)=iply(4,i)
4328 ENDIF
4329 ENDDO
4330C
4331 DO i=1,jlt
4332 IF(hh(i)==zero)cycle
4333 ig =nsvg(i)
4334 IF(ig>0) THEN
4335 IF (ig >numnod) THEN
4336 jg = ig - numnod
4337 nfic=1
4338 fici(nfic)=fxi(i)
4339 nfic = nfic + 1
4340 fici(nfic)=fyi(i)
4341 nfic = nfic + 1
4342 fici(nfic)=fzi(i)
4343 nfic = nfic + 1
4344 fici(nfic)=stif(i)
4345 nfic = nfic + 1
4346 fici(nfic)=phi(i)
4347 CALL i24forc_fic(3 ,irtse ,nsne ,is2se ,is2pt ,
4348 + jg ,nfic ,fici ,ficsth ,ixss )
4349 DO j = 1,4
4350 j1 = ixss(j)
4351 niskyl = niskyl + 1
4352 fskyi(niskyl,1)=-ficsth(j,1)
4353 fskyi(niskyl,2)=-ficsth(j,2)
4354 fskyi(niskyl,3)=-ficsth(j,3)
4355 fskyi(niskyl,4)= ficsth(j,4)
4356 ftheskyi(niskyl)=ficsth(j,5)
4357 isky(niskyl) = j1
4358 END DO
4359 ELSE
4360 niskyl = niskyl + 1
4361 fskyi(niskyl,1)=-fxi(i)
4362 fskyi(niskyl,2)=-fyi(i)
4363 fskyi(niskyl,3)=-fzi(i)
4364 fskyi(niskyl,4)= stif(i)
4365 isky(niskyl) = ig
4366 ftheskyi(niskyl)=phi(i)
4367 END IF !(IG >NUMNOD)
4368CC the second node is not plyxfem
4369cc PLYSKYI(1)%FSKYI(NISKYL,1)=-FXI(I)
4370cc PLYSKYI(1)%FSKYI(NISKYL,2)=-FYI(I)
4371cc PLYSKYI(1)%FSKYI(NISKYL,3)=-FZI(I)
4372cc PLYSKYI(I)%FSKYI(NISKYL,4)=STIF(I)
4373 ELSE
4374 ig = -ig
4375 IF(isedge_fi(nin)%P(ig)==1)THEN
4376 jg = ig
4377 nfic=1
4378 fici(nfic)=fxi(i)
4379 nfic = nfic + 1
4380 fici(nfic)=fyi(i)
4381 nfic = nfic + 1
4382 fici(nfic)=fzi(i)
4383 nfic = nfic + 1
4384 fici(nfic)=stif(i)
4385 nfic = nfic + 1
4386 fici(nfic)=phi(i)
4387 CALL i24forc_fic(3 ,irtse_fi(nin)%P(1,1) ,nsne ,is2se_fi(nin)%P(1,1) ,is2pt_fi(nin)%P(1) ,
4388 + jg ,nfic ,fici ,ficsth ,ixss )
4389 DO j = 1,4
4390 j1 = ixss(j)
4391 niskyfil = niskyfil + 1
4392 fskyfi(nin)%P(1,niskyfil)=-ficsth(j,1)
4393 fskyfi(nin)%P(2,niskyfil)=-ficsth(j,2)
4394 fskyfi(nin)%P(3,niskyfil)=-ficsth(j,3)
4395 fskyfi(nin)%P(4,niskyfil)= ficsth(j,4)
4396 ftheskyfi(nin)%P(niskyfil)=ficsth(j,5)
4397 iskyfi(nin)%P(niskyfil) = j1
4398 END DO
4399
4400 ELSE
4401 niskyfil = niskyfil + 1
4402 fskyfi(nin)%P(1,niskyfil)=-fxi(i)
4403 fskyfi(nin)%P(2,niskyfil)=-fyi(i)
4404 fskyfi(nin)%P(3,niskyfil)=-fzi(i)
4405 fskyfi(nin)%P(4,niskyfil)= stif(i)
4406 iskyfi(nin)%P(niskyfil) = ig
4407 ftheskyfi(nin)%P(niskyfil)=phi(i)
4408 ENDIF
4409 ENDIF
4410 ENDDO
4411 ENDIF
4412 ENDIF
4413C
4414 RETURN
debug_mod
Definition debug_mod.F:45
message_mod
Definition message_mod.F:1249
plyxfem_mod::plyskyi
type(ply_data), allocatable plyskyi
Definition plyxfem_mod.F:92
tri7box::ftheskyfi
type(real_pointer), dimension(:), allocatable ftheskyfi
Definition tri7box.F:449
tri7box::fskyfi
type(real_pointer2), dimension(:), allocatable fskyfi
Definition tri7box.F:459
tri7box::iskyfi
type(int_pointer), dimension(:), allocatable iskyfi
Definition tri7box.F:480
ancmsg
subroutine ancmsg(msgid, msgtype, anmode, i1, i2, i3, i4, i5, i6, i7, i8, i9, i10, i11, i12, i13, i14, i15, i16, i17, i18, i19, i20, r1, r2, r3, r4, r5, r6, r7, r8, r9, c1, c2, c3, c4, c5, c6, c7, c8, c9, prmode)
Definition message.F:889
arret
subroutine arret(nn)
Definition arret.F:87

◆ i24for3()

subroutine i24for3 ( integer jlt,
a,
v,
integer ibcc,
integer, dimension(*) icodt,
fsav,
gap,
fric,
ms,
visc,
viscf,
integer noint,
stfn,
integer, dimension(*) itab,
integer, dimension(mvsiz) cn_loc,
stiglo,
stifn,
stif,
fskyi,
integer, dimension(*) isky,
n1,
n2,
n3,
h1,
h2,
h3,
h4,
fcont,
pene,
integer, dimension(mvsiz) ix1,
integer, dimension(mvsiz) ix2,
integer, dimension(mvsiz) ix3,
integer, dimension(mvsiz) ix4,
integer, dimension(mvsiz) nsvg,
integer ivis2,
integer neltst,
integer ityptst,
dt2t,
integer, dimension(mvsiz) subtria,
gapv,
integer inacti,
integer, dimension(mvsiz) index,
integer niskyfi,
integer, dimension(*) kinet,
integer newfront,
integer isecin,
integer, dimension(*) nstrf,
secfcum,
x,
integer, dimension(4,*) irect,
integer, dimension(mvsiz) ce_loc,
integer mfrot,
integer ifq,
frot_p,
secnd_fr,
alpha0,
integer ibag,
integer, dimension(*) icontact,
integer, dimension(2,nsn) irtlm,
viscn,
vxi,
vyi,
vzi,
msi,
integer, dimension(*) kini,
integer nin,
integer nisub,
integer, dimension(*) lisub,
integer, dimension(*) addsubs,
integer, dimension(s_addsubm) addsubm,
integer, dimension(*) lisubs,
integer, dimension(*) lisubm,
fsavsub,
integer, dimension(*) cand_n,
integer ilagm,
integer, dimension(3) icurv,
fncont,
ftcont,
integer nsn,
x1,
x2,
x3,
x4,
y1,
y2,
y3,
y4,
z1,
z2,
z3,
z4,
xi,
yi,
zi,
integer iadm,
rcurvi,
rcontact,
acontact,
pcontact,
anglmi,
padm,
integer intth,
phi,
fthe,
ftheskyi,
temp,
tempi,
rstif,
integer iform,
mskyi_sms,
integer, dimension(*) iskyi_sms,
integer, dimension(*) nsms,
integer, dimension(*) cand_n_n,
pene_old,
stif_old,
integer, dimension(*) mbinflg,
integer ilev,
integer igsti,
kmin,
integer intply,
integer, dimension(4,*) iply,
integer, dimension(*) inod_pxfem,
nm1,
nm2,
nm3,
integer nrebou,
integer, dimension(5,*) irtse,
integer nsne,
integer, dimension(2,*) is2se,
integer, dimension(*) is2pt,
integer, dimension(*) msegtyp,
integer jtask,
integer, dimension(nisubmax+1) isensint,
fsavparit,
integer nft,
type(h3d_database) h3d_data,
fricc,
viscffric,
fric_coefs,
integer, dimension(6,*) t2main_sms,
integer intnitsche,
forneqsi,
integer iorthfric,
fric_coefs2,
fricc2,
viscffric2,
integer nforth,
integer nfisot,
integer, dimension(mvsiz) indexorth,
integer, dimension(mvsiz) indexisot,
dir1,
dir2,
t2fac_sms,
f_pfit,
integer, dimension(nloadp_hyd_inter,numnod) tagncont,
integer, dimension(ninter+1), intent(in) kloadpinter,
integer, dimension(s_loadpinter), intent(in) loadpinter,
integer, dimension(nloadp_hyd), intent(in) loadp_hyd_inter,
integer, dimension(s_typsub) typsub,
integer, dimension(*) inflg_subs,
integer, dimension(*) inflg_subm,
integer, intent(in) ninloadp,
dimension(s_loadpinter), intent(in) dgaploadint,
integer, intent(in) s_loadpinter,
dimension(mvsiz), intent(inout) dist,
integer, dimension(mvsiz,13), intent(in) ixx,
integer, intent(in) interefric,
integer, intent(in) intcarea,
type (parameters_), intent(in) parameters,
intent(in) penref,
intent(in) kmax,
integer s_addsubm,
integer s_lisubm,
integer s_typsub,
integer nisubmax,
integer i_stok,
integer, intent(in) nrtm,
integer, intent(in) nrtse,
integer, intent(in) nsnr )

Definition at line 46 of file i24for3.F.

83C-----------------------------------------------
84C M o d u l e s
85C-----------------------------------------------
86 USE tri7box
87 USE plyxfem_mod
88 USE h3d_mod
89 USE anim_mod
90 USE outputs_mod
91 USE parameters_mod
92 USE i24intarea_fic_mod , ONLY : i24intarea_fic
93C-----------------------------------------------
94C I m p l i c i t T y p e s
95C-----------------------------------------------
96#include "implicit_f.inc"
97#include "comlock.inc"
98C-----------------------------------------------
99C G l o b a l P a r a m e t e r s
100C-----------------------------------------------
101#include "mvsiz_p.inc"
102C-----------------------------------------------
103C C o m m o n B l o c k s
104C-----------------------------------------------
105#include "com01_c.inc"
106#include "com04_c.inc"
107#include "com06_c.inc"
108#include "com08_c.inc"
109#include "scr05_c.inc"
110#include "scr07_c.inc"
111#include "scr11_c.inc"
112#include "scr14_c.inc"
113#include "scr16_c.inc"
114#include "scr18_c.inc"
115#include "sms_c.inc"
116#include "parit_c.inc"
117#include "param_c.inc"
118#include "impl1_c.inc"
119C-----------------------------------------------
120C D u m m y A r g u m e n t s
121C-----------------------------------------------
122 INTEGER S_ADDSUBM ! Size of ADDSUBM
123 INTEGER S_LISUBM ! Size of LISUBM
124 INTEGER S_TYPSUB ! Size of TYPSUB
125 INTEGER NISUBMAX ! Size of ISENSINT
126 INTEGER I_STOK ! Number of contact candidates / Size of FSAVPARIT
127 INTEGER NELTST,ITYPTST,JLT,IBCC,IBCM,IBCS,IVIS2,INACTI,IBAG,NIN,
128 . INTNITSCHE,IORTHFRIC,NFORTH ,NFISOT,
129 . NSN,ICODT(*), ITAB(*), ISKY(*), KINET(*),IRTLM(2,NSN),
130 . MFROT, IFQ, NOINT,NEWFRONT,ISECIN, NSTRF(*),
131 . IRECT(4,*),ICONTACT(*), CAND_N(*),
132 . KINI(*),SUBTRIA(MVSIZ),MBINFLG(*),ILEV,
133 . ISET, NISKYFI,IADM,INTTH,IFORM,CAND_N_N(*),INTPLY,
134 . MSEGTYP(*),TAGNCONT(NLOADP_HYD_INTER,NUMNOD),
135 . TYPSUB(S_TYPSUB),JTASK
136 INTEGER IX1(MVSIZ), IX2(MVSIZ), IX3(MVSIZ), IX4(MVSIZ),
137 . CN_LOC(MVSIZ),CE_LOC(MVSIZ),INDEX(MVSIZ),NSVG(MVSIZ),
138 . NISUB, LISUB(*), ADDSUBS(*), ADDSUBM(S_ADDSUBM), LISUBS(*),
139 . LISUBM(*),ILAGM,ICURV(3), NREBOU,
140 . ISKYI_SMS(*), NSMS(*),IGSTI,IPLY(4,*),INOD_PXFEM(*),
141 . IRTSE(5,*),NSNE ,IS2SE(2,*),IS2PT(*), ISENSINT(NISUBMAX+1),NFT,T2MAIN_SMS(6,*),
142 . INDEXORTH(MVSIZ) ,INDEXISOT(MVSIZ),INFLG_SUBS(*), INFLG_SUBM(*)
143 INTEGER , INTENT(IN) :: NINLOADP,S_LOADPINTER
144 INTEGER , INTENT(IN) :: KLOADPINTER(NINTER+1),LOADPINTER(S_LOADPINTER),
145 . LOADP_HYD_INTER(NLOADP_HYD)
146 INTEGER , INTENT(IN) :: IXX(MVSIZ,13)
147 INTEGER , INTENT(IN) :: INTEREFRIC, INTCAREA
148 INTEGER , INTENT(IN) :: NRTM, NRTSE, NSNR
149 my_real
150 . stiglo,frot_p(*), x(3,*),
151 . a(3,*), ms(*), v(3,*), fsav(*),fcont(3,*),
152 . secnd_fr(6,*),alpha0,
153 . gap, fric,visc,viscf,vis,dt2t,stfn(*),stifn(*),
154 . fskyi(lskyi,nfskyi),fsavsub(nthvki,*), fncont(3,*),ftcont(3,*),
155 . mskyi_sms(*),kmin
156 my_real
157 . la(mvsiz), lb(mvsiz), lc(mvsiz),stif(mvsiz),
158 . gapv(mvsiz), pene(mvsiz),
159 . secfcum(7,numnod,nsect),
160 . tmp(mvsiz),
161 . stifsav(mvsiz), viscn(*),
162 . vxi(mvsiz),vyi(mvsiz),vzi(mvsiz),msi(mvsiz),
163 . x1(mvsiz),y1(mvsiz),z1(mvsiz),
164 . x2(mvsiz),y2(mvsiz),z2(mvsiz),
165 . x3(mvsiz),y3(mvsiz),z3(mvsiz),
166 . x4(mvsiz),y4(mvsiz),z4(mvsiz),
167 . xi(mvsiz),yi(mvsiz),zi(mvsiz),
168 . n1(mvsiz), n2(mvsiz), n3(mvsiz),
169 . h1(mvsiz), h2(mvsiz), h3(mvsiz), h4(mvsiz),
170 . nm1(mvsiz), nm2(mvsiz), nm3(mvsiz),
171 . phi(mvsiz), fthe(*),ftheskyi(lskyi),temp(*), tempi(mvsiz),
172 . fsavparit(nisub+1,11,*),
173 . fricc(mvsiz),viscffric(mvsiz),fric_coefs(mvsiz,10),forneqsi(mvsiz,3),
174 . fric_coefs2(mvsiz,10),fricc2(mvsiz),viscffric2(mvsiz),
175 . dir1(mvsiz,3),dir2(mvsiz,3),t2fac_sms(*)
176
177
178 my_real
179 . rcurvi(*), rcontact(*), acontact(*),
180 . pcontact(*),padm, anglmi(*),rstif,
181 . pene_old(5,nsn),stif_old(2,nsn),f_pfit
182 my_real , INTENT(IN) :: dgaploadint(s_loadpinter)
183 my_real , INTENT(INOUT) :: dist(mvsiz)
184 my_real, DIMENSION(MVSIZ), INTENT(IN) :: penref
185 my_real, INTENT(IN) :: kmax
186
187 TYPE(H3D_DATABASE) :: H3D_DATA
188 TYPE (PARAMETERS_) ,INTENT(IN):: PARAMETERS
189C-----------------------------------------------
190C L o c a l V a r i a b l e s
191C-----------------------------------------------
192 INTEGER I, J1, IG, J, JG , K0,NBINTER,K1S,K,IL,IE, NN, NI,
193 . NA1,NA2,N,IMS2,ITAG,NN1,NN2,NN3,
194 . NN4,ii,IFIT,PP,PPL,ITYPSUB,ISS1,ISS2,IMS1,IX,IGRN
195 my_real
196 . fxr(mvsiz), fyr(mvsiz), fzr(mvsiz)
197 my_real
198 . fxi(mvsiz), fyi(mvsiz), fzi(mvsiz), fni(mvsiz),
199 . fxt(mvsiz),fyt(mvsiz),fzt(mvsiz),
200 . fx1(mvsiz), fx2(mvsiz), fx3(mvsiz), fx4(mvsiz),
201 . fy1(mvsiz), fy2(mvsiz), fy3(mvsiz), fy4(mvsiz),
202 . fz1(mvsiz), fz2(mvsiz), fz3(mvsiz), fz4(mvsiz),
203 . vis2(mvsiz), dtmi(mvsiz), xmu(mvsiz),stif0(mvsiz),
204 . vx(mvsiz), vy(mvsiz), vz(mvsiz), vn(mvsiz),
205 . xp(mvsiz), yp(mvsiz), zp(mvsiz),efric_l(mvsiz),
206 . vnx, vny, vnz, aa, crit,s2,
207 . v2, fm2, dt1inv, visca, fac,ff,alphi,alpha,beta,
208 . fx, fy, fz, f2, mas2, m2sk, dtmi0,dti,ft,fn,fmax,ftn,
209 . facm1, econtt, econvt, h0, econtdt,
210 . d1,d2,d3,d4,a1,a2,a3,a4,
211 . fsav1, fsav2, fsav3, fsav4, fsav5, fsav6, fsav7, fsav8,
212 . fsav9, fsav10, fsav11, fsav12, fsav13, fsav14, fsav15,
213 . fsav22, fsav23, fsav24, ffo,fsav29,
214 . e10, h0d, s2d, sum,
215 . la1d,la2d,la3d,la4d,t1,t1d,t2,t2d,ffd,visd,facd,d1d,
216 . p1s(mvsiz),p2s(mvsiz),p3s(mvsiz),p4s(mvsiz),
217 . d2d,d3d,d4d,vnxd,vnyd,vnzd,v2d,fm2d,f2d,aad,fxd,fyd,fzd,
218 . a1d,a2d,a3d,a4d,vv,ax1,ax2,ay1,ay2,az1,az2,ax,ay,az,
219 . area,p,vv1,vv2,v21,dmu, dti2,h00 ,a0x,a0y,a0z,rx,ry,rz,
220 . anx,any,anz,aan,aax,aay,aaz ,rr,rs,aaa,bbb,stfr,visr ,tm,ts,
221 . vn1(3),vn2(3),vn3(3),vn4(4),
222 . csa ,sna ,alphaf ,ftt1 ,ftt2 ,
223 . an(nforth) ,bn(nforth) ,nep1 ,nep2 ,nep ,c11 ,c22 ,ans ,ep ,
224 . signc,dgapload,gapp,efric_ls,efric_lm
225 my_real
226 . prec,facv
227 my_real
228 . st1(mvsiz),st2(mvsiz),st3(mvsiz),st4(mvsiz),stv(mvsiz),
229 . kt(mvsiz),c(mvsiz),cf(mvsiz),dpene(mvsiz),
230 . ks(mvsiz),k1(mvsiz),k2(mvsiz),k3(mvsiz),k4(mvsiz),
231 . cs(mvsiz),c1(mvsiz),c2(mvsiz),c3(mvsiz),c4(mvsiz),
232 . cx,cy,cfi,aux, fx6(6,mvsiz), fy6(6,mvsiz), fz6(6,mvsiz),
233 . phi1(mvsiz), phi2(mvsiz), phi3(mvsiz), phi4(mvsiz),
234 . dtmini,stif0_imp(mvsiz),penn,tiny,pendr,
235 . fnnit1,fnnit2,fnnit3,fnnitsche,t3,fnitxt(mvsiz),fnityt(mvsiz),fnitzt(mvsiz),
236 . xmu2(mvsiz),facn,fnon,stifkt(mvsiz),fackt,arean_fic
237 INTEGER JSUB,KSUB,JJ,KK,IN,NSUB,ibug,ip,NCFIT
238 INTEGER ISEGTYP(MVSIZ),IXSS(4),NOD1,NOD2,NS
239 my_real
240 . fsavsub1(25,nisub),impx,impy,impz,vnm(mvsiz),sfac,stmin,
241 . fxs(4),fys(4),fzs(4),hh,finc,fa,fb,ddp,prec1
242
243 DOUBLE PRECISION RP1,RP2,STIF_N
244 INTEGER TAGIP(MVSIZ),ISPDO
245C
246 INTEGER BITGET
247 EXTERNAL bitget
248C-----------------------------------------------
249 dtmini=zero
250 fac = zero
251 IF (iresp==1) THEN
252 prec = fiveem4
253 tiny = em15
254 prec1= em05
255 ELSE
256 prec = em10
257 tiny = em30
258 prec1= zero
259 ENDIF
260 IF(dt1>zero)THEN
261 dt1inv = one/dt1
262 ELSE
263 dt1inv =zero
264 ENDIF
265 econtt = zero
266 econvt = zero
267 econtdt = zero
268 DO i=1,jlt
269 stif0(i) = stif(i)
270 stifkt(i) = stif(i)
271 ENDDO
272 efric_l(1:jlt) = zero
273C---------------------
274C COURBURE FIXE
275C---------------------
276 IF(icurv(1)==1)THEN
277 ELSEIF(icurv(1)==2)THEN
278 ELSEIF(icurv(1) == 3)THEN
279 ENDIF
280 ncfit = icurv(2)
281 IF (inacti==-1.AND.impl_s==0.AND.ncfit>0) THEN
282 ifit = 1
283 ELSE
284 ifit = 0
285 END IF
286 fnon = ep02
287 fmax = ep03
288 IF (igsti==-1) fmax = kmax
289C-------------------------------------------
290C PENE Offset removed to i24dst3.F
291C-------------------------------------------
292C
293 IF(intply == 0) THEN
294 DO i=1,jlt
295 IF(pene(i) == zero)THEN
296 h1(i) = zero
297 h2(i) = zero
298 h3(i) = zero
299 h4(i) = zero
300 fni(i)= zero
301C
302 fxi(i)=zero
303 fyi(i)=zero
304 fzi(i)=zero
305C
306 fx1(i)=zero
307 fy1(i)=zero
308 fz1(i)=zero
309C
310 fx2(i)=zero
311 fy2(i)=zero
312 fz2(i)=zero
313C
314 fx3(i)=zero
315 fy3(i)=zero
316 fz3(i)=zero
317C
318 fx4(i)=zero
319 fy4(i)=zero
320 fz4(i)=zero
321 ELSE
322 vx(i) = vxi(i) - h1(i)*v(1,ix1(i)) - h2(i)*v(1,ix2(i))
323 . - h3(i)*v(1,ix3(i)) - h4(i)*v(1,ix4(i))
324 vy(i) = vyi(i) - h1(i)*v(2,ix1(i)) - h2(i)*v(2,ix2(i))
325 . - h3(i)*v(2,ix3(i)) - h4(i)*v(2,ix4(i))
326 vz(i) = vzi(i) - h1(i)*v(3,ix1(i)) - h2(i)*v(3,ix2(i))
327 . - h3(i)*v(3,ix3(i)) - h4(i)*v(3,ix4(i))
328 vn(i) = n1(i)*vx(i) + n2(i)*vy(i) + n3(i)*vz(i)
329 vnm(i) = nm1(i)*vx(i) + nm2(i)*vy(i) + nm3(i)*vz(i)
330 ENDIF
331 ENDDO
332 ELSE
333 DO i=1,jlt
334 IF(pene(i) == zero)THEN
335 h1(i) = zero
336 h2(i) = zero
337 h3(i) = zero
338 h4(i) = zero
339 fni(i)= zero
340C
341 fxi(i)=zero
342 fyi(i)=zero
343 fzi(i)=zero
344C
345 fx1(i)=zero
346 fy1(i)=zero
347 fz1(i)=zero
348C
349 fx2(i)=zero
350 fy2(i)=zero
351 fz2(i)=zero
352C
353 fx3(i)=zero
354 fy3(i)=zero
355 fz3(i)=zero
356C
357 fx4(i)=zero
358 fy4(i)=zero
359 fz4(i)=zero
360 ELSE
361 jj = iply(1,i)
362 nn1 =inod_pxfem(ix1(i))
363 nn2 =inod_pxfem(ix2(i))
364 nn3 =inod_pxfem(ix3(i))
365 nn4 =inod_pxfem(ix4(i))
366C
367 vn1(1) = v(1,ix1(i))
368 vn1(2) = v(2,ix1(i))
369 vn1(3) = v(3,ix1(i))
370C
371 vn2(1) = v(1,ix2(i))
372 vn2(2) = v(2,ix2(i))
373 vn2(3) = v(3,ix2(i))
374C
375 vn3(1) = v(1,ix3(i))
376 vn3(2) = v(2,ix3(i))
377 vn3(3) = v(3,ix3(i))
378C
379 vn4(1) = v(1,ix4(i))
380 vn4(2) = v(2,ix4(i))
381 vn4(3) = v(3,ix4(i))
382C
383 IF(iplyxfem /= 2)THEN
384C
385C possible energy creation after full delamination
386C - seen with IPLYXFEM=1 -
387 jj = iply(1,i)
388 IF(nn1 > 0 .AND. jj > 0) THEN
389 vn1(1) = vn1(1) + ply(jj)%V(1,nn1)
390 vn1(2) = vn1(2) + ply(jj)%V(2,nn1)
391 vn1(3) = vn1(3) + ply(jj)%V(3,nn1)
392 ENDIF
393 jj = iply(2,i)
394 IF(nn2 > 0 .AND. jj > 0) THEN
395 vn2(1) = vn2(1) + ply(jj)%V(1,nn2)
396 vn2(2) = vn2(2) + ply(jj)%V(2,nn2)
397 vn2(3) = vn2(3) + ply(jj)%V(3,nn2)
398 ENDIF
399 jj = iply(3,i)
400 IF(nn3 > 0 .AND. jj > 0) THEN
401 vn3(1) = vn3(1) + ply(jj)%V(1,nn3)
402 vn3(2) = vn3(2) + ply(jj)%V(2,nn3)
403 vn3(3) = vn3(3) + ply(jj)%V(3,nn3)
404 ENDIF
405 jj = iply(4,i)
406 IF(nn4 > 0 .AND. jj > 0) THEN
407 vn4(1) = vn4(1) + ply(jj)%V(1,nn4)
408 vn4(2) = vn4(2) + ply(jj)%V(2,nn4)
409 vn4(3) = vn4(3) + ply(jj)%V(3,nn4)
410 ENDIF
411 END IF
412
413 vx(i) = vxi(i) - h1(i)*vn1(1) - h2(i)*vn2(1)
414 . - h3(i)*vn3(1) - h4(i)*vn4(1)
415 vy(i) = vyi(i) - h1(i)*vn1(2) - h2(i)*vn2(2)
416 . - h3(i)*vn3(2) - h4(i)*vn4(2)
417 vz(i) = vzi(i) - h1(i)*vn1(3) - h2(i)*vn2(3)
418 . - h3(i)*vn3(3) - h4(i)*vn4(3)
419 vn(i) = n1(i)*vx(i) + n2(i)*vy(i) + n3(i)*vz(i)
420 vnm(i) = nm1(i)*vx(i) + nm2(i)*vy(i) + nm3(i)*vz(i)
421 ENDIF
422 ENDDO
423
424 ENDIF
425C-------------------------------------------
426C stifness reduction to avoid positive
427C force jump and energy generation
428C-------------------------------------------
429 IF (ifit >0.OR.(inacti==-1.AND.impl_s==0)) THEN
430 DO i=1,jlt
431 IF(pene(i) == zero)cycle
432 stif(i) = stif0(i)
433 IF(stiglo<=zero) stif(i) = half*stif(i)
434C---- nonlinear stif
435 IF(penref(i) > zero) THEN
436 pendr = (pene(i)/penref(i))**2
437 fackt = min(fmax,(one+three*fnon*pendr))
438 facn = min(fmax,(one+fnon*pendr))
439 stifkt(i)= stifkt(i)*fackt
440 stif(i)= stif(i)*facn
441 END IF
442 fni(i)= -stif(i) * pene(i)
443 ENDDO
444 ELSEIF (igsti==-1.AND.impl_s==0) THEN
445c----- first version w/o force jump treatment
446 DO i=1,jlt
447 IF(pene(i) == zero)cycle
448 stif(i) = stif0(i)
449C---- nonlinear stif
450 IF(penref(i) > zero) THEN
451 pendr = (pene(i)/penref(i))**2
452 facn = min(fmax,(one+fnon*pendr))
453 stif(i)= min(kmax,stif(i)*facn)
454 fackt = three*(stif(i)/stif0(i)-one)
455 stifkt(i)= max(stif(i),stifkt(i)*fackt)
456 END IF
457 fni(i)= -stif(i) * pene(i)
458 ENDDO
459 ELSEIF(impl_s>0.AND.igsti==6)THEN
460C-------0.001*STIF initially (first contact), special case(rebound) 0.0001*STIF but >=KMIN------
461C--- NREBOU<0 rebound in this cycle, : -1 -> at first contact; -2 -> after first contact
462 IF (nrebou == -1) THEN
463 sfac = em04
464 ELSE
465 sfac = em03
466 END IF
467 DO i=1,jlt
468 IF(pene(i) == zero)cycle
469 jg = nsvg(i)
470 n = cand_n_n(i)
471 IF(jg > 0)THEN
472 IF (stif_old(1,n)==zero.OR.nrebou==-1) THEN
473 stif_old(1,n) = sfac*stif0(i)
474 ELSEIF (nrebou <-1) THEN
475 stmin = em04*stif0(i)
476 stif_old(1,n) = max(stmin,em01*stif_old(1,n))
477 END IF
478 stif_old(1,n) = max(kmin,stif_old(1,n))
479 stif(i) = stif_old(1,n)
480 ELSE
481 jg = -jg
482 IF (stif_oldfi(nin)%P(1,jg)==zero.OR.nrebou==-1) THEN
483 stif_oldfi(nin)%P(1,jg)= sfac*stif0(i)
484 ELSEIF (nrebou <-1) THEN
485 stmin = em04*stif0(i)
486 stif_oldfi(nin)%P(1,jg)=
487 + max(stmin,em01*stif_oldfi(nin)%P(1,jg))
488 END IF
489 stif_oldfi(nin)%P(1,jg)= max(kmin,stif_oldfi(nin)%P(1,jg))
490 stif(i) = stif_oldfi(nin)%P(1,jg)
491 END IF
492 END DO
493 IF(inconv >=0)THEN
494C---------if multi-contact in the same groupe---
495 DO i=1,jlt
496 IF(pene(i) == zero)cycle
497 jg = nsvg(i)
498 n = cand_n_n(i)
499 IF(jg > 0)THEN
500 stif0_imp(i) = stif_old(1,n)
501 ELSE
502 stif0_imp(i) = stif_oldfi(nin)%P(1,-jg)
503 END IF
504 END DO
505 DO i=1,jlt
506 IF(pene(i) == zero)cycle
507 jg = nsvg(i)
508 n = cand_n_n(i)
509 IF(jg > 0)THEN
510C---------if multi-contact in the same groupe---
511 IF(pene_old(2,n)/= zero.OR.pene_old(5,n)/= zero)THEN
512 IF(pene(i) > pene_old(2,n) )THEN
513 facm1 = pene(i)/max(em20,pene_old(2,n))
514 IF (stif0_imp(i)/stif0(i) <em02) THEN
515 facm1 = onep2*facm1
516 penn = three_half
517 ELSE
518 penn = onep05
519 END IF
520 facm1=min(penn,facm1)
521 stif(i) = facm1*stif0_imp(i)
522C-----------never higher than STIF0----------------
523 stif(i) = min(stif(i),stif0(i))
524 ENDIF
525C------------quand PENE_OLD(2,N)=PENE_OLD(1,N)
526 IF(inconv ==1 )THEN
527 pene_old(1,n) = pene(i)
528 stif_old(1,n) = stif(i)
529 ENDIF
530 ELSEIF(inconv ==1 )THEN
531c new impact, possible multi impact
532#include "lockon.inc"
533 pene_old(1,n) = max(pene_old(1,n),pene(i))
534 stif_old(1,n) = min(stif_old(1,n),stif(i))
535#include "lockoff.inc"
536 ENDIF
537 ELSE
538 jg = -jg
539 IF(pene_oldfi(nin)%P(2,jg)/= zero.OR.
540 . pene_oldfi(nin)%P(5,jg)/= zero)THEN
541 IF(pene(i) > pene_oldfi(nin)%P(2,jg) )THEN
542 facm1 = pene(i)/max(em20,pene_oldfi(nin)%P(2,jg))
543 IF (stif0_imp(i)/stif0(i) <em02) THEN
544 facm1 = onep2*facm1
545 penn = three_half
546 ELSE
547 penn = onep05
548 END IF
549 facm1=min(penn,facm1)
550 stif(i) = facm1*stif0_imp(i)
551 stif(i) = min(stif(i),stif0(i))
552 ENDIF
553 IF(inconv ==1 )THEN
554 pene_oldfi(nin)%P(1,jg) = pene(i)
555 stif_oldfi(nin)%P(1,jg) = stif(i)
556 END IF !(INCONV ==1 )THEN
557 ELSEIF(inconv ==1 )THEN
558#include "lockon.inc"
559 pene_oldfi(nin)%P(1,jg) =
560 * max(pene_oldfi(nin)%P(1,jg),pene(i))
561 stif_oldfi(nin)%P(1,jg) =
562 * min(stif_oldfi(nin)%P(1,jg),stif(i))
563#include "lockoff.inc"
564 ENDIF
565 ENDIF
566 ENDDO
567C-------to pass for implicit
568 DO i=1,jlt
569 IF(pene(i) == zero)cycle
570 jg = nsvg(i)
571 n = cand_n_n(i)
572 IF(jg > 0)THEN
573 stif_old(2,n) = stif(i)
574 ELSE
575 stif_oldfi(nin)%P(2,-jg) = stif(i)
576 END IF
577 END DO
578 END IF !(INCONV>=0)THEN
579 DO i=1,jlt
580 IF(pene(i) == zero)cycle
581 IF(stiglo<=zero)THEN
582 stif(i) = half*stif(i)
583ccc ECONTT = ECONTT + HALF*STIF(I)*PENE(I)**2
584 ELSEIF(stif(i)/=zero)THEN
585 stif(i) = stiglo
586ccc ECONTT = ECONTT + HALF*STIGLO*PENE(I)**2
587 ENDIF
588 fni(i)= -stif(i) * pene(i)
589 ENDDO
590
591C------explicit or implicit+Istif/=6
592 ELSE
593 DO i=1,jlt
594 IF(pene(i) == zero)cycle
595 dpene(i) = max(zero,-vnm(i)*dt1)
596 jg = nsvg(i)
597 n = cand_n_n(i)
598 IF(jg > 0)THEN
599 IF(tt /= zero.AND.(pene_old(2,n)/= zero.OR.pene_old(5,n)/= zero))THEN
600 ddp = pene(i)-pene_old(2,n)-dpene(i)
601 IF(pene(i) > pene_old(2,n) .AND.ddp>zero)THEN
602 rp1 = pene_old(2,n)/pene(i)
603 rp2 = dpene(i)/pene(i)
604 stif(i) = stif_old(2,n)*rp1+ stif(i)*rp2
605 ELSEIF(.NOT.(pene_old(2,n)== zero.AND.stif_old(2,n)==zero))THEN
606C-----------first contact STIF_OLD(2,N) should not be zero
607 stif(i) = stif_old(2,n)
608 END IF
609 IF(inconv ==1 )THEN
610 pene_old(1,n) = pene(i)
611 stif_old(1,n) = stif(i)
612 ENDIF
613 ELSEIF(inconv ==1 )THEN
614c new impact, possible multi impact
615#include "lockon.inc"
616 pene_old(1,n) = max(pene_old(1,n),pene(i))
617 stif_old(1,n) = max(stif_old(1,n),stif(i))
618#include "lockoff.inc"
619 ENDIF
620 ELSE
621 jg = -jg
622 IF(tt /= zero.AND.(pene_oldfi(nin)%P(2,jg)/= zero.OR.
623 . pene_oldfi(nin)%P(5,jg)/= zero))THEN
624 ddp = pene(i)-pene_oldfi(nin)%P(2,jg)-dpene(i)
625 IF(pene(i) > pene_oldfi(nin)%P(2,jg) .AND.ddp>zero)THEN
626 rp1 = pene_oldfi(nin)%P(2,jg)/pene(i)
627 rp2 = dpene(i)/pene(i)
628 stif(i) = stif_oldfi(nin)%P(2,jg)*rp1+ stif(i)*rp2
629 ELSEIF(.NOT.(pene_oldfi(nin)%P(2,jg)== zero.AND.
630 . stif_oldfi(nin)%P(2,jg)==zero))THEN
631 stif(i) = stif_oldfi(nin)%P(2,jg)
632 END IF
633 IF(inconv ==1 )THEN
634 pene_oldfi(nin)%P(1,jg) = pene(i)
635 stif_oldfi(nin)%P(1,jg) = stif(i)
636 ENDIF
637 ELSEIF(inconv ==1 )THEN
638#include "lockon.inc"
639 pene_oldfi(nin)%P(1,jg) =
640 * max(pene_oldfi(nin)%P(1,jg),pene(i))
641 stif_oldfi(nin)%P(1,jg) =
642 * max(stif_oldfi(nin)%P(1,jg),stif(i))
643#include "lockoff.inc"
644 ENDIF
645 ENDIF
646 ENDDO
647C-------------------------------------------
648 DO i=1,jlt
649 IF(pene(i) == zero)cycle
650 IF(stiglo<=zero)THEN
651 stif(i) = half*stif(i)
652 ELSEIF(stif(i)/=zero)THEN
653 stif(i) = stiglo
654 ENDIF
655 fni(i)= -stif(i) * pene(i)
656
657 ENDDO
658C-------------------------------------------
659 ENDIF !IF(IMPL_S>0.AND.IGSTI==6)
660C-------------------------------------------
661C Contact energy (elastic)
662C-------------------------------------------
663 DO i=1,jlt
664 IF(pene(i) == zero)cycle
665 econtt = econtt+half*stif(i)*pene(i)**2
666 ENDDO
667C
668 IF(intnitsche > 0 ) THEN
669 DO i=1,jlt
670 IF(pene(i) == zero)cycle
671 fnnit1 = forneqsi(i,1)*n1(i)
672 fnnit2 = forneqsi(i,2)*n2(i)
673 fnnit3 = forneqsi(i,3)*n3(i)
674 fnnitsche = fnnit1 + fnnit2 + fnnit3
675C
676 fni(i) = min(zero,fni(i) - half*fnnitsche)
677C
678 ENDDO
679 ENDIF
680C---------------------------------
681C DAMPING + FRIC
682C---------------------------------
683C goto 999
684 IF(visc/=zero)THEN
685 IF(ivis2==0)THEN
686C---------------------------------
687C VISC QUAD TYPE V227
688C---------------------------------
689 DO i=1,jlt
690 vis2(i) = two * stif(i) * msi(i)
691 ENDDO
692C---------------------------------
693C---------------------------------
694 IF(kdtint==0.AND.(idtmins/=2.AND.idtmins_int==0))THEN
695 DO i=1,jlt
696 IF(pene(i) == zero)cycle
697 fac = stif(i) / max(em30,stif(i))
698 vis = sqrt(vis2(i))
699 ff = fac * visc * vis
700 stif(i) = stif0(i) + ff * dt1inv
701 stif(i) = max(stif(i) ,fac*sqrt(viscffric(i))*vis*dt1inv)
702 ffo = ff
703 ff = ff * vn(i)
704C ECONVT = ECONVT + FF * VN(I) * DT1
705 fni(i) = fni(i) + ff
706 ENDDO
707 ELSE
708 DO i=1,jlt
709 IF(pene(i) == zero)cycle
710 fac = stif(i) / max(em30,stif(i))
711 vis = sqrt(vis2(i))
712 c(i)= fac * visc * vis
713 kt(i)= stif0(i)
714 stif(i) = stif0(i) + c(i) * dt1inv
715 ff = c(i) * vn(i)
716C ECONVT = ECONVT + FF * VN(I) * DT1
717 fni(i) = fni(i) + ff
718 cf(i) = fac*sqrt(viscffric(i))*vis
719 stif(i) = max(stif(i) ,cf(i)*dt1inv)
720 ENDDO
721 ENDIF
722 ELSEIF(ivis2==1.OR.ivis2==2)THEN
723C---------------------------------
724C VISC QUAD TYPE V227
725C---------------------------------
726 IF(ivis2==1)THEN
727 facv = two
728 ELSE
729 facv = four
730 END IF
731 DO i=1,jlt
732 IF(pene(i) == zero)cycle
733 mas2 = ms(ix1(i))*h1(i)
734 . + ms(ix2(i))*h2(i)
735 . + ms(ix3(i))*h3(i)
736 . + ms(ix4(i))*h4(i)
737 vis2(i) = facv* stif(i) * msi(i) * mas2 /
738 . max(em30,msi(i)+mas2)
739 ENDDO
740C---------------------------------
741C---------------------------------
742 IF(kdtint==0.AND.(idtmins/=2.AND.idtmins_int==0))THEN
743 DO i=1,jlt
744 IF(pene(i) == zero)cycle
745 fac = stif(i) / max(em30,stif(i))
746 vis = sqrt(vis2(i))
747 ff = fac * visc * vis
748 stif(i) = stif0(i) + ff * dt1inv
749 stif(i) = max(stif(i) ,fac*sqrt(viscffric(i))*vis*dt1inv)
750 ffo = ff
751 ff = ff * vn(i)
752C ECONTDT = ECONTDT + FF * VN(I) * DT1
753 fni(i) = fni(i) + ff
754 ENDDO
755 ELSE
756 DO i=1,jlt
757 IF(pene(i) == zero)cycle
758 fac = stif(i) / max(em30,stif(i))
759 vis = sqrt(vis2(i))
760 c(i)= fac * visc * vis
761 kt(i)= stif0(i)
762 stif(i) = stif(i) + c(i) * dt1inv
763 ff = c(i) * vn(i)
764C ECONTDT = ECONTDT + FF * VN(I) * DT1
765 fni(i) = fni(i) + ff
766 cf(i) = fac*sqrt(viscffric(i))*vis
767 stif(i) = max(stif(i) ,cf(i)*dt1inv)
768 ENDDO
769 ENDIF
770 ELSEIF(ivis2==3)THEN
771C---------------------------------
772C VISC QUAD = 0
773C---------------------------------
774 DO i=1,jlt
775 IF(pene(i) == zero)cycle
776 vis2(i) = two * stif(i) * msi(i)
777 fac = stif(i) / max(em30,stif(i))
778 vis = sqrt(vis2(i))
779 ff = fac * visc * vis
780 stif(i) = stif0(i) + two* ff * dt1inv
781 stif(i) = max(stif(i) ,fac*sqrt(viscffric(i))*vis*dt1inv)
782 ff = ff * vn(i)
783C ECONTDT = ECONTDT+ FF * VN(I) * DT1
784 fni(i) = fni(i) + ff
785 ENDDO
786 ELSEIF(ivis2==4)THEN
787C---------------------------------
788C VISC = 0
789C---------------------------------
790 DO i=1,jlt
791 IF(pene(i) == zero)cycle
792 fac = stif(i) / max(em30,stif(i))
793 vis2(i) = two* stif(i) * msi(i)
794 vis = sqrt(vis2(i))
795 stif(i) = max(stif0(i) ,fac*sqrt(viscffric(i))*vis*dt1inv)
796 ENDDO
797 ELSEIF(ivis2==5)THEN
798C---------------------------------
799C VISC = 2M/dt => pour visc < 1 , stable : dt < 2M/visc = dt
800C M=M1*M2/M1+M2 pour visc = 1 , choc elastique
801C pour visc = 0.5 , choc mou
802C---------------------------------
803 DO i=1,jlt
804 IF(pene(i) == zero)cycle
805 mas2 = ms(ix1(i))*h1(i)
806 . + ms(ix2(i))*h2(i)
807 . + ms(ix3(i))*h3(i)
808 . + ms(ix4(i))*h4(i)
809 vis2(i) = two* stif(i) * msi(i)
810 vis = 2. * visc * dt1inv * msi(i) * mas2 /
811 . max(em30,msi(i)+mas2)
812 stif(i) = max(stif0(i) ,fac*sqrt(viscffric(i)*vis2(i))*dt1inv)
813 ff = vis * vn(i)
814 econtdt = econtdt + min(zero,ff-fni(i)) * vn(i) * dt1
815 fni(i) = min(fni(i),ff)
816 ENDDO
817 ELSE
818 ENDIF
819 ELSE
820 DO i=1,jlt
821 IF(viscffric(i)/=zero)THEN
822
823 IF(ivis2==0)THEN
824C---------------------------------
825C VISC QUAD TYPE V227
826C---------------------------------
827 vis2(i) = two * stif(i) * msi(i)
828C---------------------------------
829C---------------------------------
830 IF(kdtint==0.AND.(idtmins/=2.AND.idtmins_int==0))THEN
831 IF(pene(i) == zero)cycle
832 fac = stif(i) / max(em30,stif(i))
833 vis = sqrt(vis2(i))
834 ff = fac * visc * vis
835 stif(i) = stif0(i) + ff * dt1inv
836 stif(i) = max(stif(i) ,fac*sqrt(viscffric(i))*vis*dt1inv)
837 ffo = ff
838 ff = ff * vn(i)
839C ECONTDT = ECONTDT + FF * VN(I) * DT1
840 fni(i) = fni(i) + ff
841 ELSE
842 IF(pene(i) == zero)cycle
843 fac = stif(i) / max(em30,stif(i))
844 vis = sqrt(vis2(i))
845 c(i)= fac * visc * vis
846 kt(i)= stif0(i)
847 stif(i) = stif0(i) + c(i) * dt1inv
848 ff = c(i) * vn(i)
849C ECONTDT = ECONTDT + FF * VN(I) * DT1
850 fni(i) = fni(i) + ff
851 cf(i) = fac*sqrt(viscffric(i))*vis
852 stif(i) = max(stif(i) ,cf(i)*dt1inv)
853 ENDIF
854 ELSEIF(ivis2==1.OR.ivis2==2)THEN
855C---------------------------------
856C VISC QUAD TYPE V227
857C---------------------------------
858 IF(ivis2==1)THEN
859 facv = two
860 ELSE
861 facv = four
862 END IF
863 IF(pene(i) == zero)cycle
864 mas2 = ms(ix1(i))*h1(i)
865 . + ms(ix2(i))*h2(i)
866 . + ms(ix3(i))*h3(i)
867 . + ms(ix4(i))*h4(i)
868 vis2(i) = facv* stif(i) * msi(i) * mas2 /
869 . max(em30,msi(i)+mas2)
870C---------------------------------
871C---------------------------------
872 IF(kdtint==0.AND.(idtmins/=2.AND.idtmins_int==0))THEN
873 IF(pene(i) == zero)cycle
874 fac = stif(i) / max(em30,stif(i))
875 vis = sqrt(vis2(i))
876 ff = fac * visc * vis
877 stif(i) = stif0(i) + ff * dt1inv
878 stif(i) = max(stif(i) ,fac*sqrt(viscffric(i))*vis*dt1inv)
879 ffo = ff
880 ff = ff * vn(i)
881C ECONTDT = ECONTDT + FF * VN(I) * DT1
882 fni(i) = fni(i) + ff
883 ELSE
884 IF(pene(i) == zero)cycle
885 fac = stif(i) / max(em30,stif(i))
886 vis = sqrt(vis2(i))
887 c(i)= fac * visc * vis
888 kt(i)= stif0(i)
889 stif(i) = stif(i) + c(i) * dt1inv
890 ff = c(i) * vn(i)
891C ECONTDT = ECONTDT + FF * VN(I) * DT1
892 fni(i) = fni(i) + ff
893 cf(i) = fac*sqrt(viscffric(i))*vis
894 stif(i) = max(stif(i) ,cf(i)*dt1inv)
895 ENDIF
896 ELSEIF(ivis2==2)THEN
897C---------------------------------
898C VISC QUAD TYPE
899C---------------------------------
900 IF(pene(i) == zero)cycle
901 vis2(i) = two* stif(i) * msi(i)
902 fac = stif(i) / max(em30,stif(i))
903 vis = sqrt(vis2(i))
904 ff = fac * visc * vis
905 stif(i) = stif0(i) + two * ff * dt1inv
906 stif(i) = max(stif(i) ,fac*sqrt(viscffric(i))*vis*dt1inv)
907 ff = ff * vn(i)
908C ECONTDT = ECONTDT + FF * VN(I) * DT1
909 fni(i) = fni(i) + ff
910 ELSEIF(ivis2==3)THEN
911C---------------------------------
912C VISC QUAD = 0
913C---------------------------------
914 IF(pene(i) == zero)cycle
915 vis2(i) = two * stif(i) * msi(i)
916 fac = stif(i) / max(em30,stif(i))
917 vis = sqrt(vis2(i))
918 ff = fac * visc * vis
919 stif(i) = stif0(i) + two* ff * dt1inv
920 stif(i) = max(stif(i) ,fac*sqrt(viscffric(i))*vis*dt1inv)
921 ff = ff * vn(i)
922C ECONTDT = ECONTDT + FF * VN(I) * DT1
923 fni(i) = fni(i) + ff
924 ELSEIF(ivis2==4)THEN
925C---------------------------------
926C VISC = 0
927C---------------------------------
928 IF(pene(i) == zero)cycle
929 fac = stif(i) / max(em30,stif(i))
930 vis2(i) = two* stif(i) * msi(i)
931 vis = sqrt(vis2(i))
932 stif(i) = max(stif0(i) ,fac*sqrt(viscffric(i))*vis*dt1inv)
933 ELSEIF(ivis2==5)THEN
934C---------------------------------
935C VISC = 2M/dt => pour visc < 1 , stable : dt < 2M/visc = dt
936C M=M1*M2/M1+M2 pour visc = 1 , choc elastique
937C pour visc = 0.5 , choc mou
938C---------------------------------
939 IF(pene(i) == zero)cycle
940 mas2 = ms(ix1(i))*h1(i)
941 . + ms(ix2(i))*h2(i)
942 . + ms(ix3(i))*h3(i)
943 . + ms(ix4(i))*h4(i)
944 vis2(i) = two* stif(i) * msi(i)
945 vis = 2. * visc * dt1inv * msi(i) * mas2 /
946 . max(em30,msi(i)+mas2)
947 stif(i) = max(stif0(i) ,fac*sqrt(viscffric(i)*vis2(i))*dt1inv)
948 ff = vis * vn(i)
949 econtdt = econtdt + min(zero,ff-fni(i)) * vn(i) * dt1
950 fni(i) = min(fni(i),ff)
951 ELSE
952 ENDIF
953
954
955 ELSE
956cbb initialisation du tableau VIS2 pour eviter des problemes
957 vis2(i) = zero
958 ENDIF
959 ENDDO
960 ENDIF
961
962C---------------------------------
963C SAUVEGARDE DE L'IMPULSION NORMALE
964C---------------------------------
965 fsav1 = zero
966 fsav2 = zero
967 fsav3 = zero
968 fsav8 = zero
969 fsav9 = zero
970 fsav10= zero
971 fsav11= zero
972 IF (ilev==2) THEN
973 DO i=1,jlt
974 IF(pene(i) == zero)cycle
975 ie=ce_loc(i)
976 ims2 = bitget(mbinflg(ie),1)
977 fxi(i)=n1(i)*fni(i)
978 fyi(i)=n2(i)*fni(i)
979 fzi(i)=n3(i)*fni(i)
980 impx=fxi(i)*dt12
981 impy=fyi(i)*dt12
982 impz=fzi(i)*dt12
983 fsav8 =fsav8 +abs(impx)
984 fsav9 =fsav9 +abs(impy)
985 fsav10=fsav10+abs(impz)
986 IF (ims2 >0 ) THEN
987 fsav1 =fsav1 -impx
988 fsav2 =fsav2 -impy
989 fsav3 =fsav3 -impz
990 fsav11=fsav11-fni(i)*dt12
991 ELSE
992 fsav1 =fsav1 +impx
993 fsav2 =fsav2 +impy
994 fsav3 =fsav3 +impz
995 fsav11=fsav11+fni(i)*dt12
996 END IF
997 ENDDO
998 ELSE
999 DO i=1,jlt
1000 IF(pene(i) == zero)cycle
1001 fxi(i)=n1(i)*fni(i)
1002 fyi(i)=n2(i)*fni(i)
1003 fzi(i)=n3(i)*fni(i)
1004 impx=fxi(i)*dt12
1005 impy=fyi(i)*dt12
1006 impz=fzi(i)*dt12
1007 fsav1 =fsav1 +impx
1008 fsav2 =fsav2 +impy
1009 fsav3 =fsav3 +impz
1010 fsav8 =fsav8 +abs(impx)
1011 fsav9 =fsav9 +abs(impy)
1012 fsav10=fsav10+abs(impz)
1013 fsav11=fsav11+fni(i)*dt12
1014 ENDDO
1015 END IF !(ILEV==2) THEN
1016#include "lockon.inc"
1017 fsav(1)=fsav(1)+fsav1
1018 fsav(2)=fsav(2)+fsav2
1019 fsav(3)=fsav(3)+fsav3
1020 fsav(8)=fsav(8)+fsav8
1021 fsav(9)=fsav(9)+fsav9
1022 fsav(10)=fsav(10)+fsav10
1023 fsav(11)=fsav(11)+fsav11
1024#include "lockoff.inc"
1025C
1026 IF(isensint(1)/=0) THEN
1027 DO i=1,jlt
1028 IF(pene(i) == zero)cycle
1029 fsavparit(1,1,i+nft) = fxi(i)
1030 fsavparit(1,2,i+nft) = fyi(i)
1031 fsavparit(1,3,i+nft) = fzi(i)
1032 ENDDO
1033 ENDIF
1034C---------------------------------
1035 IF((anim_v(12)+outp_v(12)+h3d_data%N_VECT_PCONT>0.AND.
1036 . ((tt>=tanim.AND.tt<=tanim_stop).OR.tt>=toutp.OR.(tt>=h3d_data%TH3D.AND.tt<=h3d_data%TH3D_STOP).OR.
1037 . (manim>=4.AND.manim<=15).OR.h3d_data%MH3D/=0))
1038 . .OR.h3d_data%N_VECT_PCONT_MAX>0)THEN
1039 IF (inconv==1) THEN
1040#include "lockon.inc"
1041 DO i=1,jlt
1042 IF(pene(i) == zero)cycle
1043 fncont(1,ix1(i)) =fncont(1,ix1(i)) + fxi(i)*h1(i)
1044 fncont(2,ix1(i)) =fncont(2,ix1(i)) + fyi(i)*h1(i)
1045 fncont(3,ix1(i)) =fncont(3,ix1(i)) + fzi(i)*h1(i)
1046 fncont(1,ix2(i)) =fncont(1,ix2(i)) + fxi(i)*h2(i)
1047 fncont(2,ix2(i)) =fncont(2,ix2(i)) + fyi(i)*h2(i)
1048 fncont(3,ix2(i)) =fncont(3,ix2(i)) + fzi(i)*h2(i)
1049 fncont(1,ix3(i)) =fncont(1,ix3(i)) + fxi(i)*h3(i)
1050 fncont(2,ix3(i)) =fncont(2,ix3(i)) + fyi(i)*h3(i)
1051 fncont(3,ix3(i)) =fncont(3,ix3(i)) + fzi(i)*h3(i)
1052 fncont(1,ix4(i)) =fncont(1,ix4(i)) + fxi(i)*h4(i)
1053 fncont(2,ix4(i)) =fncont(2,ix4(i)) + fyi(i)*h4(i)
1054 fncont(3,ix4(i)) =fncont(3,ix4(i)) + fzi(i)*h4(i)
1055 jg = nsvg(i)
1056 IF(jg>0) THEN
1057C en SPMD : traitement a refaire apres reception noeud remote si JG < 0
1058 IF (jg >numnod) THEN
1059 ig = jg - numnod
1060 CALL i24for1_fic(3 ,irtse ,nsne ,is2se ,is2pt ,
1061 + ig ,fxi(i) ,fyi(i) ,fzi(i) ,fncont ,
1062 + -1 )
1063 ELSE
1064 fncont(1,jg)=fncont(1,jg)- fxi(i)
1065 fncont(2,jg)=fncont(2,jg)- fyi(i)
1066 fncont(3,jg)=fncont(3,jg)- fzi(i)
1067 END IF
1068 ELSE ! cas noeud remote en SPMD
1069 jg = -jg
1070 IF(isedge_fi(nin)%P(jg)==1)THEN
1071 CALL i24for1_ficr(3 ,irtse_fi(nin)%P(1,1) ,nsne ,is2se_fi(nin)%P(1,1) ,is2pt_fi(nin)%P(1) ,
1072 + jg ,fxi(i) ,fyi(i) ,fzi(i) ,fnconti(nin)%P(1,1) ,
1073 + -1 ,nin )
1074 ELSE
1075 fnconti(nin)%P(1,jg)=fnconti(nin)%P(1,jg)-fxi(i)
1076 fnconti(nin)%P(2,jg)=fnconti(nin)%P(2,jg)-fyi(i)
1077 fnconti(nin)%P(3,jg)=fnconti(nin)%P(3,jg)-fzi(i)
1078 ENDIF
1079 ENDIF
1080 ENDDO
1081#include "lockoff.inc"
1082 END IF !(INCONV==1) THEN
1083 ENDIF
1084C
1085C---------------------------------
1086C SORTIES TH PAR SOUS INTERFACE
1087C---------------------------------
1088 IF(nisub/=0)THEN
1089 DO jsub=1,nisub
1090 DO j=1,25
1091 fsavsub1(j,jsub)=zero
1092 END DO
1093 ENDDO
1094 DO i=1,jlt
1095 IF(pene(i) == zero)cycle
1096 nn = nsvg(i)
1097 IF(nn>0)THEN
1098 in=cn_loc(i)
1099 IF (msegtyp(ce_loc(i)) < 0) THEN
1100 ie= - msegtyp(ce_loc(i))
1101 ELSE
1102 ie = ce_loc(i)
1103 ENDIF
1104 IF(ie > nrtm) ie=ie-nrtm
1105 jj =addsubs(in)
1106 kk =addsubm(ie)
1107 DO WHILE(jj<addsubs(in+1))
1108
1109 jsub=lisubs(jj)
1110 itypsub = typsub(jsub)
1111
1112 IF(itypsub == 1 ) THEN ! Defining specific inter
1113 iss1 = bitget(inflg_subs(jj),0)
1114 iss2 = bitget(inflg_subs(jj),1)
1115 igrn = bitget(inflg_subs(jj),2)
1116 ksub=lisubm(kk)
1117 DO WHILE((ksub<=jsub).AND.(kk<addsubm(ie+1)))
1118 ims1 = bitget(inflg_subm(kk),0)
1119 ims2 = bitget(inflg_subm(kk),1)
1120 IF(ksub==jsub)THEN
1121 IF(.NOT.(((ims1 == 1 .OR. ims2 == 1) .AND. igrn == 1).OR.
1122 . (ims1 == 1 .AND. iss2 == 1).OR.
1123 . (ims2 == 1 .AND. iss1 == 1))) THEN
1124 kk=kk+1
1125 ksub=lisubm(kk)
1126 cycle
1127 END IF
1128C
1129 impx=fxi(i)*dt12
1130 impy=fyi(i)*dt12
1131 impz=fzi(i)*dt12
1132 IF(ims2 > 0)THEN
1133 fsavsub1(1,jsub)=fsavsub1(1,jsub)-impx
1134 fsavsub1(2,jsub)=fsavsub1(2,jsub)-impy
1135 fsavsub1(3,jsub)=fsavsub1(3,jsub)-impz
1136 fsavsub1(11,jsub)=fsavsub1(11,jsub)-fni(i)*dt12
1137 ELSE
1138 fsavsub1(1,jsub)=fsavsub1(1,jsub)+impx
1139 fsavsub1(2,jsub)=fsavsub1(2,jsub)+impy
1140 fsavsub1(3,jsub)=fsavsub1(3,jsub)+impz
1141 fsavsub1(11,jsub)=fsavsub1(11,jsub)+fni(i)*dt12
1142 END IF
1143C
1144 IF(isensint(jsub+1)/=0) THEN
1145 IF(ims2 > 0)THEN
1146 fsavparit(jsub+1,1,i) = -fxi(i)
1147 fsavparit(jsub+1,2,i) = -fyi(i)
1148 fsavparit(jsub+1,3,i) = -fzi(i)
1149 ELSE
1150 fsavparit(jsub+1,1,i) = fxi(i)
1151 fsavparit(jsub+1,2,i) = fyi(i)
1152 fsavparit(jsub+1,3,i) = fzi(i)
1153 END IF
1154 ENDIF
1155C
1156 fsavsub1(8,jsub) =fsavsub1(8,jsub) +abs(impx)
1157 fsavsub1(9,jsub) =fsavsub1(9,jsub) +abs(impy)
1158 fsavsub1(10,jsub)=fsavsub1(10,jsub)+abs(impz)
1159C
1160 fsavsub1(11,jsub)=fsavsub1(11,jsub)+fni(i)*dt12
1161
1162 IF(parameters%INTCAREA > 0) THEN
1163 IF(nn <=numnod) THEN
1164 fsavsub1(25,jsub) = fsavsub1(25,jsub) + parameters%INTAREAN(nn)
1165 ELSE
1166 ig = nn - numnod
1167 CALL i24intarea_fic(irtse ,nsne ,is2se ,is2pt ,ig ,
1168 + nrtse , numnod ,parameters%INTAREAN, arean_fic)
1169 fsavsub1(25,jsub) = fsavsub1(25,jsub) + arean_fic
1170 ENDIF
1171 ENDIF
1172C
1173 END IF
1174
1175 kk=kk+1
1176 ksub=lisubm(kk)
1177 ENDDO
1178 jj=jj+1
1179
1180 ELSEIF(itypsub == 2 ) THEN ! Inter =0 : collecting forces from all inter with only 1 surf : secondary side
1181
1182 impx=fxi(i)*dt12
1183 impy=fyi(i)*dt12
1184 impz=fzi(i)*dt12
1185
1186 fsavsub1(1,jsub)=fsavsub1(1,jsub)+impx
1187 fsavsub1(2,jsub)=fsavsub1(2,jsub)+impy
1188 fsavsub1(3,jsub)=fsavsub1(3,jsub)+impz
1189
1190 fsavsub1(8,jsub) =fsavsub1(8,jsub) +abs(impx)
1191 fsavsub1(9,jsub) =fsavsub1(9,jsub) +abs(impy)
1192 fsavsub1(10,jsub)=fsavsub1(10,jsub)+abs(impz)
1193
1194 fsavsub1(11,jsub)=fsavsub1(11,jsub)+fni(i)*dt12
1195
1196 IF(isensint(jsub+1)/=0) THEN
1197 fsavparit(jsub+1,1,i+nft) = fxi(i)
1198 fsavparit(jsub+1,2,i+nft) = fyi(i)
1199 fsavparit(jsub+1,3,i+nft) = fzi(i)
1200 ENDIF
1201
1202 IF(parameters%INTCAREA > 0) THEN
1203 IF(nn <=numnod) THEN
1204 fsavsub1(25,jsub) = fsavsub1(25,jsub) + parameters%INTAREAN(nn)
1205 ELSE
1206 ig = nn - numnod
1207 CALL i24intarea_fic(irtse ,nsne ,is2se ,is2pt ,ig ,
1208 + nrtse , numnod , parameters%INTAREAN, arean_fic)
1209 fsavsub1(25,jsub) = fsavsub1(25,jsub) + arean_fic
1210 ENDIF
1211 ENDIF
1212
1213 jj=jj+1
1214
1215 ELSEIF(itypsub == 3) THEN !Inter =0 : collecting forces from all inter with 2 surfs
1216
1217
1218
1219 iss2 = bitget(inflg_subs(jj),0)
1220 iss1 = bitget(inflg_subs(jj),1)
1221 ksub=lisubm(kk)
1222
1223 DO WHILE((ksub<=jsub).AND.(kk<addsubm(ie+1)))
1224 ims2 = bitget(inflg_subm(kk),0)
1225 ims1 = bitget(inflg_subm(kk),1)
1226
1227 IF(ksub==jsub)THEN
1228 IF(.NOT.((ims1 == 1 .AND. iss2 == 1).OR.
1229 . (ims2 == 1 .AND. iss1 == 1))) THEN
1230 kk=kk+1
1231 ksub=lisubm(kk)
1232 cycle
1233 END IF
1234
1235 impx=fxi(i)*dt12
1236 impy=fyi(i)*dt12
1237 impz=fzi(i)*dt12
1238
1239 IF(ims2 > 0)THEN
1240 fsavsub1(1,jsub)=fsavsub1(1,jsub)-impx
1241 fsavsub1(2,jsub)=fsavsub1(2,jsub)-impy
1242 fsavsub1(3,jsub)=fsavsub1(3,jsub)-impz
1243 fsavsub1(11,jsub)=fsavsub1(11,jsub)-fni(i)*dt12
1244 ELSE
1245 fsavsub1(1,jsub)=fsavsub1(1,jsub)+impx
1246 fsavsub1(2,jsub)=fsavsub1(2,jsub)+impy
1247 fsavsub1(3,jsub)=fsavsub1(3,jsub)+impz
1248 fsavsub1(11,jsub)=fsavsub1(11,jsub)+fni(i)*dt12
1249 END IF
1250C
1251 IF(isensint(jsub+1)/=0) THEN
1252 IF(ims2 > 0)THEN
1253 fsavparit(jsub+1,1,i+nft) = -fxi(i)
1254 fsavparit(jsub+1,2,i+nft) = -fyi(i)
1255 fsavparit(jsub+1,3,i+nft) = -fzi(i)
1256 ELSE
1257 fsavparit(jsub+1,1,i+nft) = fxi(i)
1258 fsavparit(jsub+1,2,i+nft) = fyi(i)
1259 fsavparit(jsub+1,3,i+nft) = fzi(i)
1260 END IF
1261 ENDIF
1262C
1263 fsavsub1(8,jsub) =fsavsub1(8,jsub) +abs(impx)
1264 fsavsub1(9,jsub) =fsavsub1(9,jsub) +abs(impy)
1265 fsavsub1(10,jsub)=fsavsub1(10,jsub)+abs(impz)
1266C
1267 IF(parameters%INTCAREA > 0) THEN
1268 IF(nn <=numnod) THEN
1269 fsavsub1(25,jsub) = fsavsub1(25,jsub) + parameters%INTAREAN(nn)
1270 ELSE
1271 ig = nn - numnod
1272 CALL i24intarea_fic(irtse ,nsne ,is2se ,is2pt ,ig ,
1273 + nrtse , numnod , parameters%INTAREAN, arean_fic)
1274 fsavsub1(25,jsub) = fsavsub1(25,jsub) + arean_fic
1275 ENDIF
1276 ENDIF
1277 END IF
1278
1279 kk=kk+1
1280 ksub=lisubm(kk)
1281 ENDDO
1282 jj=jj+1
1283 ENDIF
1284
1285 END DO
1286 END IF
1287
1288
1289 IF (msegtyp(ce_loc(i)) < 0) THEN
1290 ie= - msegtyp(ce_loc(i))
1291 ELSE
1292 ie = ce_loc(i)
1293 ENDIF
1294 IF(ie > nrtm) ie=ie-nrtm
1295
1296 ! Compiler issue workaround,
1297 ! Loop packed in subroutine to avoid optimization breaking the code.
1298 CALL i24_save_sub(numnod,mvsiz,nisub,s_addsubm,s_lisubm,s_typsub,nisubmax,i_stok,
1299 * ie,itypsub,nin,i,nn,nft,
1300 * addsubm,lisubm,typsub,
1301 * parameters%INTAREAN,parameters%INTCAREA,isensint,
1302 * fxi,fyi,fzi,fni,dt12,
1303 * fsavsub1,fsavparit ,nrtse,
1304 * irtse,nsne,is2se ,is2pt,nsnr)
1305
1306 END DO
1307
1308
1309 IF(nspmd>1) THEN
1310 DO i=1,jlt
1311 IF(pene(i) == zero)cycle
1312 nn = nsvg(i)
1313 IF(nn<0)THEN
1314 nn = -nn
1315 IF (msegtyp(ce_loc(i)) < 0) THEN
1316 ie= - msegtyp(ce_loc(i))
1317 ELSE
1318 ie = ce_loc(i)
1319 ENDIF
1320 jj =addsubsfi(nin)%P(nn)
1321 kk =addsubm(ie)
1322 DO WHILE(jj<addsubsfi(nin)%P(nn+1))
1323 jsub=lisubsfi(nin)%P(jj)
1324 itypsub = typsub(jsub)
1325
1326 IF(itypsub == 1 ) THEN ! Defining specific inter
1327
1328 iss1 = bitget(inflg_subsfi(nin)%P(jj),0)
1329 iss2 = bitget(inflg_subsfi(nin)%P(jj),1)
1330 igrn = bitget(inflg_subsfi(nin)%P(jj),2)
1331 ksub=lisubm(kk)
1332 DO WHILE((ksub<=jsub).AND.(kk<addsubm(ie+1)))
1333 ims1 = bitget(inflg_subm(kk),0)
1334 ims2 = bitget(inflg_subm(kk),1)
1335 IF(ksub==jsub)THEN
1336 IF(.NOT.(((ims1 == 1 .OR. ims2 == 1) .AND. igrn == 1).OR.
1337 . (ims1 == 1 .AND. iss2 == 1).OR.
1338 . (ims2 == 1 .AND. iss1 == 1))) THEN
1339 kk=kk+1
1340 ksub=lisubm(kk)
1341 cycle
1342 END IF
1343 impx=fxi(i)*dt12
1344 impy=fyi(i)*dt12
1345 impz=fzi(i)*dt12
1346C
1347 IF(ims2 > 0)THEN
1348 fsavsub1(1,jsub)=fsavsub1(1,jsub)-impx
1349 fsavsub1(2,jsub)=fsavsub1(2,jsub)-impy
1350 fsavsub1(3,jsub)=fsavsub1(3,jsub)-impz
1351 fsavsub1(11,jsub)=fsavsub1(11,jsub)-fni(i)*dt12
1352 ELSE
1353 fsavsub1(1,jsub)=fsavsub1(1,jsub)+impx
1354 fsavsub1(2,jsub)=fsavsub1(2,jsub)+impy
1355 fsavsub1(3,jsub)=fsavsub1(3,jsub)+impz
1356 fsavsub1(11,jsub)=fsavsub1(11,jsub)+fni(i)*dt12
1357 END IF
1358C
1359 IF(isensint(jsub+1)/=0) THEN
1360 IF(ims2 > 0)THEN
1361 fsavparit(jsub+1,1,i) = -fxi(i)
1362 fsavparit(jsub+1,2,i) = -fyi(i)
1363 fsavparit(jsub+1,3,i) = -fzi(i)
1364 ELSE
1365 fsavparit(jsub+1,1,i) = fxi(i)
1366 fsavparit(jsub+1,2,i) = fyi(i)
1367 fsavparit(jsub+1,3,i) = fzi(i)
1368 END IF
1369 ENDIF
1370C
1371 fsavsub1(8,jsub) =fsavsub1(8,jsub) +abs(impx)
1372 fsavsub1(9,jsub) =fsavsub1(9,jsub) +abs(impy)
1373 fsavsub1(10,jsub)=fsavsub1(10,jsub)+abs(impz)
1374C
1375 fsavsub1(11,jsub)=fsavsub1(11,jsub)+fni(i)*dt12
1376C
1377 IF(parameters%INTCAREA > 0) THEN
1378 IF(isedge_fi(nin)%P(nn)==1)THEN
1379 CALL i24intarea_fic(irtse_fi(nin)%P ,nsnr ,is2se_fi(nin)%P ,is2pt_fi(nin)%P ,nn ,
1380 + nsnr , nsnr , intareanfi(nin)%P, arean_fic)
1381 fsavsub1(25,jsub) = fsavsub1(25,jsub) + arean_fic
1382 ELSE ! cas noeud remote en SPMD non edge
1383 fsavsub1(25,jsub) = fsavsub1(25,jsub) + intareanfi(nin)%P(nn)
1384 ENDIF
1385 ENDIF
1386 END IF
1387
1388 kk=kk+1
1389 ksub=lisubm(kk)
1390 ENDDO
1391 jj=jj+1
1392
1393 ELSEIF(itypsub == 2 ) THEN ! Inter =0 : collecting forces from all inter with only 1 surf : secondary side
1394
1395 impx=fxi(i)*dt12
1396 impy=fyi(i)*dt12
1397 impz=fzi(i)*dt12
1398
1399 fsavsub1(1,jsub)=fsavsub1(1,jsub)+impx
1400 fsavsub1(2,jsub)=fsavsub1(2,jsub)+impy
1401 fsavsub1(3,jsub)=fsavsub1(3,jsub)+impz
1402
1403 fsavsub1(8,jsub) =fsavsub1(8,jsub) +abs(impx)
1404 fsavsub1(9,jsub) =fsavsub1(9,jsub) +abs(impy)
1405 fsavsub1(10,jsub)=fsavsub1(10,jsub)+abs(impz)
1406
1407 fsavsub1(11,jsub)=fsavsub1(11,jsub)+fni(i)*dt12
1408
1409 IF(isensint(jsub+1)/=0) THEN
1410 fsavparit(jsub+1,1,i+nft) = fxi(i)
1411 fsavparit(jsub+1,2,i+nft) = fyi(i)
1412 fsavparit(jsub+1,3,i+nft) = fzi(i)
1413 ENDIF
1414
1415 IF(parameters%INTCAREA > 0) THEN
1416 IF(isedge_fi(nin)%P(nn)==1)THEN
1417 CALL i24intarea_fic(irtse_fi(nin)%P ,nsnr ,is2se_fi(nin)%P ,is2pt_fi(nin)%P ,nn ,
1418 + nsnr , nsnr , intareanfi(nin)%P, arean_fic)
1419 fsavsub1(25,jsub) = fsavsub1(25,jsub) + arean_fic
1420 ELSE ! cas noeud remote en SPMD non edge
1421 fsavsub1(25,jsub) = fsavsub1(25,jsub) + intareanfi(nin)%P(nn)
1422 ENDIF
1423 ENDIF
1424
1425 jj=jj+1
1426
1427 ELSEIF(itypsub == 3 ) THEN ! Inter =0 : collecting forces from all inter with 2 surfaces
1428
1429 iss2 = bitget(inflg_subsfi(nin)%P(jj),0)
1430 iss1 = bitget(inflg_subsfi(nin)%P(jj),1)
1431 ksub=lisubm(kk)
1432 DO WHILE((ksub<=jsub).AND.(kk<addsubm(ie+1)))
1433 ims2 = bitget(inflg_subm(kk),0)
1434 ims1 = bitget(inflg_subm(kk),1)
1435 IF(ksub==jsub)THEN
1436 IF(.NOT.((ims1 == 1 .AND. iss2 == 1).OR.
1437 . (ims2 == 1 .AND. iss1 == 1))) THEN
1438 kk=kk+1
1439 ksub=lisubm(kk)
1440 cycle
1441 END IF
1442
1443 impx=fxi(i)*dt12
1444 impy=fyi(i)*dt12
1445 impz=fzi(i)*dt12
1446
1447 IF(ims2 > 0)THEN
1448 fsavsub1(1,jsub)=fsavsub1(1,jsub)-impx
1449 fsavsub1(2,jsub)=fsavsub1(2,jsub)-impy
1450 fsavsub1(3,jsub)=fsavsub1(3,jsub)-impz
1451 fsavsub1(11,jsub)=fsavsub1(11,jsub)-fni(i)*dt12
1452 ELSE
1453 fsavsub1(1,jsub)=fsavsub1(1,jsub)+impx
1454 fsavsub1(2,jsub)=fsavsub1(2,jsub)+impy
1455 fsavsub1(3,jsub)=fsavsub1(3,jsub)+impz
1456 fsavsub1(11,jsub)=fsavsub1(11,jsub)+fni(i)*dt12
1457 END IF
1458C
1459 IF(isensint(jsub+1)/=0) THEN
1460 IF(ims2 > 0)THEN
1461 fsavparit(jsub+1,1,i+nft) = -fxi(i)
1462 fsavparit(jsub+1,2,i+nft) = -fyi(i)
1463 fsavparit(jsub+1,3,i+nft) = -fzi(i)
1464 ELSE
1465 fsavparit(jsub+1,1,i+nft) = fxi(i)
1466 fsavparit(jsub+1,2,i+nft) = fyi(i)
1467 fsavparit(jsub+1,3,i+nft) = fzi(i)
1468 END IF
1469 ENDIF
1470
1471 fsavsub1(8,jsub) =fsavsub1(8,jsub) +abs(impx)
1472 fsavsub1(9,jsub) =fsavsub1(9,jsub) +abs(impy)
1473 fsavsub1(10,jsub)=fsavsub1(10,jsub)+abs(impz)
1474
1475 IF(parameters%INTCAREA > 0) THEN
1476 IF(isedge_fi(nin)%P(nn)==1)THEN
1477 CALL i24intarea_fic(irtse_fi(nin)%P ,nsnr ,is2se_fi(nin)%P ,is2pt_fi(nin)%P ,nn ,
1478 + nsnr , nsnr , intareanfi(nin)%P, arean_fic)
1479 fsavsub1(25,jsub) = fsavsub1(25,jsub) + arean_fic
1480 ELSE ! cas noeud remote en SPMD non edge
1481 fsavsub1(25,jsub) = fsavsub1(25,jsub) + intareanfi(nin)%P(nn)
1482 ENDIF
1483 ENDIF
1484
1485 ENDIF
1486 kk=kk+1
1487 ksub=lisubm(kk)
1488 ENDDO
1489 jj=jj+1
1490
1491 ENDIF
1492
1493 END DO
1494 END IF
1495 END DO
1496 END IF
1497 END IF
1498
1499C---------------------------------
1500C NEW FRICTION MODELS
1501C---------------------------------
1502
1503C Friction coefficient computation
1504C++++++++++++++++++++++++++++++++++++++++
1505
1506 IF(iorthfric == 0) THEN
1507C++ Isotropic Friction
1508
1509 IF (mfrot==0) THEN
1510C--- Coulomb friction
1511 DO i=1,jlt
1512 IF(pene(i) == 0) THEN
1513 xmu(i) = zero
1514 cycle
1515 ENDIF
1516 xmu(i) = fricc(i)
1517 ENDDO
1518 ELSEIF (mfrot==1) THEN
1519C--- Viscous friction
1520 DO i=1,jlt
1521 IF(pene(i) == 0) THEN
1522 xmu(i) = zero
1523 cycle
1524 ENDIF
1525 ie=ce_loc(i)
1526 aa = n1(i)*vx(i) + n2(i)*vy(i) + n3(i)*vz(i)
1527 v2 = (vx(i) - n1(i)*aa)**2
1528 . + (vy(i) - n2(i)*aa)**2
1529 . + (vz(i) - n3(i)*aa)**2
1530 vv = sqrt(max(em30,v2))
1531 ax1 = x(1,irect(3,ie)) - x(1,irect(1,ie))
1532 ay1 = x(2,irect(3,ie)) - x(2,irect(1,ie))
1533 az1 = x(3,irect(3,ie)) - x(3,irect(1,ie))
1534 ax2 = x(1,irect(4,ie)) - x(1,irect(2,ie))
1535 ay2 = x(2,irect(4,ie)) - x(2,irect(2,ie))
1536 az2 = x(3,irect(4,ie)) - x(3,irect(2,ie))
1537 ax = ay1*az2 - az1*ay2
1538 ay = az1*ax2 - ax1*az2
1539 az = ax1*ay2 - ay1*ax2
1540 area = half*sqrt(ax*ax+ay*ay+az*az)
1541 p = -fni(i)/area
1542 xmu(i) = fricc(i) + (fric_coefs(i,1) + fric_coefs(i,4)*p ) * p
1543 . +(fric_coefs(i,2) + fric_coefs(i,3)*p) * vv + fric_coefs(i,5)*v2
1544 xmu(i) = max(xmu(i),em30)
1545 ENDDO
1546 ELSEIF(mfrot==2)THEN
1547C--- Darmstad LAW
1548 DO i=1,jlt
1549 IF(pene(i) == 0) THEN
1550 xmu(i) = zero
1551 cycle
1552 ENDIF
1553 ie=ce_loc(i)
1554 aa = n1(i)*vx(i) + n2(i)*vy(i) + n3(i)*vz(i)
1555 v2 = (vx(i) - n1(i)*aa)**2
1556 . + (vy(i) - n2(i)*aa)**2
1557 . + (vz(i) - n3(i)*aa)**2
1558 vv = sqrt(max(em30,v2))
1559 ax1 = x(1,irect(3,ie)) - x(1,irect(1,ie))
1560 ay1 = x(2,irect(3,ie)) - x(2,irect(1,ie))
1561 az1 = x(3,irect(3,ie)) - x(3,irect(1,ie))
1562 ax2 = x(1,irect(4,ie)) - x(1,irect(2,ie))
1563 ay2 = x(2,irect(4,ie)) - x(2,irect(2,ie))
1564 az2 = x(3,irect(4,ie)) - x(3,irect(2,ie))
1565 ax = ay1*az2 - az1*ay2
1566 ay = az1*ax2 - ax1*az2
1567 az = ax1*ay2 - ay1*ax2
1568 area = half*sqrt(ax*ax+ay*ay+az*az)
1569 p = -fni(i)/area
1570 xmu(i) = fricc(i)
1571 . + fric_coefs(i,1)*exp(fric_coefs(i,2)*vv)*p*p
1572 . + fric_coefs(i,3)*exp(fric_coefs(i,4)*vv)*p
1573 . + fric_coefs(i,5)*exp(fric_coefs(i,6)*vv)
1574 xmu(i) = max(xmu(i),em30)
1575 ENDDO
1576 ELSEIF (mfrot==3) THEN
1577C--- Renard LAW
1578 DO i=1,jlt
1579 IF(pene(i) == 0) THEN
1580 xmu(i) = zero
1581 cycle
1582 ENDIF
1583 aa = n1(i)*vx(i) + n2(i)*vy(i) + n3(i)*vz(i)
1584 v2 = (vx(i) - n1(i)*aa)**2
1585 . + (vy(i) - n2(i)*aa)**2
1586 . + (vz(i) - n3(i)*aa)**2
1587 vv = sqrt(max(em30,v2))
1588 IF(vv>=0.AND.vv<=fric_coefs(i,5)) THEN
1589 dmu = fric_coefs(i,3)-fric_coefs(i,1)
1590 vv1 = vv / fric_coefs(i,5)
1591 xmu(i) = fric_coefs(i,1)+ dmu*vv1*(two-vv1)
1592 ELSEIF(vv>fric_coefs(i,5).AND.vv<fric_coefs(i,6)) THEN
1593 dmu = fric_coefs(i,4)-fric_coefs(i,3)
1594 vv1 = (vv - fric_coefs(i,5))/(fric_coefs(i,6)-fric_coefs(i,5))
1595 xmu(i) = fric_coefs(i,3)+ dmu * (three-two*vv1)*vv1**2
1596 ELSE
1597 dmu = fric_coefs(i,2)-fric_coefs(i,4)
1598 vv2 = (vv - fric_coefs(i,6))**2
1599 xmu(i) = fric_coefs(i,2) - dmu / (one + dmu*vv2)
1600 ENDIF
1601 xmu(i) = max(xmu(i),em30)
1602 ENDDO
1603 ELSEIF(mfrot==4)THEN
1604C--- Exponential decay model
1605 DO i=1,jlt
1606 aa = n1(i)*vx(i) + n2(i)*vy(i) + n3(i)*vz(i)
1607 v2 = (vx(i) - n1(i)*aa)**2
1608 . + (vy(i) - n2(i)*aa)**2
1609 . + (vz(i) - n3(i)*aa)**2
1610 vv = sqrt(max(em30,v2))
1611 xmu(i) = fric_coefs(i,1)
1612 . + (fricc(i)-fric_coefs(i,1))*exp(-fric_coefs(i,2)*vv)
1613 xmu(i) = max(xmu(i),em30)
1614 ENDDO
1615 ENDIF
1616
1617 ELSE
1618C++ Orthotropic Friction
1619
1620 IF (mfrot==0) THEN
1621C--- Coulomb friction
1622#include "vectorize.inc"
1623 DO k=1,nfisot
1624 i = indexisot(k)
1625 xmu(i) = fricc(i)
1626 ENDDO
1627
1628#include "vectorize.inc"
1629 DO k=1,nforth
1630 i = indexorth(k)
1631 xmu(i) = fricc(i)
1632 xmu2(i) = fricc2(i)
1633 IF(xmu(i)<=em30) THEN
1634 xmu(i) = em30
1635 dir1(i,1:3) = zero
1636 an(k) = zero
1637 ELSE
1638 an(k)=xmu(i)**2
1639 an(k)=one/an(k)
1640 ENDIF
1641 IF(xmu2(i)<=em30) THEN
1642 xmu2(i) = em30
1643 dir2(i,1:3) = zero
1644 bn(k) = zero
1645 ELSE
1646 bn(k)=xmu2(i)**2
1647 bn(k)=one/bn(k)
1648 ENDIF
1649 ENDDO
1650 ELSEIF (mfrot==1) THEN
1651C--- Viscous friction
1652#include "vectorize.inc"
1653 DO k=1,nfisot ! isotropic friction couples
1654 i = indexisot(k)
1655 aa = n1(i)*vx(i) + n2(i)*vy(i) + n3(i)*vz(i)
1656 v2 = (vx(i) - n1(i)*aa)**2
1657 . + (vy(i) - n2(i)*aa)**2
1658 . + (vz(i) - n3(i)*aa)**2
1659 vv = sqrt(max(em30,v2))
1660
1661 ax1 = x(1,irect(3,ie)) - x(1,irect(1,ie))
1662 ay1 = x(2,irect(3,ie)) - x(2,irect(1,ie))
1663 az1 = x(3,irect(3,ie)) - x(3,irect(1,ie))
1664 ax2 = x(1,irect(4,ie)) - x(1,irect(2,ie))
1665 ay2 = x(2,irect(4,ie)) - x(2,irect(2,ie))
1666 az2 = x(3,irect(4,ie)) - x(3,irect(2,ie))
1667 ax = ay1*az2 - az1*ay2
1668 ay = az1*ax2 - ax1*az2
1669 az = ax1*ay2 - ay1*ax2
1670 area = half*sqrt(ax*ax+ay*ay+az*az)
1671 p = -fni(i)/area
1672 xmu(i) = fricc(i) + (fric_coefs(i,1) + fric_coefs(i,4)*p ) * p
1673 . +(fric_coefs(i,2) + fric_coefs(i,3)*p) * vv + fric_coefs(i,5)*v2
1674 xmu(i) = max(xmu(i),em30)
1675 ENDDO
1676
1677#include "vectorize.inc"
1678 DO k=1,nforth ! Orthotropic friction couples
1679 i = indexorth(k)
1680 ie=ce_loc(i)
1681c
1682 ax1 = x(1,irect(3,ie)) - x(1,irect(1,ie))
1683 ay1 = x(2,irect(3,ie)) - x(2,irect(1,ie))
1684 az1 = x(3,irect(3,ie)) - x(3,irect(1,ie))
1685 ax2 = x(1,irect(4,ie)) - x(1,irect(2,ie))
1686 ay2 = x(2,irect(4,ie)) - x(2,irect(2,ie))
1687 az2 = x(3,irect(4,ie)) - x(3,irect(2,ie))
1688 ax = ay1*az2 - az1*ay2
1689 ay = az1*ax2 - ax1*az2
1690 az = ax1*ay2 - ay1*ax2
1691 area = half*sqrt(ax*ax+ay*ay+az*az)
1692 p = -fni(i)/area
1693c
1694 v2 = vx(i)*dir1(k,1) +vy(i)*dir1(k,2)+vz(i)*dir1(k,3)
1695 vv = max(em30,v2)
1696 xmu(i) = fricc(i) + (fric_coefs(i,1) + fric_coefs(i,4)*p ) * p
1697 . +(fric_coefs(i,2) + fric_coefs(i,3)*p) * vv + fric_coefs(i,5)*vv*vv
1698
1699 v2 = vx(i)*dir2(k,1) +vy(i)*dir2(k,2)+vz(i)*dir2(k,3)
1700 vv = max(em30,v2)
1701 xmu2(i) = fricc2(i) + (fric_coefs2(i,1) + fric_coefs2(i,4)*p ) * p
1702 . +(fric_coefs2(i,2) + fric_coefs2(i,3)*p) * vv + fric_coefs2(i,5)*vv*vv
1703
1704 xmu(i) = max(xmu(i),em30)
1705 xmu2(i) = max(xmu2(i),em30)
1706
1707 ENDDO
1708
1709#include "vectorize.inc"
1710 DO k=1,nforth ! Orthotropic friction couples : treat differently cases where mu1=0 AND/OR mu2=0
1711 i = indexorth(k)
1712 IF(xmu(i)<=em30) THEN
1713 xmu(i) = em30
1714 dir1(i,1:3) = zero
1715 an(k) = zero
1716 ELSE
1717 an(k)=xmu(i)**2
1718 an(k)=one/an(k)
1719 ENDIF
1720 IF(xmu2(i)<=em30) THEN
1721 xmu2(i) = em30
1722 dir2(i,1:3) = zero
1723 bn(k) = zero
1724 ELSE
1725 bn(k)=xmu2(i)**2
1726 bn(k)=one/bn(k)
1727 ENDIF
1728 ENDDO
1729
1730
1731 ELSEIF(mfrot==2)THEN
1732C--- Loi Darmstad
1733#include "vectorize.inc"
1734 DO k=1,nfisot ! isotropic friction couples
1735 i = indexisot(k)
1736 ie=ce_loc(i)
1737 aa = n1(i)*vx(i) + n2(i)*vy(i) + n3(i)*vz(i)
1738 v2 = (vx(i) - n1(i)*aa)**2
1739 . + (vy(i) - n2(i)*aa)**2
1740 . + (vz(i) - n3(i)*aa)**2
1741 vv = sqrt(max(em30,v2))
1742 ax1 = x(1,irect(3,ie)) - x(1,irect(1,ie))
1743 ay1 = x(2,irect(3,ie)) - x(2,irect(1,ie))
1744 az1 = x(3,irect(3,ie)) - x(3,irect(1,ie))
1745 ax2 = x(1,irect(4,ie)) - x(1,irect(2,ie))
1746 ay2 = x(2,irect(4,ie)) - x(2,irect(2,ie))
1747 az2 = x(3,irect(4,ie)) - x(3,irect(2,ie))
1748 ax = ay1*az2 - az1*ay2
1749 ay = az1*ax2 - ax1*az2
1750 az = ax1*ay2 - ay1*ax2
1751 area = half*sqrt(ax*ax+ay*ay+az*az)
1752 p = -fni(i)/area
1753 xmu(i) = fricc(i)
1754 . + fric_coefs(i,1)*exp(fric_coefs(i,2)*vv)*p*p
1755 . + fric_coefs(i,3)*exp(fric_coefs(i,4)*vv)*p
1756 . + fric_coefs(i,5)*exp(fric_coefs(i,6)*vv)
1757 xmu(i) = max(xmu(i),em30)
1758 ENDDO
1759c
1760#include "vectorize.inc"
1761 DO k=1,nforth ! Orthotropic friction couples
1762 i = indexorth(k)
1763 ie=ce_loc(i)
1764c
1765 ax1 = x(1,irect(3,ie)) - x(1,irect(1,ie))
1766 ay1 = x(2,irect(3,ie)) - x(2,irect(1,ie))
1767 az1 = x(3,irect(3,ie)) - x(3,irect(1,ie))
1768 ax2 = x(1,irect(4,ie)) - x(1,irect(2,ie))
1769 ay2 = x(2,irect(4,ie)) - x(2,irect(2,ie))
1770 az2 = x(3,irect(4,ie)) - x(3,irect(2,ie))
1771 ax = ay1*az2 - az1*ay2
1772 ay = az1*ax2 - ax1*az2
1773 az = ax1*ay2 - ay1*ax2
1774 area = half*sqrt(ax*ax+ay*ay+az*az)
1775 p = -fni(i)/area
1776c
1777 v2 = vx(i)*dir1(k,1) +vy(i)*dir1(k,2)+vz(i)*dir1(k,3)
1778 vv = max(em30,v2)
1779 xmu(i) = fricc(i)
1780 . + fric_coefs(i,1)*exp(fric_coefs(i,2)*vv)*p*p
1781 . + fric_coefs(i,3)*exp(fric_coefs(i,4)*vv)*p
1782 . + fric_coefs(i,5)*exp(fric_coefs(i,6)*vv)
1783c
1784 v2 = vx(i)*dir2(k,1) +vy(i)*dir2(k,2)+vz(i)*dir2(k,3)
1785 vv = max(em30,v2)
1786 xmu2(i) = fricc2(i)
1787 . + fric_coefs2(i,1)*exp(fric_coefs2(i,2)*vv)*p*p
1788 . + fric_coefs2(i,3)*exp(fric_coefs2(i,4)*vv)*p
1789 . + fric_coefs2(i,5)*exp(fric_coefs2(i,6)*vv)
1790
1791 xmu(i) = max(xmu(i),em30)
1792 xmu2(i) = max(xmu2(i),em30)
1793
1794 ENDDO
1795
1796#include "vectorize.inc"
1797 DO k=1,nforth ! Orthotropic friction couples : treat differently cases where mu1=0 AND/OR mu2=0
1798 i = indexorth(k)
1799 IF(xmu(i)<=em30) THEN
1800 xmu(i) = em30
1801 dir1(i,1:3) = zero
1802 an(k) = zero
1803 ELSE
1804 an(k)=xmu(i)**2
1805 an(k)=one/an(k)
1806 ENDIF
1807 IF(xmu2(i)<=em30) THEN
1808 xmu2(i) = em30
1809 dir2(i,1:3) = zero
1810 bn(k) = zero
1811 ELSE
1812 bn(k)=xmu2(i)**2
1813 bn(k)=one/bn(k)
1814 ENDIF
1815 ENDDO
1816
1817 ELSEIF (mfrot==3) THEN
1818C--- Loi Renard
1819#include "vectorize.inc"
1820 DO k=1,nfisot ! isotropic friction couples
1821 i = indexisot(k)
1822 aa = n1(i)*vx(i) + n2(i)*vy(i) + n3(i)*vz(i)
1823 v2 = (vx(i) - n1(i)*aa)**2
1824 . + (vy(i) - n2(i)*aa)**2
1825 . + (vz(i) - n3(i)*aa)**2
1826 vv = sqrt(max(em30,v2))
1827 IF(vv>=0.AND.vv<=fric_coefs(i,5)) THEN
1828 dmu = fric_coefs(i,3)-fric_coefs(i,1)
1829 vv1 = vv / fric_coefs(i,5)
1830 xmu(i) = fric_coefs(i,1)+ dmu*vv1*(two-vv1)
1831 ELSEIF(vv>fric_coefs(i,5).AND.vv<fric_coefs(i,6)) THEN
1832 dmu = fric_coefs(i,4)-fric_coefs(i,3)
1833 vv1 = (vv - fric_coefs(i,5))/(fric_coefs(i,6)-fric_coefs(i,5))
1834 xmu(i) = fric_coefs(i,3)+ dmu * (three-two*vv1)*vv1**2
1835 ELSE
1836 dmu = fric_coefs(i,2)-fric_coefs(i,4)
1837 vv2 = (vv - fric_coefs(i,6))**2
1838 xmu(i) = fric_coefs(i,2) - dmu / (one + dmu*vv2)
1839 ENDIF
1840 xmu(i) = max(xmu(i),em30)
1841 ENDDO
1842
1843#include "vectorize.inc"
1844 DO k=1,nforth ! Orthotropic friction couples
1845 i = indexorth(k)
1846c
1847 v2 = vx(i)*dir1(k,1) +vy(i)*dir1(k,2)+vz(i)*dir1(k,3)
1848 vv = max(em30,v2)
1849 IF(vv>=0.AND.vv<=fric_coefs(i,5)) THEN
1850 dmu = fric_coefs(i,3)-fric_coefs(i,1)
1851 vv1 = vv / fric_coefs(i,5)
1852 xmu(i) = fric_coefs(i,1)+ dmu*vv1*(two-vv1)
1853 ELSEIF(vv>fric_coefs(i,5).AND.vv<fric_coefs(i,6)) THEN
1854 dmu = fric_coefs(i,4)-fric_coefs(i,3)
1855 vv1 = (vv - fric_coefs(i,5))/(fric_coefs(i,6)-fric_coefs(i,5))
1856 xmu(i) = fric_coefs(i,3)+ dmu * (three-two*vv1)*vv1**2
1857 ELSE
1858 dmu = fric_coefs(i,2)-fric_coefs(i,4)
1859 vv2 = (vv - fric_coefs(i,6))**2
1860 xmu(i) = fric_coefs(i,2) - dmu / (one + dmu*vv2)
1861 ENDIF
1862
1863 v2 = vx(i)*dir2(k,1) +vy(i)*dir2(k,2)+vz(i)*dir2(k,3)
1864 vv = max(em30,v2)
1865 IF(vv>=0.AND.vv<=fric_coefs2(i,5)) THEN
1866 dmu = fric_coefs2(i,3)-fric_coefs2(i,1)
1867 vv1 = vv / fric_coefs2(i,5)
1868 xmu2(i) = fric_coefs2(i,1)+ dmu*vv1*(two-vv1)
1869 ELSEIF(vv>fric_coefs2(i,5).AND.vv<fric_coefs2(i,6)) THEN
1870 dmu = fric_coefs2(i,4)-fric_coefs2(i,3)
1871 vv1 = (vv - fric_coefs2(i,5))/(fric_coefs2(i,6)-fric_coefs2(i,5))
1872 xmu2(i) = fric_coefs2(i,3)+ dmu * (three-two*vv1)*vv1**2
1873 ELSE
1874 dmu = fric_coefs2(i,2)-fric_coefs2(i,4)
1875 vv2 = (vv - fric_coefs2(i,6))**2
1876 xmu2(i) = fric_coefs2(i,2) - dmu / (one + dmu*vv2)
1877 ENDIF
1878 xmu(i) = max(xmu(i),em30)
1879 xmu2(i) = max(xmu2(i),em30)
1880 ENDDO
1881
1882#include "vectorize.inc"
1883 DO k=1,nforth ! Orthotropic friction couples : treat differently cases where mu1=0 AND/OR mu2=0
1884 i = indexorth(k)
1885 IF(xmu(i)<=em30) THEN
1886 xmu(i) = em30
1887 dir1(i,1:3) = zero
1888 an(k) = zero
1889 ELSE
1890 an(k)=xmu(i)**2
1891 an(k)=one/an(k)
1892 ENDIF
1893 IF(xmu2(i)<=em30) THEN
1894 xmu2(i) = em30
1895 dir2(i,1:3) = zero
1896 bn(k) = zero
1897 ELSE
1898 bn(k)=xmu2(i)**2
1899 bn(k)=one/bn(k)
1900 ENDIF
1901 ENDDO
1902
1903
1904 ELSEIF(mfrot==4)THEN
1905C--- Exponential decay model
1906#include "vectorize.inc"
1907 DO k=1,nfisot ! isotropic friction couples
1908 i = indexisot(k)
1909 IF(pene(i) == 0) THEN
1910 xmu(i) = zero
1911 xmu2(i) = zero
1912 cycle
1913 ENDIF
1914 aa = n1(i)*vx(i) + n2(i)*vy(i) + n3(i)*vz(i)
1915 v2 = (vx(i) - n1(i)*aa)**2
1916 . + (vy(i) - n2(i)*aa)**2
1917 . + (vz(i) - n3(i)*aa)**2
1918 vv = sqrt(max(em30,v2))
1919 xmu(i) = fric_coefs(i,1)
1920 . + (fricc(i)-fric_coefs(i,1))*exp(-fric_coefs(i,2)*vv)
1921 xmu(i) = max(xmu(i),em30)
1922 ENDDO
1923c
1924#include "vectorize.inc"
1925 DO k=1,nforth ! Orthotropic friction couples
1926 i = indexorth(k)
1927 IF(pene(i) == 0) THEN
1928 xmu(i) = zero
1929 xmu2(i) = zero
1930 cycle
1931 ENDIF
1932c
1933 v2 = vx(i)*dir1(i,1) +vy(i)*dir1(i,2)+vz(i)*dir1(i,3)
1934 vv = max(em30,v2)
1935 xmu(i) = fricc(i)
1936 . + (fric_coefs(i,1)-fricc(i))*exp(-fric_coefs(i,2)*vv)
1937c
1938 v2 = vx(i)*dir2(i,1) +vy(i)*dir2(i,2)+vz(i)*dir2(i,3)
1939 vv = max(em30,v2)
1940 xmu2(i) = fric_coefs2(i,1)
1941 . + (fricc2(i)-fric_coefs2(i,1))*exp(-fric_coefs2(i,2)*vv)
1942 xmu(i) = max(xmu(i),em30)
1943 xmu2(i) = max(xmu2(i),em30)
1944 ENDDO
1945
1946#include "vectorize.inc"
1947 DO k=1,nforth ! Orthotropic friction couples : treat differently cases where mu1=0 AND/OR mu2=0
1948 i = indexorth(k)
1949 IF(xmu(i)<=em30) THEN
1950 xmu(i) = em30
1951 dir1(i,1:3) = zero
1952 an(k) = zero
1953 ELSE
1954 an(k)=xmu(i)**2
1955 an(k)=one/an(k)
1956 ENDIF
1957 IF(xmu2(i)<=em30) THEN
1958 xmu2(i) = em30
1959 dir2(i,1:3) = zero
1960 bn(k) = zero
1961 ELSE
1962 bn(k)=xmu2(i)**2
1963 bn(k)=one/bn(k)
1964 ENDIF
1965 ENDDO
1966
1967 ENDIF
1968
1969 ENDIF ! IORTHFRIC
1970C------------------
1971C TANGENT FORCE CALCULATION
1972C------------------
1973 fsav4 = zero
1974 fsav5 = zero
1975 fsav6 = zero
1976 fsav12= zero
1977 fsav13= zero
1978 fsav14= zero
1979 fsav15= zero
1980 IF (ifq /= 0) THEN
1981C---------------------------------
1982C INCREMENTAL (STIFFNESS) FORMULATION
1983C---------------------------------
1984 IF (ifq==13) THEN
1985 alpha = max(one,alpha0*dt12)
1986 ELSE
1987 alpha = alpha0
1988 ENDIF
1989
1990 IF(intnitsche > 0 ) THEN
1991
1992 DO i=1,jlt
1993 IF(pene(i) == zero)cycle
1994 ftn = forneqsi(i,1)*n1(i) + forneqsi(i,2)*n2(i) + forneqsi(i,3)*n3(i)
1995 fnitxt(i) = half*(forneqsi(i,1) - ftn*n1(i))
1996 fnityt(i) = half*(forneqsi(i,2) - ftn*n2(i))
1997 fnitzt(i) = half*(forneqsi(i,3) - ftn*n3(i))
1998 ENDDO
1999 ELSE
2000 DO i=1,jlt
2001 fnitxt(i) = zero
2002 fnityt(i) = zero
2003 fnitzt(i) = zero
2004 ENDDO
2005 ENDIF
2006
2007 IF(iorthfric ==0 ) THEN
2008C++ Isotropic Friction
2009
2010 IF (inconv==1) THEN
2011 DO i=1,jlt
2012 IF(pene(i) == zero)cycle
2013 fx = stif0(i)*vx(i)*dt12
2014 fy = stif0(i)*vy(i)*dt12
2015 fz = stif0(i)*vz(i)*dt12
2016 jg = nsvg(i)
2017 IF(jg>0) THEN
2018 n = cand_n_n(i)
2019c SECND_FR(4:6,N) is the old friction force
2020 fx = secnd_fr(4,n) + alpha*fx
2021 fy = secnd_fr(5,n) + alpha*fy
2022 fz = secnd_fr(6,n) + alpha*fz
2023 ftn = fx*n1(i) + fy*n2(i) + fz*n3(i)
2024 fx = fx - ftn*n1(i)
2025 fy = fy - ftn*n2(i)
2026 fz = fz - ftn*n3(i)
2027 fx = fx + fnitxt(i)
2028 fy = fy + fnityt(i)
2029 fz = fz + fnitzt(i)
2030 ft = fx*fx + fy*fy + fz*fz
2031 ft = max(ft,tiny)
2032 fn = fxi(i)**2+fyi(i)**2+fzi(i)**2
2033 beta = min(one,xmu(i)*sqrt(fn/ft))
2034 fxt(i) = fx * beta
2035 fyt(i) = fy * beta
2036 fzt(i) = fz * beta
2037
2038#include "lockon.inc"
2039 IF (abs(fxt(i)- fnitxt(i))>abs(secnd_fr(1,n))) secnd_fr(1,n) = fxt(i)- fnitxt(i)
2040 IF (abs(fyt(i)- fnityt(i))>abs(secnd_fr(2,n))) secnd_fr(2,n) = fyt(i)- fnityt(i)
2041 IF (abs(fzt(i)- fnitzt(i))>abs(secnd_fr(3,n))) secnd_fr(3,n) = fzt(i)- fnitzt(i)
2042C case equal abs but opposite sign
2043 IF ((fxt(i)- fnitxt(i))==-secnd_fr(1,n)) secnd_fr(1,n) = abs(fxt(i)- fnitxt(i))
2044 IF ((fyt(i)- fnityt(i))==-secnd_fr(2,n)) secnd_fr(2,n) = abs(fyt(i)- fnityt(i))
2045 IF ((fzt(i)- fnitzt(i))==-secnd_fr(3,n)) secnd_fr(3,n) = abs(fzt(i)- fnitzt(i))
2046
2047#include "lockoff.inc"
2048 ELSE ! cas noeud remote en SPMD
2049 jg = -jg
2050c SECND_FRFI(NIN)%P(4:6,JG) is the old friction force
2051 fx = secnd_frfi(nin)%P(4,jg) + alpha*fx
2052 fy = secnd_frfi(nin)%P(5,jg) + alpha*fy
2053 fz = secnd_frfi(nin)%P(6,jg) + alpha*fz
2054 ftn = fx*n1(i) + fy*n2(i) + fz*n3(i)
2055 fx = fx - ftn*n1(i)
2056 fy = fy - ftn*n2(i)
2057 fz = fz - ftn*n3(i)
2058 fx = fx + fnitxt(i)
2059 fy = fy + fnityt(i)
2060 fz = fz + fnitzt(i)
2061 ft = fx*fx + fy*fy + fz*fz
2062 ft = max(ft,tiny)
2063 fn = fxi(i)**2+fyi(i)**2+fzi(i)**2
2064 beta = min(one,xmu(i)*sqrt(fn/ft))
2065 fxt(i) = fx * beta
2066 fyt(i) = fy * beta
2067 fzt(i) = fz * beta
2068#include "lockon.inc"
2069 IF ( abs(fxt(i)- fnitxt(i)) > abs(secnd_frfi(nin)%P(1,jg)) )
2070 * secnd_frfi(nin)%P(1,jg) = fxt(i)- fnitxt(i)
2071 IF ( abs(fyt(i)- fnityt(i)) > abs(secnd_frfi(nin)%P(2,jg)) )
2072 * secnd_frfi(nin)%P(2,jg) = fyt(i)- fnityt(i)
2073 IF ( abs(fzt(i)- fnitzt(i)) > abs(secnd_frfi(nin)%P(3,jg)) )
2074 * secnd_frfi(nin)%P(3,jg) = fzt(i)- fnitzt(i)
2075C
2076 IF ((fxt(i)- fnitxt(i))== - secnd_frfi(nin)%P(1,jg) )
2077 * secnd_frfi(nin)%P(1,jg) = abs(fxt(i)- fnitxt(i))
2078 IF ((fyt(i)- fnityt(i))== - secnd_frfi(nin)%P(2,jg) )
2079 * secnd_frfi(nin)%P(2,jg) = abs(fyt(i)- fnityt(i))
2080 IF ((fzt(i)- fnitzt(i))== - secnd_frfi(nin)%P(3,jg) )
2081 * secnd_frfi(nin)%P(3,jg) = abs(fzt(i)- fnitzt(i))
2082#include "lockoff.inc"
2083 ENDIF
2084C------- total force
2085 fxi(i)=fxi(i) + fxt(i)
2086 fyi(i)=fyi(i) + fyt(i)
2087 fzi(i)=fzi(i) + fzt(i)
2088 efric_l(i) = dt1*(vx(i)*fxt(i)+vy(i)*fyt(i)+vz(i)*fzt(i))
2089 econvt = econvt + efric_l(i)
2090 ENDDO
2091C--------implicit non converge---
2092 ELSE
2093 DO i=1,jlt
2094 IF(pene(i) == zero)cycle
2095 fx = stif0(i)*vx(i)*dt12
2096 fy = stif0(i)*vy(i)*dt12
2097 fz = stif0(i)*vz(i)*dt12
2098 jg = nsvg(i)
2099 n = cand_n_n(i)
2100 IF(jg>0) THEN
2101c SECND_FR(4:6,N) is the old friction force
2102 fx = secnd_fr(4,n) + alpha*fx
2103 fy = secnd_fr(5,n) + alpha*fy
2104 fz = secnd_fr(6,n) + alpha*fz
2105 ftn = fx*n1(i) + fy*n2(i) + fz*n3(i)
2106 fx = fx - ftn*n1(i)
2107 fy = fy - ftn*n2(i)
2108 fz = fz - ftn*n3(i)
2109 fx = fx + fnitxt(i)
2110 fy = fy + fnityt(i)
2111 fz = fz + fnitzt(i)
2112 ft = fx*fx + fy*fy + fz*fz
2113 ft = max(ft,tiny)
2114 fn = fxi(i)**2+fyi(i)**2+fzi(i)**2
2115 beta = min(one,xmu(i)*sqrt(fn/ft))
2116 fxt(i) = fx * beta
2117 fyt(i) = fy * beta
2118 fzt(i) = fz * beta
2119 fxi(i)=fxi(i) + fxt(i)
2120 fyi(i)=fyi(i) + fyt(i)
2121 fzi(i)=fzi(i) + fzt(i)
2122 ELSE ! cas noeud remote en SPMD
2123 jg = -jg
2124 fx = secnd_frfi(nin)%P(4,jg) + alpha*fx
2125 fy = secnd_frfi(nin)%P(5,jg) + alpha*fy
2126 fz = secnd_frfi(nin)%P(6,jg) + alpha*fz
2127 ftn = fx*n1(i) + fy*n2(i) + fz*n3(i)
2128 fx = fx - ftn*n1(i)
2129 fy = fy - ftn*n2(i)
2130 fz = fz - ftn*n3(i)
2131 fx = fx + fnitxt(i)
2132 fy = fy + fnityt(i)
2133 fz = fz + fnitzt(i)
2134 ft = fx*fx + fy*fy + fz*fz
2135 ft = max(ft,tiny)
2136 fn = fxi(i)**2+fyi(i)**2+fzi(i)**2
2137 beta = min(one,xmu(i)*sqrt(fn/ft))
2138 fxt(i) = fx * beta
2139 fyt(i) = fy * beta
2140 fzt(i) = fz * beta
2141 fxi(i)=fxi(i) + fxt(i)
2142 fyi(i)=fyi(i) + fyt(i)
2143 fzi(i)=fzi(i) + fzt(i)
2144 ENDIF
2145C IFPEN(INDEX(I)) = 1
2146 ENDDO
2147 ENDIF
2148
2149 ELSE
2150C++ Orthotropic Friction
2151
2152 IF (inconv==1) THEN
2153#include "vectorize.inc"
2154 DO k=1,nfisot ! isotropic friction couples
2155 i = indexisot(k)
2156 IF(pene(i) == zero)cycle
2157 fx = stif0(i)*vx(i)*dt12
2158 fy = stif0(i)*vy(i)*dt12
2159 fz = stif0(i)*vz(i)*dt12
2160 jg = nsvg(i)
2161 IF(jg>0) THEN
2162 n = cand_n_n(i)
2163c SECND_FR(4:6,N) is the old friction force
2164 fx = secnd_fr(4,n) + alpha*fx
2165 fy = secnd_fr(5,n) + alpha*fy
2166 fz = secnd_fr(6,n) + alpha*fz
2167 ftn = fx*n1(i) + fy*n2(i) + fz*n3(i)
2168 fx = fx - ftn*n1(i)
2169 fy = fy - ftn*n2(i)
2170 fz = fz - ftn*n3(i)
2171 fx = fx + fnitxt(i)
2172 fy = fy + fnityt(i)
2173 fz = fz + fnitzt(i)
2174 ft = fx*fx + fy*fy + fz*fz
2175 ft = max(ft,tiny)
2176 fn = fxi(i)**2+fyi(i)**2+fzi(i)**2
2177 beta = min(one,xmu(i)*sqrt(fn/ft))
2178 fxt(i) = fx * beta
2179 fyt(i) = fy * beta
2180 fzt(i) = fz * beta
2181#include "lockon.inc"
2182 IF (abs(fxt(i)- fnitxt(i))>abs(secnd_fr(1,n))) secnd_fr(1,n) = fxt(i)- fnitxt(i)
2183 IF (abs(fyt(i)- fnityt(i))>abs(secnd_fr(2,n))) secnd_fr(2,n) = fyt(i)- fnityt(i)
2184 IF (abs(fzt(i)- fnitzt(i))>abs(secnd_fr(3,n))) secnd_fr(3,n) = fzt(i)- fnitzt(i)
2185C case equal abs but opposite sign
2186 IF ((fxt(i)- fnitxt(i))==-secnd_fr(1,n)) secnd_fr(1,n) = abs(fxt(i)- fnitxt(i))
2187 IF ((fyt(i)- fnityt(i))==-secnd_fr(2,n)) secnd_fr(2,n) = abs(fyt(i)- fnityt(i))
2188 IF ((fzt(i)- fnitzt(i))==-secnd_fr(3,n)) secnd_fr(3,n) = abs(fzt(i)- fnitzt(i))
2189#include "lockoff.inc"
2190 ELSE ! cas noeud remote en SPMD
2191 jg = -jg
2192c SECND_FRFI(NIN)%P(4:6,JG) is the old friction force
2193 fx = secnd_frfi(nin)%P(4,jg) + alpha*fx
2194 fy = secnd_frfi(nin)%P(5,jg) + alpha*fy
2195 fz = secnd_frfi(nin)%P(6,jg) + alpha*fz
2196 ftn = fx*n1(i) + fy*n2(i) + fz*n3(i)
2197 fx = fx - ftn*n1(i)
2198 fy = fy - ftn*n2(i)
2199 fz = fz - ftn*n3(i)
2200 fx = fx + fnitxt(i)
2201 fy = fy + fnityt(i)
2202 fz = fz + fnitzt(i)
2203 ft = fx*fx + fy*fy + fz*fz
2204 ft = max(ft,tiny)
2205 fn = fxi(i)**2+fyi(i)**2+fzi(i)**2
2206 beta = min(one,xmu(i)*sqrt(fn/ft))
2207 fxt(i) = fx * beta
2208 fyt(i) = fy * beta
2209 fzt(i) = fz * beta
2210#include "lockon.inc"
2211 IF ( abs(fxt(i)- fnitxt(i)) > abs(secnd_frfi(nin)%P(1,jg)) )
2212 * secnd_frfi(nin)%P(1,jg) = fxt(i)- fnitxt(i)
2213 IF ( abs(fyt(i)- fnityt(i)) > abs(secnd_frfi(nin)%P(2,jg)) )
2214 * secnd_frfi(nin)%P(2,jg) = fyt(i)- fnityt(i)
2215 IF ( abs(fzt(i)- fnitzt(i)) > abs(secnd_frfi(nin)%P(3,jg)) )
2216 * secnd_frfi(nin)%P(3,jg) = fzt(i)- fnitzt(i)
2217C
2218 IF ((fxt(i)- fnitxt(i))== - secnd_frfi(nin)%P(1,jg) )
2219 * secnd_frfi(nin)%P(1,jg) = abs(fxt(i)- fnitxt(i))
2220 IF ((fyt(i)- fnityt(i))== - secnd_frfi(nin)%P(2,jg) )
2221 * secnd_frfi(nin)%P(2,jg) = abs(fyt(i)- fnityt(i))
2222 IF ((fzt(i)- fnitzt(i))== - secnd_frfi(nin)%P(3,jg) )
2223 * secnd_frfi(nin)%P(3,jg) = abs(fzt(i)- fnitzt(i))
2224#include "lockoff.inc"
2225 ENDIF
2226C------- total force
2227 fxi(i)=fxi(i) + fxt(i)
2228 fyi(i)=fyi(i) + fyt(i)
2229 fzi(i)=fzi(i) + fzt(i)
2230 efric_l(i) = dt1*(vx(i)*fxt(i)+vy(i)*fyt(i)+vz(i)*fzt(i))
2231 econvt = econvt + efric_l(i)
2232 ENDDO
2233
2234#include "vectorize.inc"
2235 DO k=1,nforth ! Orthotropic friction couples
2236 i = indexorth(k)
2237 IF(pene(i) == zero)cycle
2238 fx = stif0(i)*vx(i)*dt12
2239 fy = stif0(i)*vy(i)*dt12
2240 fz = stif0(i)*vz(i)*dt12
2241 jg = nsvg(i)
2242 IF(jg>0) THEN
2243 n = cand_n_n(i)
2244c SECND_FR(4:6,N) is the old friction force
2245 fx = secnd_fr(4,n) + alpha*fx
2246 fy = secnd_fr(5,n) + alpha*fy
2247 fz = secnd_fr(6,n) + alpha*fz
2248 ftn = fx*n1(i) + fy*n2(i) + fz*n3(i)
2249 fx = fx - ftn*n1(i)
2250 fy = fy - ftn*n2(i)
2251 fz = fz - ftn*n3(i)
2252
2253 fx = fx + fnitxt(i)
2254 fy = fy + fnityt(i)
2255 fz = fz + fnitzt(i)
2256
2257 ftt1= fx*dir1(i,1) + fy*dir1(i,2) + fz*dir1(i,3)
2258 ftt2= fx*dir2(i,1) + fy*dir2(i,2) + fz*dir2(i,3)
2259 ft = ftt1*ftt1*an(k) + ftt2*ftt2*bn(k)
2260 ft = max(ft,em30)
2261 fn = fxi(i)**2+fyi(i)**2+fzi(i)**2
2262
2263 beta = fn/ft
2264
2265 IF(beta == 0 ) THEN
2266 fxt(i) = zero
2267 fyt(i) = zero
2268 fzt(i) = zero
2269 ELSEIF(beta > 1) THEN ! Inside the ellipse
2270 fxt(i) = fx
2271 fyt(i) = fy
2272 fzt(i) = fz
2273 ELSE ! outside the ellipse
2274
2275! Projection on local tangent of ellipse (outside of ellipse)
2276! ANS = (Fadh-Fproj).n
2277 nep1 =ftt1*an(k)/fn
2278 nep2 =ftt2*bn(k)/fn
2279 nep =nep1*nep1+nep2*nep2
2280 nep =sqrt(nep)
2281 ep=nep1*ftt1+nep2*ftt2
2282
2283 ans=(ep-sqrt(ep))/max(em20,nep)
2284 nep1 =nep1/max(em20,nep)
2285 nep2 =nep2/max(em20,nep)
2286! Projection on ellipse
2287 c11 =ftt1-ans*nep1
2288 c22 =ftt2-ans*nep2
2289 alphaf = atan(c22/c11)
2290 signc = ftt1/max(em20,abs(ftt1))
2291 csa = signc*abs(cos(alphaf))
2292 signc = ftt2/max(em20,abs(ftt2))
2293 sna = signc*abs(sin(alphaf))
2294! Ft computation
2295 ft = sqrt(fn / (csa*csa*an(k) + sna*sna*bn(k)))
2296 ftt1 = ft * csa
2297 ftt2 = ft * sna
2298
2299 fxt(i) = ftt1 * dir1(i,1) + ftt2 * dir2(i,1)
2300 fyt(i) = ftt1 * dir1(i,2) + ftt2 * dir2(i,2)
2301 fzt(i) = ftt1 * dir1(i,3) + ftt2 * dir2(i,3)
2302
2303 ENDIF
2304
2305#include "lockon.inc"
2306 IF (abs(fxt(i))>abs(secnd_fr(1,n))) secnd_fr(1,n) = fxt(i)
2307 IF (abs(fyt(i))>abs(secnd_fr(2,n))) secnd_fr(2,n) = fyt(i)
2308 IF (abs(fzt(i))>abs(secnd_fr(3,n))) secnd_fr(3,n) = fzt(i)
2309C case equal abs but opposite sign
2310 IF (fxt(i)==-secnd_fr(1,n)) secnd_fr(1,n) = abs(fxt(i))
2311 IF (fyt(i)==-secnd_fr(2,n)) secnd_fr(2,n) = abs(fyt(i))
2312 IF (fzt(i)==-secnd_fr(3,n)) secnd_fr(3,n) = abs(fzt(i))
2313#include "lockoff.inc"
2314 ELSE ! cas noeud remote en SPMD
2315 jg = -jg
2316c SECND_FRFI(NIN)%P(4:6,JG) is the old friction force
2317 fx = secnd_frfi(nin)%P(4,jg) + alpha*fx
2318 fy = secnd_frfi(nin)%P(5,jg) + alpha*fy
2319 fz = secnd_frfi(nin)%P(6,jg) + alpha*fz
2320 ftn = fx*n1(i) + fy*n2(i) + fz*n3(i)
2321 fx = fx - ftn*n1(i)
2322 fy = fy - ftn*n2(i)
2323 fz = fz - ftn*n3(i)
2324 fx = fx + fnitxt(i)
2325 fy = fy + fnityt(i)
2326 fz = fz + fnitzt(i)
2327 ftt1= fx*dir1(i,1) + fy*dir1(i,2) + fz*dir1(i,3)
2328 ftt2= fx*dir2(i,1) + fy*dir2(i,2) + fz*dir2(i,3)
2329 ft = ftt1*ftt1*an(k) + ftt2*ftt2*bn(k)
2330 ft = max(ft,em30)
2331 fn = fxi(i)**2+fyi(i)**2+fzi(i)**2
2332
2333 beta = fn/ft
2334
2335 IF(beta == 0 ) THEN
2336 fxt(i) = zero
2337 fyt(i) = zero
2338 fzt(i) = zero
2339 ELSEIF(beta > 1) THEN ! Inside the ellipse
2340 fxt(i) = fx
2341 fyt(i) = fy
2342 fzt(i) = fz
2343 ELSE ! outside the ellipse
2344
2345! Projection on local tangent of ellipse (outside of ellipse)
2346! ANS = (Fadh-Fproj).n
2347 nep1 =ftt1*an(k)/fn
2348 nep2 =ftt2*bn(k)/fn
2349 nep =nep1*nep1+nep2*nep2
2350 nep =sqrt(nep)
2351
2352 ep=nep1*ftt1+nep2*ftt2
2353
2354 ans=(ep-sqrt(ep))/max(em20,nep)
2355 nep1 =nep1/max(em20,nep)
2356 nep2 =nep2/max(em20,nep)
2357
2358! Projection on ellipse
2359 c11 =ftt1-ans*nep1
2360 c22 =ftt2-ans*nep2
2361
2362 alphaf = atan(c22/c11)
2363
2364 signc = ftt1/max(em20,abs(ftt1))
2365 csa = signc*abs(cos(alphaf))
2366 signc = ftt2/max(em20,abs(ftt2))
2367 sna = signc*abs(sin(alphaf))
2368! Ft computation
2369 ft = sqrt(fn / (csa*csa*an(k) + sna*sna*bn(k)))
2370 ftt1 = ft * csa
2371 ftt2 = ft * sna
2372
2373 fxt(i) = ftt1 * dir1(i,1) + ftt2 * dir2(i,1)
2374 fyt(i) = ftt1 * dir1(i,2) + ftt2 * dir2(i,2)
2375 fzt(i) = ftt1 * dir1(i,3) + ftt2 * dir2(i,3)
2376
2377 ENDIF
2378#include "lockon.inc"
2379 IF ( abs(fxt(i)) > abs(secnd_frfi(nin)%P(1,jg)) )
2380 * secnd_frfi(nin)%P(1,jg) = fxt(i)
2381 IF ( abs(fyt(i)) > abs(secnd_frfi(nin)%P(2,jg)) )
2382 * secnd_frfi(nin)%P(2,jg) = fyt(i)
2383 IF ( abs(fzt(i)) > abs(secnd_frfi(nin)%P(3,jg)) )
2384 * secnd_frfi(nin)%P(3,jg) = fzt(i)
2385C
2386 IF (fxt(i)== - secnd_frfi(nin)%P(1,jg) )
2387 * secnd_frfi(nin)%P(1,jg) = abs(fxt(i))
2388 IF (fyt(i)== - secnd_frfi(nin)%P(2,jg) )
2389 * secnd_frfi(nin)%P(2,jg) = abs(fyt(i))
2390 IF (fzt(i)== - secnd_frfi(nin)%P(3,jg) )
2391 * secnd_frfi(nin)%P(3,jg) = abs(fzt(i))
2392#include "lockoff.inc"
2393 ENDIF
2394C------- total force
2395 fxi(i)=fxi(i) + fxt(i)
2396 fyi(i)=fyi(i) + fyt(i)
2397 fzi(i)=fzi(i) + fzt(i)
2398 efric_l(i) = dt1*(vx(i)*fxt(i)+vy(i)*fyt(i)+vz(i)*fzt(i))
2399 econvt = econvt + efric_l(i)
2400 ENDDO
2401
2402
2403C--------implicit non converge---
2404 ELSE
2405#include "vectorize.inc"
2406 DO k=1,nfisot ! isotropic friction couples
2407 i = indexisot(k)
2408 IF(pene(i) == zero)cycle
2409 fx = stif0(i)*vx(i)*dt12
2410 fy = stif0(i)*vy(i)*dt12
2411 fz = stif0(i)*vz(i)*dt12
2412 jg = nsvg(i)
2413 n = cand_n_n(i)
2414 IF(jg>0) THEN
2415c SECND_FR(4:6,N) is the old friction force
2416 fx = secnd_fr(4,n) + alpha*fx
2417 fy = secnd_fr(5,n) + alpha*fy
2418 fz = secnd_fr(6,n) + alpha*fz
2419 ftn = fx*n1(i) + fy*n2(i) + fz*n3(i)
2420 fx = fx - ftn*n1(i)
2421 fy = fy - ftn*n2(i)
2422 fz = fz - ftn*n3(i)
2423 fx = fx + fnitxt(i)
2424 fy = fy + fnityt(i)
2425 fz = fz + fnitzt(i)
2426 ft = fx*fx + fy*fy + fz*fz
2427 ft = max(ft,tiny)
2428 fn = fxi(i)**2+fyi(i)**2+fzi(i)**2
2429 beta = min(one,xmu(i)*sqrt(fn/ft))
2430 fxt(i) = fx * beta
2431 fyt(i) = fy * beta
2432 fzt(i) = fz * beta
2433 fxi(i)=fxi(i) + fxt(i)
2434 fyi(i)=fyi(i) + fyt(i)
2435 fzi(i)=fzi(i) + fzt(i)
2436 ELSE ! cas noeud remote en SPMD
2437 jg = -jg
2438 fx = secnd_frfi(nin)%P(4,jg) + alpha*fx
2439 fy = secnd_frfi(nin)%P(5,jg) + alpha*fy
2440 fz = secnd_frfi(nin)%P(6,jg) + alpha*fz
2441 ftn = fx*n1(i) + fy*n2(i) + fz*n3(i)
2442 fx = fx - ftn*n1(i)
2443 fy = fy - ftn*n2(i)
2444 fz = fz - ftn*n3(i)
2445 fx = fx + fnitxt(i)
2446 fy = fy + fnityt(i)
2447 fz = fz + fnitzt(i)
2448 ft = fx*fx + fy*fy + fz*fz
2449 ft = max(ft,tiny)
2450 fn = fxi(i)**2+fyi(i)**2+fzi(i)**2
2451 beta = min(one,xmu(i)*sqrt(fn/ft))
2452 fxt(i) = fx * beta
2453 fyt(i) = fy * beta
2454 fzt(i) = fz * beta
2455 fxi(i)=fxi(i) + fxt(i)
2456 fyi(i)=fyi(i) + fyt(i)
2457 fzi(i)=fzi(i) + fzt(i)
2458 ENDIF
2459C IFPEN(INDEX(I)) = 1
2460 ENDDO
2461
2462#include "vectorize.inc"
2463 DO k=1,nforth ! Orthotropic friction couples
2464 i = indexorth(k)
2465 IF(pene(i) == zero)cycle
2466 fx = stif0(i)*vx(i)*dt12
2467 fy = stif0(i)*vy(i)*dt12
2468 fz = stif0(i)*vz(i)*dt12
2469 jg = nsvg(i)
2470 n = cand_n_n(i)
2471 IF(jg>0) THEN
2472c SECND_FR(4:6,N) is the old friction force
2473 fx = secnd_fr(4,n) + alpha*fx
2474 fy = secnd_fr(5,n) + alpha*fy
2475 fz = secnd_fr(6,n) + alpha*fz
2476 ftn = fx*n1(i) + fy*n2(i) + fz*n3(i)
2477 fx = fx - ftn*n1(i)
2478 fy = fy - ftn*n2(i)
2479 fz = fz - ftn*n3(i)
2480 fx = fx + fnitxt(i)
2481 fy = fy + fnityt(i)
2482 fz = fz + fnitzt(i)
2483
2484 ftt1= fx*dir1(i,1) + fy*dir1(i,2) + fz*dir1(i,3)
2485 ftt2= fx*dir2(i,1) + fy*dir2(i,2) + fz*dir2(i,3)
2486 ft = ftt1*ftt1*an(k) + ftt2*ftt2*bn(k)
2487 ft = max(ft,em30)
2488 fn = fxi(i)**2+fyi(i)**2+fzi(i)**2
2489
2490 beta = fn/ft
2491
2492 IF(beta == 0 ) THEN
2493 fxt(i) = zero
2494 fyt(i) = zero
2495 fzt(i) = zero
2496 ELSEIF(beta > 1) THEN ! Inside the ellipse
2497 fxt(i) = fx
2498 fyt(i) = fy
2499 fzt(i) = fz
2500 ELSE ! outside the ellipse
2501
2502! Projection on local tangent of ellipse (outside of ellipse)
2503! ANS = (Fadh-Fproj).n
2504 nep1 =ftt1*an(k)/fn
2505 nep2 =ftt2*bn(k)/fn
2506 nep =nep1*nep1+nep2*nep2
2507 nep =sqrt(nep)
2508
2509 ep=nep1*ftt1+nep2*ftt2
2510
2511 ans=(ep-sqrt(ep))/max(em20,nep)
2512 nep1 =nep1/max(em20,nep)
2513 nep2 =nep2/max(em20,nep)
2514
2515! Projection on ellipse
2516 c11 =ftt1-ans*nep1
2517 c22 =ftt2-ans*nep2
2518
2519 alphaf = atan(c22/c11)
2520
2521 signc = ftt1/max(em20,abs(ftt1))
2522 csa = signc*abs(cos(alphaf))
2523 signc = ftt2/max(em20,abs(ftt2))
2524 sna = signc*abs(sin(alphaf))
2525! Ft computation
2526 ft = sqrt(fn / (csa*csa*an(k) + sna*sna*bn(k)))
2527 ftt1 = ft * csa
2528 ftt2 = ft * sna
2529
2530 fxt(i) = ftt1 * dir1(i,1) + ftt2 * dir2(i,1)
2531 fyt(i) = ftt1 * dir1(i,2) + ftt2 * dir2(i,2)
2532 fzt(i) = ftt1 * dir1(i,3) + ftt2 * dir2(i,3)
2533
2534 ENDIF
2535 fxi(i)=fxi(i) + fxt(i)
2536 fyi(i)=fyi(i) + fyt(i)
2537 fzi(i)=fzi(i) + fzt(i)
2538 ELSE ! cas noeud remote en SPMD
2539 jg = -jg
2540 fx = secnd_frfi(nin)%P(4,jg) + alpha*fx
2541 fy = secnd_frfi(nin)%P(5,jg) + alpha*fy
2542 fz = secnd_frfi(nin)%P(6,jg) + alpha*fz
2543 ftn = fx*n1(i) + fy*n2(i) + fz*n3(i)
2544 fx = fx - ftn*n1(i)
2545 fy = fy - ftn*n2(i)
2546 fz = fz - ftn*n3(i)
2547 fx = fx + fnitxt(i)
2548 fy = fy + fnityt(i)
2549 fz = fz + fnitzt(i)
2550
2551 ftt1= fx*dir1(i,1) + fy*dir1(i,2) + fz*dir1(i,3)
2552 ftt2= fx*dir2(i,1) + fy*dir2(i,2) + fz*dir2(i,3)
2553 ft = ftt1*ftt1*an(k) + ftt2*ftt2*bn(k)
2554 ft = max(ft,em30)
2555 fn = fxi(i)**2+fyi(i)**2+fzi(i)**2
2556
2557 beta = fn/ft
2558
2559 IF(beta == 0 ) THEN
2560 fxt(i) = zero
2561 fyt(i) = zero
2562 fzt(i) = zero
2563 ELSEIF(beta > 1) THEN ! Inside the ellipse
2564 fxt(i) = fx
2565 fyt(i) = fy
2566 fzt(i) = fz
2567 ELSE ! outside the ellipse
2568
2569! Projection on local tangent of ellipse (outside of ellipse)
2570! ANS = (Fadh-Fproj).n
2571 nep1 =ftt1*an(k)/fn
2572 nep2 =ftt2*bn(k)/fn
2573 nep =nep1*nep1+nep2*nep2
2574 nep =sqrt(nep)
2575
2576 ep=nep1*ftt1+nep2*ftt2
2577
2578 ans=(ep-sqrt(ep))/max(em20,nep)
2579 nep1 =nep1/max(em20,nep)
2580 nep2 =nep2/max(em20,nep)
2581
2582! Projection on ellipse
2583 c11 =ftt1-ans*nep1
2584 c22 =ftt2-ans*nep2
2585
2586 alphaf = atan(c22/c11)
2587
2588 signc = ftt1/max(em20,abs(ftt1))
2589 csa = signc*abs(cos(alphaf))
2590 signc = ftt2/max(em20,abs(ftt2))
2591 sna = signc*abs(sin(alphaf))
2592! Ft computation
2593 ft = sqrt(fn / (csa*csa*an(k) + sna*sna*bn(k)))
2594 ftt1 = ft * csa
2595 ftt2 = ft * sna
2596
2597 fxt(i) = ftt1 * dir1(i,1) + ftt2 * dir2(i,1)
2598 fyt(i) = ftt1 * dir1(i,2) + ftt2 * dir2(i,2)
2599 fzt(i) = ftt1 * dir1(i,3) + ftt2 * dir2(i,3)
2600
2601 ENDIF
2602
2603 fxi(i)=fxi(i) + fxt(i)
2604 fyi(i)=fyi(i) + fyt(i)
2605 fzi(i)=fzi(i) + fzt(i)
2606 ENDIF
2607C IFPEN(INDEX(I)) = 1
2608 ENDDO
2609
2610 ENDIF !imp
2611
2612 ENDIF ! iorth
2613
2614
2615
2616 ENDIF !ifq
2617C
2618C---------------------------------
2619 IF((anim_v(12)+outp_v(12)+h3d_data%N_VECT_PCONT>0.AND.
2620 . ((tt>=tanim.AND.tt<=tanim_stop).OR.tt>=toutp.OR.(tt>=h3d_data%TH3D.AND.tt<=h3d_data%TH3D_STOP).OR.
2621 . (manim>=4.AND.manim<=15).OR.h3d_data%MH3D/=0))
2622 . .OR.h3d_data%N_VECT_PCONT_MAX>0)THEN
2623 IF (inconv==1) THEN
2624#include "lockon.inc"
2625 DO i=1,jlt
2626 IF(pene(i) == zero)cycle
2627 ftcont(1,ix1(i)) =ftcont(1,ix1(i)) + fxt(i)*h1(i)
2628 ftcont(2,ix1(i)) =ftcont(2,ix1(i)) + fyt(i)*h1(i)
2629 ftcont(3,ix1(i)) =ftcont(3,ix1(i)) + fzt(i)*h1(i)
2630 ftcont(1,ix2(i)) =ftcont(1,ix2(i)) + fxt(i)*h2(i)
2631 ftcont(2,ix2(i)) =ftcont(2,ix2(i)) + fyt(i)*h2(i)
2632 ftcont(3,ix2(i)) =ftcont(3,ix2(i)) + fzt(i)*h2(i)
2633 ftcont(1,ix3(i)) =ftcont(1,ix3(i)) + fxt(i)*h3(i)
2634 ftcont(2,ix3(i)) =ftcont(2,ix3(i)) + fyt(i)*h3(i)
2635 ftcont(3,ix3(i)) =ftcont(3,ix3(i)) + fzt(i)*h3(i)
2636 ftcont(1,ix4(i)) =ftcont(1,ix4(i)) + fxt(i)*h4(i)
2637 ftcont(2,ix4(i)) =ftcont(2,ix4(i)) + fyt(i)*h4(i)
2638 ftcont(3,ix4(i)) =ftcont(3,ix4(i)) + fzt(i)*h4(i)
2639 jg = nsvg(i)
2640 IF(jg>0) THEN
2641C en SPMD : traitement a refaire apres reception noeud remote si JG < 0
2642 IF (jg >numnod) THEN
2643 ig = jg - numnod
2644 CALL i24for1_fic(3 ,irtse ,nsne ,is2se ,is2pt ,
2645 + ig ,fxt(i) ,fyt(i) ,fzt(i) ,ftcont ,
2646 + -1 )
2647 ELSE
2648 ftcont(1,jg)=ftcont(1,jg)- fxt(i)
2649 ftcont(2,jg)=ftcont(2,jg)- fyt(i)
2650 ftcont(3,jg)=ftcont(3,jg)- fzt(i)
2651 ENDIF
2652 ELSE
2653 jg=-jg
2654 IF(isedge_fi(nin)%P(jg)==1)THEN
2655 CALL i24for1_fic(3 ,irtse_fi(nin)%P(1,1) ,nsne ,is2se_fi(nin)%P(1,1) ,is2pt_fi(nin)%P(1) ,
2656 + jg ,fxt(i) ,fyt(i) ,fzt(i) ,ftconti(nin)%P(1,1) ,
2657 + -1 )
2658 ELSE ! cas noeud remote en SPMD non edge
2659 ftconti(nin)%P(1,jg)=ftconti(nin)%P(1,jg)-fxt(i)
2660 ftconti(nin)%P(2,jg)=ftconti(nin)%P(2,jg)-fyt(i)
2661 ftconti(nin)%P(3,jg)=ftconti(nin)%P(3,jg)-fzt(i)
2662 ENDIF
2663 END IF
2664 ENDDO
2665#include "lockoff.inc"
2666 END IF !(INCONV==1) THEN
2667 ENDIF
2668C
2669C---------------------------------
2670 fsav22= zero
2671 fsav23= zero
2672 fsav24= zero
2673 fsav29= zero
2674 DO i=1,jlt
2675 IF(pene(i) == zero)cycle
2676 impx=fxt(i)*dt12
2677 impy=fyt(i)*dt12
2678 impz=fzt(i)*dt12
2679 fsav4 =fsav4 +impx
2680 fsav5 =fsav5 +impy
2681 fsav6 =fsav6 +impz
2682 impx=fxi(i)*dt12
2683 impy=fyi(i)*dt12
2684 impz=fzi(i)*dt12
2685 fsav12=fsav12+abs(impx)
2686 fsav13=fsav13+abs(impy)
2687 fsav14=fsav14+abs(impz)
2688 fsav15=fsav15+sqrt(impx*impx+impy*impy+impz*impz)
2689 xp(i)=xi(i)+pene(i)*n1(i)
2690 yp(i)=yi(i)+pene(i)*n2(i)
2691 zp(i)=zi(i)+pene(i)*n3(i)
2692 fsav22=fsav22+yp(i)*impz-zp(i)*impy
2693 fsav23=fsav23+zp(i)*impx-xp(i)*impz
2694 fsav24=fsav24+xp(i)*impy-yp(i)*impx
2695 ENDDO
2696
2697 IF(intcarea > 0) THEN
2698 DO i=1,jlt
2699 IF(pene(i) == zero)cycle
2700 jg = nsvg(i)
2701 IF(jg>0) THEN ! Only local nodes in i25for3, remote nodes will be taken into account later
2702 IF(jg <= numnod) THEN
2703 fsav29=fsav29+parameters%INTAREAN(jg)
2704 ELSE
2705 ig = jg - numnod
2706 CALL i24intarea_fic(irtse ,nsne ,is2se ,is2pt ,ig ,
2707 + nrtse , numnod , parameters%INTAREAN, arean_fic)
2708 fsav29=fsav29 + arean_fic
2709 ENDIF
2710 ELSE
2711 jg=-jg
2712 IF(isedge_fi(nin)%P(jg)==1)THEN
2713 CALL i24intarea_fic(irtse_fi(nin)%P ,nsnr ,is2se_fi(nin)%P ,is2pt_fi(nin)%P ,jg ,
2714 + nsnr , nsnr ,intareanfi(nin)%P, arean_fic)
2715 fsav29=fsav29 + arean_fic
2716 ELSE ! cas noeud remote en SPMD non edge
2717 fsav29=fsav29 + intareanfi(nin)%P(jg)
2718 ENDIF
2719 ENDIF
2720 ENDDO
2721 ENDIF
2722#include "lockon.inc"
2723 fsav(4) = fsav(4) + fsav4
2724 fsav(5) = fsav(5) + fsav5
2725 fsav(6) = fsav(6) + fsav6
2726 fsav(12) = fsav(12) + fsav12
2727 fsav(13) = fsav(13) + fsav13
2728 fsav(14) = fsav(14) + fsav14
2729 fsav(15) = fsav(15) + fsav15
2730 fsav(22) = fsav(22) + fsav22
2731 fsav(23) = fsav(23) + fsav23
2732 fsav(24) = fsav(24) + fsav24
2733 fsav(26) = fsav(26) + econtt
2734 fsav(27) = fsav(27) + econvt
2735 fsav(28) = fsav(28) + econtdt
2736 fsav(29) = fsav(29) + fsav29
2737#include "lockoff.inc"
2738C
2739 IF(isensint(1)/=0) THEN
2740 DO i=1,jlt
2741 fsavparit(1,4,i+nft) = fxi(i)
2742 fsavparit(1,5,i+nft) = fyi(i)
2743 fsavparit(1,6,i+nft) = fzi(i)
2744 ENDDO
2745 ENDIF
2746C---------------------------------
2747C SORTIES TH PAR SOUS INTERFACE
2748C---------------------------------
2749 IF(nisub/=0)THEN
2750 DO i=1,jlt
2751 IF(pene(i) == zero)cycle
2752 nn = nsvg(i)
2753 IF(nn>0)THEN
2754 in=cn_loc(i)
2755 IF (msegtyp(ce_loc(i)) < 0) THEN
2756 ie= - msegtyp(ce_loc(i))
2757 ELSE
2758 ie = ce_loc(i)
2759 ENDIF
2760 jj =addsubs(in)
2761 kk =addsubm(ie)
2762 DO WHILE(jj<addsubs(in+1))
2763 jsub=lisubs(jj)
2764 itypsub = typsub(jsub)
2765
2766 IF(itypsub == 1 ) THEN
2767 iss1 = bitget(inflg_subs(jj),0)
2768 iss2 = bitget(inflg_subs(jj),1)
2769 igrn = bitget(inflg_subs(jj),2)
2770 ksub=lisubm(kk)
2771 DO WHILE((ksub<=jsub).AND.(kk<addsubm(ie+1)))
2772 ims1 = bitget(inflg_subm(kk),0)
2773 ims2 = bitget(inflg_subm(kk),1)
2774 IF(ksub==jsub)THEN
2775! S and M candidates on the same sub_interface
2776 IF(.NOT.(((ims1 == 1 .OR. ims2 == 1) .AND. igrn == 1).OR.
2777 . (ims1 == 1 .AND. iss2 == 1).OR.
2778 . (ims2 == 1 .AND. iss1 == 1))) THEN
2779 kk=kk+1
2780 ksub=lisubm(kk)
2781 cycle
2782 END IF
2783 impx=fxt(i)*dt12
2784 impy=fyt(i)*dt12
2785 impz=fzt(i)*dt12
2786C main side :
2787 fsavsub1(4,jsub)=fsavsub1(4,jsub)+impx
2788 fsavsub1(5,jsub)=fsavsub1(5,jsub)+impy
2789 fsavsub1(6,jsub)=fsavsub1(6,jsub)+impz
2790C
2791 impx=fxi(i)*dt12
2792 impy=fyi(i)*dt12
2793 impz=fzi(i)*dt12
2794 fsavsub1(12,jsub)=fsavsub1(12,jsub)+abs(impx)
2795 fsavsub1(13,jsub)=fsavsub1(13,jsub)+abs(impy)
2796 fsavsub1(14,jsub)=fsavsub1(14,jsub)+abs(impz)
2797C
2798 IF(isensint(jsub+1)/=0) THEN
2799 fsavparit(jsub+1,4,i+nft) = fxt(i)
2800 fsavparit(jsub+1,5,i+nft) = fyt(i)
2801 fsavparit(jsub+1,6,i+nft) = fzt(i)
2802 ENDIF
2803C
2804 fsavsub1(15,jsub)= fsavsub1(15,jsub)
2805 . +sqrt(impx*impx+impy*impy+impz*impz)
2806 fsavsub1(22,jsub)=fsavsub1(22,jsub)
2807 . +yp(i)*impz-zp(i)*impy
2808 fsavsub1(23,jsub)=fsavsub1(23,jsub)
2809 . +zp(i)*impx-xp(i)*impz
2810 fsavsub1(24,jsub)=fsavsub1(24,jsub)
2811 . +xp(i)*impy-yp(i)*impx
2812C
2813 END IF
2814
2815 kk=kk+1
2816 ksub=lisubm(kk)
2817 ENDDO
2818 jj=jj+1
2819
2820 ELSEIF(itypsub == 2 ) THEN ! Inter =0 : collecting forces from all inter with only 1 surface : second side
2821
2822 impx=fxt(i)*dt12
2823 impy=fyt(i)*dt12
2824 impz=fzt(i)*dt12
2825C main side :
2826 fsavsub1(4,jsub)=fsavsub1(4,jsub)+impx
2827 fsavsub1(5,jsub)=fsavsub1(5,jsub)+impy
2828 fsavsub1(6,jsub)=fsavsub1(6,jsub)+impz
2829C
2830 impx=fxi(i)*dt12
2831 impy=fyi(i)*dt12
2832 impz=fzi(i)*dt12
2833 fsavsub1(12,jsub)=fsavsub1(12,jsub)+abs(impx)
2834 fsavsub1(13,jsub)=fsavsub1(13,jsub)+abs(impy)
2835 fsavsub1(14,jsub)=fsavsub1(14,jsub)+abs(impz)
2836C
2837 IF(isensint(jsub+1)/=0) THEN
2838 fsavparit(jsub+1,4,i+nft) = fxt(i)
2839 fsavparit(jsub+1,5,i+nft) = fyt(i)
2840 fsavparit(jsub+1,6,i+nft) = fzt(i)
2841 ENDIF
2842C
2843 fsavsub1(15,jsub)= fsavsub1(15,jsub)
2844 . +sqrt(impx*impx+impy*impy+impz*impz)
2845 fsavsub1(22,jsub)=fsavsub1(22,jsub)
2846 . +yp(i)*impz-zp(i)*impy
2847 fsavsub1(23,jsub)=fsavsub1(23,jsub)
2848 . +zp(i)*impx-xp(i)*impz
2849 fsavsub1(24,jsub)=fsavsub1(24,jsub)
2850 . +xp(i)*impy-yp(i)*impx
2851
2852 jj=jj+1
2853
2854
2855 ELSEIF(itypsub == 3 ) THEN ! Inter =0 : collecting forces from all inter with 2 surfaces
2856C
2857C Find if node is on Surface S1 -> S2 or S2 ->S1
2858 iss2 = bitget(inflg_subs(jj),0)
2859 iss1 = bitget(inflg_subs(jj),1)
2860 ksub=lisubm(kk)
2861 DO WHILE((ksub<=jsub).AND.(kk<addsubm(ie+1)))
2862 ims2 = bitget(inflg_subm(kk),0)
2863 ims1 = bitget(inflg_subm(kk),1)
2864 IF(ksub==jsub)THEN
2865! S and M candidates on the same sub_interface
2866 IF(.NOT.((ims1 == 1 .AND. iss2 == 1).OR.
2867 . (ims2 == 1 .AND. iss1 == 1))) THEN
2868 kk=kk+1
2869 ksub=lisubm(kk)
2870 cycle
2871 END IF
2872
2873 impx=fxt(i)*dt12
2874 impy=fyt(i)*dt12
2875 impz=fzt(i)*dt12
2876 IF(ims2 > 0) THEN
2877C main side :
2878 fsavsub1(4,jsub)=fsavsub1(4,jsub)-impx
2879 fsavsub1(5,jsub)=fsavsub1(5,jsub)-impy
2880 fsavsub1(6,jsub)=fsavsub1(6,jsub)-impz
2881 ELSE
2882 fsavsub1(4,jsub)=fsavsub1(4,jsub)+impx
2883 fsavsub1(5,jsub)=fsavsub1(5,jsub)+impy
2884 fsavsub1(6,jsub)=fsavsub1(6,jsub)+impz
2885 ENDIF
2886C
2887 impx=fxi(i)*dt12
2888 impy=fyi(i)*dt12
2889 impz=fzi(i)*dt12
2890 fsavsub1(12,jsub)=fsavsub1(12,jsub)+abs(impx)
2891 fsavsub1(13,jsub)=fsavsub1(13,jsub)+abs(impy)
2892 fsavsub1(14,jsub)=fsavsub1(14,jsub)+abs(impz)
2893C
2894 IF(isensint(jsub+1)/=0) THEN
2895 IF(ims2 > 0) THEN
2896 fsavparit(jsub+1,4,i+nft) = -fxt(i)
2897 fsavparit(jsub+1,5,i+nft) = -fyt(i)
2898 fsavparit(jsub+1,6,i+nft) = -fzt(i)
2899 ELSE
2900 fsavparit(jsub+1,4,i+nft) = fxt(i)
2901 fsavparit(jsub+1,5,i+nft) = fyt(i)
2902 fsavparit(jsub+1,6,i+nft) = fzt(i)
2903 ENDIF
2904 ENDIF
2905C
2906 fsavsub1(15,jsub)= fsavsub1(15,jsub)
2907 . +sqrt(impx*impx+impy*impy+impz*impz)
2908 fsavsub1(22,jsub)=fsavsub1(22,jsub)
2909 . +yp(i)*impz-zp(i)*impy
2910 fsavsub1(23,jsub)=fsavsub1(23,jsub)
2911 . +zp(i)*impx-xp(i)*impz
2912 fsavsub1(24,jsub)=fsavsub1(24,jsub)
2913 . +xp(i)*impy-yp(i)*impx
2914C
2915 ENDIF
2916 kk=kk+1
2917 ksub=lisubm(kk)
2918 ENDDO
2919 jj=jj+1
2920
2921
2922 ENDIF
2923
2924 END DO
2925 END IF
2926
2927 IF (msegtyp(ce_loc(i)) < 0) THEN
2928 ie= - msegtyp(ce_loc(i))
2929 ELSE
2930 ie = ce_loc(i)
2931 ENDIF
2932 IF(ie > nrtm) ie=ie-nrtm
2933
2934 kk =addsubm(ie) ! Addresse du M
2935 DO WHILE(kk<addsubm(ie+1))
2936! all sub interfaces of S
2937 ksub=lisubm(kk)
2938! KSUB Croissant
2939 itypsub = typsub(ksub)
2940 IF(itypsub == 2 ) THEN ! Inter =0 : collecting forces from all inter with only 1 surface
2941 impx=-fxt(i)*dt12
2942 impy=-fyt(i)*dt12
2943 impz=-fzt(i)*dt12
2944C main side :
2945 fsavsub1(4,ksub)=fsavsub1(4,ksub)+impx
2946 fsavsub1(5,ksub)=fsavsub1(5,ksub)+impy
2947 fsavsub1(6,ksub)=fsavsub1(6,ksub)+impz
2948C
2949 impx=fxi(i)*dt12
2950 impy=fyi(i)*dt12
2951 impz=fzi(i)*dt12
2952 fsavsub1(12,ksub)=fsavsub1(12,ksub)+abs(impx)
2953 fsavsub1(13,ksub)=fsavsub1(13,ksub)+abs(impy)
2954 fsavsub1(14,ksub)=fsavsub1(14,ksub)+abs(impz)
2955C
2956 IF(isensint(ksub+1)/=0) THEN
2957 fsavparit(ksub+1,4,i+nft) = -fxt(i)
2958 fsavparit(ksub+1,5,i+nft) = -fyt(i)
2959 fsavparit(ksub+1,6,i+nft) = -fzt(i)
2960 ENDIF
2961C
2962 fsavsub1(15,ksub)= fsavsub1(15,ksub)
2963 . +sqrt(impx*impx+impy*impy+impz*impz)
2964 fsavsub1(22,ksub)=fsavsub1(22,ksub)
2965 . +yp(i)*impz-zp(i)*impy
2966 fsavsub1(23,ksub)=fsavsub1(23,ksub)
2967 . +zp(i)*impx-xp(i)*impz
2968 fsavsub1(24,ksub)=fsavsub1(24,ksub)
2969 . +xp(i)*impy-yp(i)*impx
2970
2971 ENDIF
2972 kk=kk+1
2973 ENDDO
2974
2975 END DO
2976
2977
2978 IF(nspmd>1) THEN
2979C loop split because of a PGI bug
2980 DO i=1,jlt
2981 IF(pene(i) == zero)cycle
2982 nn = nsvg(i)
2983 IF(nn<0)THEN
2984 nn = -nn
2985 IF (msegtyp(ce_loc(i)) < 0) THEN
2986 ie= - msegtyp(ce_loc(i))
2987 ELSE
2988 ie = ce_loc(i)
2989 ENDIF
2990 jj =addsubsfi(nin)%P(nn)
2991 kk =addsubm(ie)
2992 DO WHILE(jj<addsubsfi(nin)%P(nn+1))
2993 jsub=lisubsfi(nin)%P(jj)
2994 itypsub = typsub(jsub)
2995
2996 IF(itypsub == 1 ) THEN
2997 iss1 = bitget(inflg_subsfi(nin)%P(jj),0)
2998 iss2 = bitget(inflg_subsfi(nin)%P(jj),1)
2999 igrn = bitget(inflg_subsfi(nin)%P(jj),2)
3000 ksub=lisubm(kk)
3001 DO WHILE((ksub<=jsub).AND.(kk<addsubm(ie+1)))
3002 ims1 = bitget(inflg_subm(kk),0)
3003 ims2 = bitget(inflg_subm(kk),1)
3004 IF(ksub==jsub)THEN
3005 IF(.NOT.(((ims1 == 1 .OR. ims2 == 1) .AND. igrn == 1).OR.
3006 . (ims1 == 1 .AND. iss2 == 1).OR.
3007 . (ims2 == 1 .AND. iss1 == 1))) THEN
3008 kk=kk+1
3009 ksub=lisubm(kk)
3010 cycle
3011 END IF
3012 impx=fxt(i)*dt12
3013 impy=fyt(i)*dt12
3014 impz=fzt(i)*dt12
3015C main side :
3016 fsavsub1(4,jsub)=fsavsub1(4,jsub)+impx
3017 fsavsub1(5,jsub)=fsavsub1(5,jsub)+impy
3018 fsavsub1(6,jsub)=fsavsub1(6,jsub)+impz
3019C
3020 impx=fxi(i)*dt12
3021 impy=fyi(i)*dt12
3022 impz=fzi(i)*dt12
3023 fsavsub1(12,jsub)=fsavsub1(12,jsub)+abs(impx)
3024 fsavsub1(13,jsub)=fsavsub1(13,jsub)+abs(impy)
3025 fsavsub1(14,jsub)=fsavsub1(14,jsub)+abs(impz)
3026C
3027 IF(isensint(jsub+1)/=0) THEN
3028 fsavparit(jsub+1,4,i+nft) = fxt(i)
3029 fsavparit(jsub+1,5,i+nft) = fyt(i)
3030 fsavparit(jsub+1,6,i+nft) = fzt(i)
3031 ENDIF
3032C
3033 fsavsub1(15,jsub)= fsavsub1(15,jsub)
3034 . +sqrt(impx*impx+impy*impy+impz*impz)
3035 fsavsub1(22,jsub)=fsavsub1(22,jsub)
3036 . +yp(i)*impz-zp(i)*impy
3037 fsavsub1(23,jsub)=fsavsub1(23,jsub)
3038 . +zp(i)*impx-xp(i)*impz
3039 fsavsub1(24,jsub)=fsavsub1(24,jsub)
3040 . +xp(i)*impy-yp(i)*impx
3041
3042C
3043 END IF
3044
3045 kk=kk+1
3046 ksub=lisubm(kk)
3047 ENDDO
3048 jj=jj+1
3049
3050 ELSEIF(itypsub == 2 ) THEN ! Inter =0 : collecting forces from all inter with only 1 surface
3051
3052 impx=fxt(i)*dt12
3053 impy=fyt(i)*dt12
3054 impz=fzt(i)*dt12
3055C main side :
3056 fsavsub1(4,jsub)=fsavsub1(4,jsub)+impx
3057 fsavsub1(5,jsub)=fsavsub1(5,jsub)+impy
3058 fsavsub1(6,jsub)=fsavsub1(6,jsub)+impz
3059C
3060 impx=fxi(i)*dt12
3061 impy=fyi(i)*dt12
3062 impz=fzi(i)*dt12
3063 fsavsub1(12,jsub)=fsavsub1(12,jsub)+abs(impx)
3064
3065 fsavsub1(13,jsub)=fsavsub1(13,jsub)+abs(impy)
3066 fsavsub1(14,jsub)=fsavsub1(14,jsub)+abs(impz)
3067C
3068 IF(isensint(jsub+1)/=0) THEN
3069 fsavparit(jsub+1,4,i+nft) = fxt(i)
3070 fsavparit(jsub+1,5,i+nft) = fyt(i)
3071 fsavparit(jsub+1,6,i+nft) = fzt(i)
3072 ENDIF
3073C
3074 fsavsub1(15,jsub)= fsavsub1(15,jsub)
3075 . +sqrt(impx*impx+impy*impy+impz*impz)
3076 fsavsub1(22,jsub)=fsavsub1(22,jsub)
3077 . +yp(i)*impz-zp(i)*impy
3078 fsavsub1(23,jsub)=fsavsub1(23,jsub)
3079 . +zp(i)*impx-xp(i)*impz
3080 fsavsub1(24,jsub)=fsavsub1(24,jsub)
3081 . +xp(i)*impy-yp(i)*impx
3082
3083 jj=jj+1
3084
3085 ELSEIF(itypsub == 3 ) THEN ! Inter =0 : collecting forces from all inter with 2 surfaces
3086
3087 iss2 = bitget(inflg_subsfi(nin)%P(jj),0)
3088 iss1 = bitget(inflg_subsfi(nin)%P(jj),1)
3089 ksub=lisubm(kk)
3090 DO WHILE((ksub<=jsub).AND.(kk<addsubm(ie+1)))
3091 ims2 = bitget(inflg_subm(kk),0)
3092 ims1 = bitget(inflg_subm(kk),1)
3093 IF(ksub==jsub)THEN
3094 IF(.NOT.((ims1 == 1 .AND. iss2 == 1).OR.
3095 . (ims2 == 1 .AND. iss1 == 1))) THEN
3096 kk=kk+1
3097 ksub=lisubm(kk)
3098 cycle
3099 END IF
3100
3101 impx=fxi(i)*dt12
3102 impy=fyi(i)*dt12
3103 impz=fzi(i)*dt12
3104
3105 IF(ims2 > 0)THEN
3106 fsavsub1(1,jsub)=fsavsub1(1,jsub)-impx
3107 fsavsub1(2,jsub)=fsavsub1(2,jsub)-impy
3108 fsavsub1(3,jsub)=fsavsub1(3,jsub)-impz
3109 fsavsub1(11,jsub)=fsavsub1(11,jsub)-fni(i)*dt12
3110 ELSE
3111 fsavsub1(1,jsub)=fsavsub1(1,jsub)+impx
3112 fsavsub1(2,jsub)=fsavsub1(2,jsub)+impy
3113 fsavsub1(3,jsub)=fsavsub1(3,jsub)+impz
3114 fsavsub1(11,jsub)=fsavsub1(11,jsub)+fni(i)*dt12
3115 END IF
3116C
3117 IF(isensint(jsub+1)/=0) THEN
3118 IF(ims2 > 0 ) THEN
3119 fsavparit(jsub+1,4,i+nft) = -fxt(i)
3120 fsavparit(jsub+1,5,i+nft) = -fyt(i)
3121 fsavparit(jsub+1,6,i+nft) = -fzt(i)
3122 ELSE
3123 fsavparit(jsub+1,4,i+nft) = fxt(i)
3124 fsavparit(jsub+1,5,i+nft) = fyt(i)
3125 fsavparit(jsub+1,6,i+nft) = fzt(i)
3126 ENDIF
3127 ENDIF
3128
3129 fsavsub1(8,jsub) =fsavsub1(8,jsub) +abs(impx)
3130 fsavsub1(9,jsub) =fsavsub1(9,jsub) +abs(impy)
3131 fsavsub1(10,jsub)=fsavsub1(10,jsub)+abs(impz)
3132
3133 ENDIF
3134 kk=kk+1
3135 ksub=lisubm(kk)
3136 ENDDO
3137 jj=jj+1
3138
3139 ENDIF
3140 END DO
3141 END IF
3142 END DO
3143 END IF
3144#include "lockon.inc"
3145 DO jsub=1,nisub
3146 nsub=lisub(jsub)
3147 DO j=1,15
3148 fsavsub(j,nsub)=fsavsub(j,nsub)+fsavsub1(j,jsub)
3149 END DO
3150 fsavsub(22,nsub)=fsavsub(22,nsub)+fsavsub1(22,jsub)
3151 fsavsub(23,nsub)=fsavsub(23,nsub)+fsavsub1(23,jsub)
3152 fsavsub(24,nsub)=fsavsub(24,nsub)+fsavsub1(24,jsub)
3153 IF(parameters%INTCAREA > 0) fsavsub(29,nsub)=fsavsub(29,nsub)+fsavsub1(25,jsub)
3154 END DO
3155#include "lockoff.inc"
3156 END IF
3157C---------------------------------
3158#include "lockon.inc"
3159 IF (inconv==1) THEN
3160 econtv = econtv + econvt ! Frictional Energy
3161 econt = econt + econtt ! Elastic Energy
3162 econtd = econtd + econtdt ! Damping Energy
3163 END IF !(INCONV==1) THEN
3164#include "lockoff.inc"
3165C---------------------------------
3166C---------------------------------
3167 IF(interefric >0)THEN
3168 IF (inconv==1) THEN
3169#include "lockon.inc"
3170 DO i=1,jlt
3171 IF(pene(i) == zero)cycle
3172 efric_lm = half*efric_l(i)
3173 efric(interefric,ix1(i)) = efric(interefric,ix1(i)) + h1(i)*efric_lm
3174 efric(interefric,ix2(i)) = efric(interefric,ix2(i)) + h2(i)*efric_lm
3175 efric(interefric,ix3(i)) = efric(interefric,ix3(i)) + h3(i)*efric_lm
3176 efric(interefric,ix4(i)) = efric(interefric,ix4(i)) + h4(i)*efric_lm
3177 jg = nsvg(i)
3178 efric_ls = half*efric_l(i)
3179 IF(jg>0) THEN
3180 efric(interefric,jg) = efric(interefric,jg) + efric_ls
3181 ELSE ! node remote
3182 jg = -jg
3183 efricfi(nin)%P(jg)=efricfi(nin)%P(jg)+ efric_ls
3184 ENDIF
3185 ENDDO
3186#include "lockoff.inc"
3187 END IF !(INCONV==1) THEN
3188 ENDIF
3189C---------------------------------
3190 IF(h3d_data%N_SCAL_CSE_FRIC >0)THEN
3191 IF (inconv==1) THEN
3192#include "lockon.inc"
3193 DO i=1,jlt
3194 IF(pene(i) == zero)cycle
3195 efric_lm = half*efric_l(i)
3196 efricg(ix1(i)) = efricg(ix1(i)) + h1(i)*efric_lm
3197 efricg(ix2(i)) = efricg(ix2(i)) + h2(i)*efric_lm
3198 efricg(ix3(i)) = efricg(ix3(i)) + h3(i)*efric_lm
3199 efricg(ix4(i)) = efricg(ix4(i)) + h4(i)*efric_lm
3200 jg = nsvg(i)
3201 efric_ls = half*efric_l(i)
3202 IF(jg>0) THEN
3203 efricg(jg) = efricg(jg) + efric_ls
3204 ELSE ! node remote
3205 jg = -jg
3206 efricgfi(nin)%P(jg)=efricgfi(nin)%P(jg)+ efric_ls
3207 ENDIF
3208 ENDDO
3209#include "lockoff.inc"
3210 END IF !(INCONV==1) THEN
3211 ENDIF
3212C---------------------------------
3213 IF(kdtint==1)THEN
3214 IF( (visc/=zero)
3215 . .AND.(ivis2==0.OR.ivis2==1))THEN
3216 DO i=1,jlt
3217 IF(pene(i) == zero)cycle
3218C C(I)=2.*C(I)
3219C
3220 IF(msi(i)==zero)THEN
3221 ks(i) =zero
3222 cs(i) =zero
3223 stv(i)=zero
3224 ELSE
3225 cx = four*c(i)*c(i)
3226 cy = eight*msi(i)*kt(i)
3227 aux = sqrt(cx+cy)+two*c(i)
3228 stv(i)= kt(i)*aux*aux/max(cy,em30)
3229 aux = two*cf(i)*cf(i)/max(msi(i),em20)
3230 IF(aux>stv(i))THEN
3231 ks(i) =zero
3232 cs(i) =cf(i)
3233 stv(i)=aux
3234 ELSE
3235 ks(i)= kt(i)
3236 cs(i) =c(i)
3237 ENDIF
3238 ENDIF
3239C
3240 j1=ix1(i)
3241 IF(ms(j1)==zero)THEN
3242 k1(i) =zero
3243 c1(i) =zero
3244 st1(i)=zero
3245 ELSE
3246 k1(i)=kt(i)*abs(h1(i))
3247 c1(i)=c(i)*abs(h1(i))
3248 cx =four*c1(i)*c1(i)
3249 cy =eight*ms(j1)*k1(i)
3250 aux = sqrt(cx+cy)+two*c1(i)
3251 st1(i)= k1(i)*aux*aux/max(cy,em30)
3252 cfi = cf(i)*abs(h1(i))
3253 aux = two*cfi*cfi/max(ms(j1),em20)
3254 IF(aux>st1(i))THEN
3255 k1(i) =zero
3256 c1(i) =cfi
3257 st1(i)=aux
3258 ENDIF
3259 ENDIF
3260C
3261 j1=ix2(i)
3262 IF(ms(j1)==zero)THEN
3263 k2(i) =zero
3264 c2(i) =zero
3265 st2(i)=zero
3266 ELSE
3267 k2(i)=kt(i)*abs(h2(i))
3268 c2(i)=c(i)*abs(h2(i))
3269 cx =four*c2(i)*c2(i)
3270 cy =eight*ms(j1)*k2(i)
3271 aux = sqrt(cx+cy)+two*c2(i)
3272 st2(i)= k2(i)*aux*aux/max(cy,em30)
3273 cfi = cf(i)*abs(h2(i))
3274 aux = two*cfi*cfi/max(ms(j1),em20)
3275 IF(aux>st2(i))THEN
3276 k2(i) =zero
3277 c2(i) =cfi
3278 st2(i)=aux
3279 ENDIF
3280 ENDIF
3281C
3282 j1=ix3(i)
3283 IF(ms(j1)==zero)THEN
3284 k3(i) =zero
3285 c3(i) =zero
3286 st3(i)=zero
3287 ELSE
3288 k3(i)=kt(i)*abs(h3(i))
3289 c3(i)=c(i)*abs(h3(i))
3290 cx =four*c3(i)*c3(i)
3291 cy =eight*ms(j1)*k3(i)
3292 aux = sqrt(cx+cy)+two*c3(i)
3293 st3(i)= k3(i)*aux*aux/max(cy,em30)
3294 cfi = cf(i)*abs(h3(i))
3295 aux = two*cfi*cfi/max(ms(j1),em20)
3296 IF(aux>st3(i))THEN
3297 k3(i) =zero
3298 c3(i) =cfi
3299 st3(i)=aux
3300 ENDIF
3301 ENDIF
3302C
3303 j1=ix4(i)
3304 IF(ms(j1)==zero)THEN
3305 k4(i) =zero
3306 c4(i) =zero
3307 st4(i)=zero
3308 ELSE
3309 k4(i)=kt(i)*abs(h4(i))
3310 c4(i)=c(i)*abs(h4(i))
3311 cx =four*c4(i)*c4(i)
3312 cy =eight*ms(j1)*k4(i)
3313 aux = sqrt(cx+cy)+two*c4(i)
3314 st4(i)= k4(i)*aux*aux/max(cy,em30)
3315 cfi = cf(i)*abs(h4(i))
3316 aux = two*cfi*cfi/max(ms(j1),em20)
3317 IF(aux>st4(i))THEN
3318 k4(i) =zero
3319 c4(i) =cfi
3320 st4(i)=aux
3321 ENDIF
3322 ENDIF
3323 ENDDO
3324C
3325 ELSE
3326 DO i=1,jlt
3327 IF(viscffric(i)/=zero
3328 . .AND.(ivis2==0.OR.ivis2==1))THEN
3329 IF(pene(i) == zero)cycle
3330C C(I)=2.*C(I)
3331C
3332 IF(msi(i)==zero)THEN
3333 ks(i) =zero
3334 cs(i) =zero
3335 stv(i)=zero
3336 ELSE
3337 cx = four*c(i)*c(i)
3338 cy = eight*msi(i)*kt(i)
3339 aux = sqrt(cx+cy)+two*c(i)
3340 stv(i)= kt(i)*aux*aux/max(cy,em30)
3341 aux = two*cf(i)*cf(i)/max(msi(i),em20)
3342 IF(aux>stv(i))THEN
3343 ks(i) =zero
3344 cs(i) =cf(i)
3345 stv(i)=aux
3346 ELSE
3347 ks(i)= kt(i)
3348 cs(i) =c(i)
3349 ENDIF
3350 ENDIF
3351C
3352 j1=ix1(i)
3353 IF(ms(j1)==zero)THEN
3354 k1(i) =zero
3355 c1(i) =zero
3356 st1(i)=zero
3357 ELSE
3358 k1(i)=kt(i)*abs(h1(i))
3359 c1(i)=c(i)*abs(h1(i))
3360 cx =four*c1(i)*c1(i)
3361 cy =eight*ms(j1)*k1(i)
3362 aux = sqrt(cx+cy)+two*c1(i)
3363 st1(i)= k1(i)*aux*aux/max(cy,em30)
3364 cfi = cf(i)*abs(h1(i))
3365 aux = two*cfi*cfi/max(ms(j1),em20)
3366 IF(aux>st1(i))THEN
3367 k1(i) =zero
3368 c1(i) =cfi
3369 st1(i)=aux
3370 ENDIF
3371 ENDIF
3372C
3373 j1=ix2(i)
3374 IF(ms(j1)==zero)THEN
3375 k2(i) =zero
3376 c2(i) =zero
3377 st2(i)=zero
3378 ELSE
3379 k2(i)=kt(i)*abs(h2(i))
3380 c2(i)=c(i)*abs(h2(i))
3381 cx =four*c2(i)*c2(i)
3382 cy =eight*ms(j1)*k2(i)
3383 aux = sqrt(cx+cy)+two*c2(i)
3384 st2(i)= k2(i)*aux*aux/max(cy,em30)
3385 cfi = cf(i)*abs(h2(i))
3386 aux = two*cfi*cfi/max(ms(j1),em20)
3387 IF(aux>st2(i))THEN
3388 k2(i) =zero
3389 c2(i) =cfi
3390 st2(i)=aux
3391 ENDIF
3392 ENDIF
3393C
3394 j1=ix3(i)
3395 IF(ms(j1)==zero)THEN
3396 k3(i) =zero
3397 c3(i) =zero
3398 st3(i)=zero
3399 ELSE
3400 k3(i)=kt(i)*abs(h3(i))
3401 c3(i)=c(i)*abs(h3(i))
3402 cx =four*c3(i)*c3(i)
3403 cy =eight*ms(j1)*k3(i)
3404 aux = sqrt(cx+cy)+two*c3(i)
3405 st3(i)= k3(i)*aux*aux/max(cy,em30)
3406 cfi = cf(i)*abs(h3(i))
3407 aux = two*cfi*cfi/max(ms(j1),em20)
3408 IF(aux>st3(i))THEN
3409 k3(i) =zero
3410 c3(i) =cfi
3411 st3(i)=aux
3412 ENDIF
3413 ENDIF
3414C
3415 j1=ix4(i)
3416 IF(ms(j1)==zero)THEN
3417 k4(i) =zero
3418 c4(i) =zero
3419 st4(i)=zero
3420 ELSE
3421 k4(i)=kt(i)*abs(h4(i))
3422 c4(i)=c(i)*abs(h4(i))
3423 cx =four*c4(i)*c4(i)
3424 cy =eight*ms(j1)*k4(i)
3425 aux = sqrt(cx+cy)+two*c4(i)
3426 st4(i)= k4(i)*aux*aux/max(cy,em30)
3427 cfi = cf(i)*abs(h4(i))
3428 aux = two*cfi*cfi/max(ms(j1),em20)
3429 IF(aux>st4(i))THEN
3430 k4(i) =zero
3431 c4(i) =cfi
3432 st4(i)=aux
3433 ENDIF
3434 ENDIF
3435 ELSE
3436 IF(pene(i) == zero)cycle
3437 ks(i) =stif(i)
3438 cs(i) =zero
3439 stv(i)=ks(i)
3440 k1(i) =stif(i)*abs(h1(i))
3441 c1(i) =zero
3442 st1(i)=k1(i)
3443 k2(i) =stif(i)*abs(h2(i))
3444 c2(i) =zero
3445 st2(i)=k2(i)
3446 k3(i) =stif(i)*abs(h3(i))
3447 c3(i) =zero
3448 st3(i)=k3(i)
3449 k4(i) =stif(i)*abs(h4(i))
3450 c4(i) =zero
3451 st4(i)=k4(i)
3452 ENDIF
3453 ENDDO
3454 ENDIF
3455 ENDIF
3456C
3457C=======================================================================
3458C---------------------------------
3459 itag = 0
3460 DO i=1,jlt
3461 IF(pene(i) == zero)cycle
3462!!
3463 fx1(i)=fxi(i)*h1(i)
3464 fy1(i)=fyi(i)*h1(i)
3465 fz1(i)=fzi(i)*h1(i)
3466C
3467 fx2(i)=fxi(i)*h2(i)
3468 fy2(i)=fyi(i)*h2(i)
3469 fz2(i)=fzi(i)*h2(i)
3470C
3471 fx3(i)=fxi(i)*h3(i)
3472 fy3(i)=fyi(i)*h3(i)
3473 fz3(i)=fzi(i)*h3(i)
3474C
3475 fx4(i)=fxi(i)*h4(i)
3476 fy4(i)=fyi(i)*h4(i)
3477 fz4(i)=fzi(i)*h4(i)
3478C
3479 phi1(i) = zero
3480 phi2(i) = zero
3481 phi3(i) = zero
3482 phi4(i) = zero
3483 ENDDO
3484
3485
3486 IF(idtmins==2.OR.idtmins_int/=0)THEN
3487 dti=dt2t
3488 CALL i24sms2(jlt ,ix1 ,ix2 ,ix3 ,ix4 ,
3489 2 nsvg ,h1 ,h2 ,h3 ,h4 ,stif ,
3490 3 nin ,noint ,mskyi_sms, iskyi_sms,nsms ,
3491 4 kt ,c ,cf ,dtmini,dti ,
3492 5 irtse ,nsne ,is2se,is2pt ,t2main_sms ,
3493 6 t2fac_sms)
3494 IF(dti<dt2t)THEN
3495 dt2t = dti
3496 neltst = noint
3497 ityptst = 10
3498 ENDIF
3499 ENDIF
3500C
3501 IF(idtmins_int/=0)THEN
3502 stif(1:jlt)=zero
3503 END IF
3504C------------For /LOAD/PRESSURE tag nodes in contact-------------
3505 tagip(1:jlt) = 0
3506 IF(ninloadp > 0) THEN
3507 DO k = kloadpinter(nin)+1, kloadpinter(nin+1)
3508 pp = loadpinter(k)
3509 ppl = loadp_hyd_inter(pp)
3510 dgapload = dgaploadint(k)
3511 DO i=1,jlt
3512 gapp= gapv(i) + dgapload
3513 IF((pene(i)/=zero).OR.(dist(i) <= gapp.AND.dist(i) >= zero)) THEN
3514 tagip(i) = 1
3515 IF(pene(i)==zero) THEN
3516 ix1(i) = ixx(i,1)
3517 ix2(i) = ixx(i,2)
3518 ix3(i) = ixx(i,3)
3519 ix4(i) = ixx(i,4)
3520 ENDIF
3521 tagncont(ppl,ixx(i,1)) = 1
3522 tagncont(ppl,ixx(i,2)) = 1
3523 tagncont(ppl,ixx(i,3)) = 1
3524 tagncont(ppl,ixx(i,4)) = 1
3525 jg = nsvg(i)
3526 IF(jg>0) THEN
3527C SPMD : do same after reception of forces for remote nodes
3528 IF (jg <= numnod) THEN
3529 tagncont(ppl,jg) = 1
3530 ELSE
3531 ig = jg - numnod
3532 IF (ig > 0) THEN
3533 IF (is2se(1,ig) > 0) THEN
3534 ie = is2se(1,ig)
3535 ELSEIF(is2se(2,ig) > 0) THEN
3536 ie = is2se(2,ig)
3537 END IF
3538 DO j =1,4
3539 tagncont(ppl,irtse(j,ie)) = 1
3540 END DO
3541 ENDIF
3542 ENDIF
3543 ENDIF
3544 ENDIF
3545 ENDDO
3546 ENDDO
3547
3548 ENDIF
3549C
3550C
3551C spmd : identification des noeuds interf. utiles a envoyer
3552 IF (nspmd>1) THEN
3553ctmp+1 mic only
3554#include "mic_lockon.inc"
3555 DO i = 1,jlt
3556 nn = nsvg(i)
3557 IF(nn<0)THEN
3558C tag temporaire de NSVFI a -
3559 hh = h1(i)+h2(i)+h3(i)+h4(i)
3560 IF(hh /= zero.OR.tagip(i)==1)THEN
3561 IF(isedge_fi(nin)%P(-nn)==0)THEN
3562 nsvfi(nin)%P(-nn) = -abs(nsvfi(nin)%P(-nn))
3563 ELSE
3564 j1=irtse_fi(nin)%P(1,-nn)
3565 nsvfi(nin)%P(j1) = -abs(nsvfi(nin)%P(j1))
3566 j1=irtse_fi(nin)%P(2,-nn)
3567 nsvfi(nin)%P(j1) = -abs(nsvfi(nin)%P(j1))
3568 j1=irtse_fi(nin)%P(3,-nn)
3569 nsvfi(nin)%P(j1) = -abs(nsvfi(nin)%P(j1))
3570 j1=irtse_fi(nin)%P(4,-nn)
3571 nsvfi(nin)%P(j1) = -abs(nsvfi(nin)%P(j1))
3572 ENDIF
3573 ENDIF
3574 ENDIF
3575 ENDDO
3576ctmp+1 mic only
3577#include "mic_lockoff.inc"
3578 ENDIF
3579
3580C
3581 IF (impl_s ==0.AND.(inacti==-1.OR.igsti==-1)) stif(1:jlt)=max(stif(1:jlt),stifkt(1:jlt))
3582 IF(iparit==3)THEN
3583 stop 789
3584 ELSEIF(iparit==0)THEN
3585 CALL i24ass0(jlt ,ix1 ,ix2 ,ix3 ,ix4 ,
3586 2 nsvg ,h1 ,h2 ,h3 ,h4 ,stif ,
3587 3 fx1 ,fy1 ,fz1 ,fx2 ,fy2 ,fz2 ,
3588 4 fx3 ,fy3 ,fz3 ,fx4 ,fy4 ,fz4 ,
3589 5 fxi ,fyi ,fzi ,a ,stifn ,nin ,
3590 6 intth ,phi ,fthe ,phi1 , phi2 ,phi3 ,
3591 7 phi4 ,intply ,iply ,inod_pxfem ,
3592 8 irtse ,nsne ,is2se ,is2pt,jtask )
3593
3594 ELSE
3595 CALL i24ass2(jlt ,ix1 ,ix2 ,ix3 ,ix4 ,
3596 2 nsvg ,h1 ,h2 ,h3 ,h4 ,stif ,
3597 3 fx1 ,fy1 ,fz1 ,fx2 ,fy2 ,fz2 ,
3598 4 fx3 ,fy3 ,fz3 ,fx4 ,fy4 ,fz4 ,
3599 5 fxi ,fyi ,fzi ,fskyi,isky ,niskyfi,
3600 6 nin ,noint ,intth,phi ,ftheskyi ,phi1,
3601 7 phi2 ,phi3 , phi4 ,itab ,intply,iply ,
3602 8 inod_pxfem,irtse,nsne,is2se,is2pt,tagip )
3603 ENDIF
3604C
3605 IF(anim_v(4)+outp_v(4)+h3d_data%N_VECT_CONT>0)THEN
3606 IF (inconv==1) THEN
3607#include "lockon.inc"
3608 DO i=1,jlt
3609 IF(pene(i) == zero)cycle
3610 fcont(1,ix1(i)) =fcont(1,ix1(i)) + fx1(i)
3611 fcont(2,ix1(i)) =fcont(2,ix1(i)) + fy1(i)
3612 fcont(3,ix1(i)) =fcont(3,ix1(i)) + fz1(i)
3613 fcont(1,ix2(i)) =fcont(1,ix2(i)) + fx2(i)
3614 fcont(2,ix2(i)) =fcont(2,ix2(i)) + fy2(i)
3615 fcont(3,ix2(i)) =fcont(3,ix2(i)) + fz2(i)
3616 fcont(1,ix3(i)) =fcont(1,ix3(i)) + fx3(i)
3617 fcont(2,ix3(i)) =fcont(2,ix3(i)) + fy3(i)
3618 fcont(3,ix3(i)) =fcont(3,ix3(i)) + fz3(i)
3619 fcont(1,ix4(i)) =fcont(1,ix4(i)) + fx4(i)
3620 fcont(2,ix4(i)) =fcont(2,ix4(i)) + fy4(i)
3621 fcont(3,ix4(i)) =fcont(3,ix4(i)) + fz4(i)
3622 jg = nsvg(i)
3623 IF(jg>0) THEN
3624C en SPMD : traitement a refaire apres reception noeud remote si JG < 0
3625 IF (jg >numnod) THEN
3626 ig = jg - numnod
3627 CALL i24for1_fic(3 ,irtse ,nsne ,is2se ,is2pt ,
3628 + ig ,fxi(i) ,fyi(i) ,fzi(i) ,fcont ,
3629 + -1 )
3630 ELSE
3631 fcont(1,jg)=fcont(1,jg)- fxi(i)
3632 fcont(2,jg)=fcont(2,jg)- fyi(i)
3633 fcont(3,jg)=fcont(3,jg)- fzi(i)
3634 ENDIF
3635
3636 ELSE
3637! On QA test ../miniqa/INTERF/INT_24/implicit_Int24_Ipen0/data/CONVI7TU_0001.rad 2 x 2 JG = -1
3638C IF(ISEDGE_FI(NIN)%P(JG)==1)THEN
3639C ELSE
3640C ENDIF
3641
3642 ENDIF
3643 ENDDO
3644#include "lockoff.inc"
3645 END IF !(INCONV==1) THEN
3646 ENDIF
3647C-----------------------------------------------------
3648 IF(isecin>0.AND.inconv==1)THEN
3649 k0=nstrf(25)
3650 IF(nstrf(1)+nstrf(2)/=0)THEN
3651 DO i=1,nsect
3652 nbinter=nstrf(k0+14)
3653 k1s=k0+30
3654 DO j=1,nbinter
3655 IF(nstrf(k1s)==noint)THEN
3656 IF(isecut/=0)THEN
3657#include "lockon.inc"
3658 DO k=1,jlt
3659 IF(pene(k) == zero)cycle
3660C attention aux signes pour le cumul des efforts
3661C a rendre conforme avec CFORC3
3662 IF(secfcum(4,ix1(k),i)==1.)THEN
3663 secfcum(1,ix1(k),i)=secfcum(1,ix1(k),i)-fx1(k)
3664 secfcum(2,ix1(k),i)=secfcum(2,ix1(k),i)-fy1(k)
3665 secfcum(3,ix1(k),i)=secfcum(3,ix1(k),i)-fz1(k)
3666 ENDIF
3667 IF(secfcum(4,ix2(k),i)==1.)THEN
3668 secfcum(1,ix2(k),i)=secfcum(1,ix2(k),i)-fx2(k)
3669 secfcum(2,ix2(k),i)=secfcum(2,ix2(k),i)-fy2(k)
3670 secfcum(3,ix2(k),i)=secfcum(3,ix2(k),i)-fz2(k)
3671 ENDIF
3672 IF(secfcum(4,ix3(k),i)==1.)THEN
3673 secfcum(1,ix3(k),i)=secfcum(1,ix3(k),i)-fx3(k)
3674 secfcum(2,ix3(k),i)=secfcum(2,ix3(k),i)-fy3(k)
3675 secfcum(3,ix3(k),i)=secfcum(3,ix3(k),i)-fz3(k)
3676 ENDIF
3677 IF(secfcum(4,ix4(k),i)==1.)THEN
3678 secfcum(1,ix4(k),i)=secfcum(1,ix4(k),i)-fx4(k)
3679 secfcum(2,ix4(k),i)=secfcum(2,ix4(k),i)-fy4(k)
3680 secfcum(3,ix4(k),i)=secfcum(3,ix4(k),i)-fz4(k)
3681 ENDIF
3682 jg = nsvg(k)
3683 IF(jg>0) THEN
3684C en SPMD : traitement a refaire apres reception noeud remote si JG < 0
3685 IF (jg >numnod) THEN
3686 ig = jg - numnod
3687 IF(secfcum(4,jg,i)==1.)THEN
3688 CALL i24for1_fic(3 ,irtse ,nsne ,is2se ,is2pt ,
3689 + ig ,fxi(i) ,fyi(i) ,fzi(i) ,secfcum(1,1,i),
3690 + 1 )
3691 ENDIF
3692 ELSE
3693 IF(secfcum(4,jg,i)==1.)THEN
3694 secfcum(1,jg,i)=secfcum(1,jg,i)+fxi(k)
3695 secfcum(2,jg,i)=secfcum(2,jg,i)+fyi(k)
3696 secfcum(3,jg,i)=secfcum(3,jg,i)+fzi(k)
3697 ENDIF
3698 ENDIF
3699 ELSE
3700C Noeud Remote
3701 jg=-jg
3702 IF(isedge_fi(nin)%P(jg)==1)THEN
3703 IF(secfcum(4,jg,i)==1.)THEN
3704 CALL i24for1_fic(3 ,irtse_fi(nin)%P(1,1) ,nsne ,is2se_fi(nin)%P(1,1) ,is2pt_fi(nin)%P(1) ,
3705 + jg ,fxi(i) ,fyi(i) ,fzi(i) ,secfcum(1,1,i),
3706 + 1 )
3707 ENDIF
3708 ELSE
3709 IF(secfcum(4,jg,i)==1.)THEN
3710 secfcum(1,jg,i)=secfcum(1,jg,i)+fxi(k)
3711 secfcum(2,jg,i)=secfcum(2,jg,i)+fyi(k)
3712 secfcum(3,jg,i)=secfcum(3,jg,i)+fzi(k)
3713 ENDIF
3714 ENDIF
3715 ENDIF
3716 ENDDO
3717#include "lockoff.inc"
3718 ENDIF
3719C +fsav(section)
3720 ENDIF
3721 k1s=k1s+1
3722 ENDDO
3723 k0=nstrf(k0+24)
3724 ENDDO
3725 ENDIF
3726 ENDIF
3727C-----------------------------------------------------
3728 IF(ibag/=0.OR.iadm/=0)THEN
3729 DO i=1,jlt
3730 IF(pene(i) == zero)cycle
3731C test modifie pour coherence avec communication SPMD (spmd_i7tools)
3732c IF(FXI(I)/=ZERO.OR.FYI(I)/=ZERO.OR.FZI(I)/=ZERO)THEN
3733 jg = nsvg(i)
3734 IF(jg>0) THEN
3735C en SPMD : traitement a refaire apres reception noeud remote si JG < 0
3736 icontact(jg)=1
3737 ENDIF
3738 icontact(ix1(i))=1
3739 icontact(ix2(i))=1
3740 icontact(ix3(i))=1
3741 icontact(ix4(i))=1
3742 ENDDO
3743 ENDIF
3744C
3745 IF(iadm/=0)THEN
3746 END IF
3747 IF(iadm>=2)THEN
3748 END IF
3749C
3750 IF(ibcc==0) RETURN
3751C
3752 DO i=1,jlt
3753 IF(pene(i) == zero)cycle
3754 ibcm = ibcc / 8
3755 ibcs = ibcc - 8 * ibcm
3756 IF(ibcs>0) THEN
3757 ig=nsvg(i)
3758C---------obsolate option, no need to update
3759 IF(ig>0.AND.ig<=numnod) THEN
3760C en SPMD : traitement a refaire apres reception noeud remote si JG < 0
3761 CALL ibcoff(ibcs,icodt(ig))
3762 ENDIF
3763 ENDIF
3764 IF(ibcm>0) THEN
3765 ig=ix1(i)
3766 CALL ibcoff(ibcm,icodt(ig))
3767 ig=ix2(i)
3768 CALL ibcoff(ibcm,icodt(ig))
3769 ig=ix3(i)
3770 CALL ibcoff(ibcm,icodt(ig))
3771 ig=ix4(i)
3772 CALL ibcoff(ibcm,icodt(ig))
3773 ENDIF
3774 ENDDO
3775C
3776 RETURN
bitget
integer function bitget(i, n)
Definition bitget.F:37
alpha
#define alpha
Definition eval.h:35
area
subroutine area(d1, x, x2, y, y2, eint, stif0)
Definition hm_read_prop32.F:567
i24_save_sub
subroutine i24_save_sub(numnod, mvsiz, nisub, s_addsubm, s_lisubm, s_typsub, nisubmax, i_stok, ie, itypsub, nin, i, nn, nft, addsubm, lisubm, typsub, intarean, intcarea, isensint, fxi, fyi, fzi, fni, dt12, fsavsub1, fsavparit, nrtse, irtse, nsne, is2se, is2pt, nsnr)
Definition i24_save_sub.F:40
i24sms2
subroutine i24sms2(jlt, ix1, ix2, ix3, ix4, nsvg, h1, h2, h3, h4, stif, nin, noint, mskyi_sms, iskyi_sms, nsms, kt, c, cf, dtmini, dti, irtse, nsne, is2se, is2pt, t2main_sms, t2fac_sms)
Definition i24for3.F:4758
i24ass2
subroutine i24ass2(jlt, ix1, ix2, ix3, ix4, nsvg, h1, h2, h3, h4, stif, fx1, fy1, fz1, fx2, fy2, fz2, fx3, fy3, fz3, fx4, fy4, fz4, fxi, fyi, fzi, fskyi, isky, niskyfi, nin, noint, intth, phi, ftheskyi, phi1, phi2, phi3, phi4, itab, intply, iply, inod, irtse, nsne, is2se, is2pt, tagip)
Definition i24for3.F:3800
i24ass0
subroutine i24ass0(jlt, ix1, ix2, ix3, ix4, nsvg, h1, h2, h3, h4, stif, fx1, fy1, fz1, fx2, fy2, fz2, fx3, fy3, fz3, fx4, fy4, fz4, fxi, fyi, fzi, a, stifn, nin, intth, phi, fthe, phi1, phi2, phi3, phi4, intply, iply, inod, irtse, nsne, is2se, is2pt, jtask)
Definition i24for3.F:4434
i24for1_fic
subroutine i24for1_fic(npt, irtse, nsne, is2se, is2pt, ns, fxi, fyi, fzi, for1, inega)
Definition i24for3e.F:148
i24for1_ficr
subroutine i24for1_ficr(npt, irtse, nsne, is2se, is2pt, ns, fxi, fyi, fzi, for1, inega, ni)
Definition i24for3e.F:285
ibcoff
subroutine ibcoff(ibc, icodt)
Definition ibcoff.F:44
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
anim_mod
Definition anim_mod.F:96
h3d_mod
Definition h3d_mod.F:308
outputs_mod
Definition outputs_mod.F:59
parameters_mod
Definition parameters_mod.F:49
tri7box::stif_oldfi
type(real_pointer2), dimension(:), allocatable stif_oldfi
Definition tri7box.F:545
tri7box::secnd_frfi
type(real_pointer2), dimension(:), allocatable secnd_frfi
Definition tri7box.F:543
tri7box::inflg_subsfi
type(int_pointer), dimension(:), allocatable inflg_subsfi
Definition tri7box.F:505
tri7box::fnconti
type(real_pointer2), dimension(:), allocatable fnconti
Definition tri7box.F:510
tri7box::efricgfi
type(real_pointer), dimension(:), allocatable efricgfi
Definition tri7box.F:511
tri7box::lisubsfi
type(int_pointer), dimension(:), allocatable lisubsfi
Definition tri7box.F:501
tri7box::nsvfi
type(int_pointer), dimension(:), allocatable nsvfi
Definition tri7box.F:431
tri7box::intareanfi
type(real_pointer), dimension(:), allocatable intareanfi
Definition tri7box.F:554
tri7box::efricfi
type(real_pointer), dimension(:), allocatable efricfi
Definition tri7box.F:511
tri7box::addsubsfi
type(int_pointer), dimension(:), allocatable addsubsfi
Definition tri7box.F:509
tri7box::ftconti
type(real_pointer2), dimension(:), allocatable ftconti
Definition tri7box.F:510
tri7box::pene_oldfi
type(real_pointer2), dimension(:), allocatable pene_oldfi
Definition tri7box.F:544

◆ i24sms2()

subroutine i24sms2 ( integer jlt,
integer, dimension(mvsiz) ix1,
integer, dimension(mvsiz) ix2,
integer, dimension(mvsiz) ix3,
integer, dimension(mvsiz) ix4,
integer, dimension(mvsiz) nsvg,
h1,
h2,
h3,
h4,
stif,
integer nin,
integer noint,
mskyi_sms,
integer, dimension(lskyi_sms,*) iskyi_sms,
integer, dimension(*) nsms,
kt,
c,
cf,
dtmini,
dti,
integer, dimension(5,*) irtse,
integer nsne,
integer, dimension(2,*) is2se,
integer, dimension(*) is2pt,
integer, dimension(6,*) t2main_sms,
t2fac_sms )

Definition at line 4752 of file i24for3.F.

4758C-----------------------------------------------
4759C M o d u l e s
4760C-----------------------------------------------
4761 USE tri7box
4762 USE message_mod
4763C-----------------------------------------------
4764C I m p l i c i t T y p e s
4765C-----------------------------------------------
4766#include "implicit_f.inc"
4767#include "comlock.inc"
4768C-----------------------------------------------
4769C G l o b a l P a r a m e t e r s
4770C-----------------------------------------------
4771#include "mvsiz_p.inc"
4772C-----------------------------------------------
4773C C o m m o n B l o c k s
4774C-----------------------------------------------
4775#include "parit_c.inc"
4776#include "task_c.inc"
4777#include "sms_c.inc"
4778#include "com04_c.inc"
4779C-----------------------------------------------
4780C D u m m y A r g u m e n t s
4781C-----------------------------------------------
4782 INTEGER JLT,NIN,NOINT,
4783 . IX1(MVSIZ),IX2(MVSIZ),IX3(MVSIZ),IX4(MVSIZ),NSVG(MVSIZ),
4784 . NSMS(*), ISKYI_SMS(LSKYI_SMS,*),IRTSE(5,*),NSNE,IS2SE(2,*),
4785 . IS2PT(*),T2MAIN_SMS(6,*)
4786 my_real
4787 . h1(mvsiz),h2(mvsiz),h3(mvsiz),h4(mvsiz),stif(mvsiz),
4788 . mskyi_sms(*), kt(mvsiz), c(mvsiz), cf(mvsiz), dtmini, dti,t2fac_sms(*)
4789C-----------------------------------------------
4790C L o c a l V a r i a b l e s
4791C-----------------------------------------------
4792 INTEGER I, IG, NISKYL1, NISKYL, NN,JG,IXSS(4),NFIC,J,NS,NL, K,KK,IE
4793 my_real
4794 . hh(mvsiz),mas, dts,dtm_int,fici,fics(4),fac
4795
4796C-----------------------------------------------
4797C IF contact on nodes of TYPE2 additional connections must be created (T2MAIN_SMS(1) > 1)
4798C-----------------------------------------------
4799C
4800 niskyl1 = 0
4801C
4802 DO i=1,jlt
4803
4804 hh(i)=h1(i)+h2(i)+h3(i)+h4(i)
4805 IF(nsms(i)==0.OR.stif(i)==zero.OR.hh(i)==zero) cycle
4806C
4807 ig =nsvg(i)
4808 IF(ig > 0) THEN
4809 IF (ig > numnod) THEN
4810 nl = 0
4811 ns = ig - numnod
4812 fici = one
4813 CALL i24forc_fic(3,irtse,nsne,is2se,is2pt,ns,1,fici,fics,ixss)
4814 DO j=1,4
4815 nl = nl + t2main_sms(1,ixss(j))*ceiling(fics(j))
4816 ENDDO
4817 ELSE
4818 nl = t2main_sms(1,ig)
4819 ENDIF
4820 ELSE
4821 IF(isedge_fi(nin)%P(-ig)==1) THEN
4822 nl = 0
4823 fici = one
4824 CALL i24forc_fic(3,irtse_fi(nin)%P(1,1),nsne,is2se_fi(nin)%P(1,1), is2pt_fi(nin)%P(1),-ig,1,fici,fics,ixss)
4825 DO j=1,4
4826 nl = nl + t2main_sms_fi(nin)%P(1,ixss(j))*ceiling(fics(j))
4827 ENDDO
4828 ELSE
4829 nl = t2main_sms_fi(nin)%P(1,-ig)
4830 ENDIF
4831 ENDIF
4832C
4833 IF (h1(i)/=zero) niskyl1 = niskyl1 + nl*t2main_sms(1,ix1(i))
4834 IF (h2(i)/=zero) niskyl1 = niskyl1 + nl*t2main_sms(1,ix2(i))
4835 IF (h3(i)/=zero) niskyl1 = niskyl1 + nl*t2main_sms(1,ix3(i))
4836 IF (h4(i)/=zero) niskyl1 = niskyl1 + nl*t2main_sms(1,ix4(i))
4837C
4838 ENDDO
4839C
4840#include "lockon.inc"
4841 niskyl = nisky_sms
4842 nisky_sms = nisky_sms + niskyl1
4843#include "lockoff.inc"
4844C
4845 IF (niskyl+niskyl1 > lskyi_sms) THEN
4846 CALL ancmsg(msgid=26,anmode=aninfo)
4847 CALL arret(2)
4848 ENDIF
4849C
4850 IF (dtmini>zero) THEN
4851 dtm_int=dtmini
4852 ELSE
4853 dtm_int=dtmins_int
4854 ENDIF
4855
4856 DO i=1,jlt
4857 IF(nsms(i)==0.OR.stif(i)==zero.OR.hh(i)==zero) cycle
4858C
4859 IF(nsms(i)>0)THEN
4860 dts = dtmins/dtfacs
4861 dti=min(dti,dtmins)
4862 ELSE
4863 dts = dtm_int/dtfacs_int
4864 dti=min(dti,dtm_int)
4865 END IF
4866
4867 mas= dts * ( dts * kt(i) + c(i) )
4868 mas = half * max( mas, dts * cf(i) )
4869C
4870C-----------------------------------------------
4871C FAC = 1.0 for all cases
4872C FAC = 0.5 for type2 crossed connections - T2MAIN(i)*T2MAIN(j)=16
4873C-----------------------------------------------
4874C
4875 ig =nsvg(i)
4876 IF(ig > 0)THEN
4877 IF (ig > numnod) THEN
4878C-----------------------------------------------
4879C-- Edge 2 Edge contact
4880C-----------------------------------------------
4881 ns = ig - numnod
4882 fici = mas
4883 CALL i24forc_fic(3 ,irtse ,nsne ,is2se ,is2pt ,
4884 + ns ,1 ,fici ,fics ,ixss )
4885C
4886 DO j = 1, 4
4887 IF (fics(j) > zero) THEN
4888 IF(h1(i)/=zero)THEN
4889 fac = t2fac_sms(ixss(j))
4890 IF (t2main_sms(1,ixss(j))*t2main_sms(1,ix1(i))==16) fac = half*max(t2fac_sms(ixss(j)),t2fac_sms(ix1(i)))
4891 DO k=1,t2main_sms(1,ixss(j))
4892 DO kk=1,t2main_sms(1,ix1(i))
4893 niskyl=niskyl+1
4894 mskyi_sms(niskyl)=abs(h1(i))*fics(j)*fac
4895 iskyi_sms(niskyl,1)=t2main_sms(k+1,ixss(j))
4896 iskyi_sms(niskyl,2)=t2main_sms(kk+1,ix1(i))
4897 iskyi_sms(niskyl,3)=ispmd+1
4898 ENDDO
4899 ENDDO
4900 ENDIF
4901 IF(h2(i)/=zero)THEN
4902 fac = t2fac_sms(ixss(j))
4903 IF (t2main_sms(1,ixss(j))*t2main_sms(1,ix2(i))==16) fac = half*max(t2fac_sms(ixss(j)),t2fac_sms(ix2(i)))
4904 DO k=1,t2main_sms(1,ixss(j))
4905 DO kk=1,t2main_sms(1,ix2(i))
4906 niskyl=niskyl+1
4907 mskyi_sms(niskyl)=abs(h2(i))*fics(j)*fac
4908 iskyi_sms(niskyl,1)=t2main_sms(k+1,ixss(j))
4909 iskyi_sms(niskyl,2)=t2main_sms(kk+1,ix2(i))
4910 iskyi_sms(niskyl,3)=ispmd+1
4911 ENDDO
4912 ENDDO
4913 ENDIF
4914 IF(h3(i)/=zero)THEN
4915 fac = t2fac_sms(ixss(j))
4916 IF (t2main_sms(1,ixss(j))*t2main_sms(1,ix3(i))==16) fac = half*max(t2fac_sms(ixss(j)),t2fac_sms(ix3(i)))
4917 DO k=1,t2main_sms(1,ixss(j))
4918 DO kk=1,t2main_sms(1,ix3(i))
4919 niskyl=niskyl+1
4920 mskyi_sms(niskyl)=abs(h3(i))*fics(j)*fac
4921 iskyi_sms(niskyl,1)=t2main_sms(k+1,ixss(j))
4922 iskyi_sms(niskyl,2)=t2main_sms(kk+1,ix3(i))
4923 iskyi_sms(niskyl,3)=ispmd+1
4924 ENDDO
4925 ENDDO
4926 ENDIF
4927 IF(h4(i)/=zero)THEN
4928 fac = t2fac_sms(ixss(j))
4929 IF (t2main_sms(1,ixss(j))*t2main_sms(1,ix4(i))==16) fac = half*max(t2fac_sms(ixss(j)),t2fac_sms(ix4(i)))
4930 DO k=1,t2main_sms(1,ixss(j))
4931 DO kk=1,t2main_sms(1,ix4(i))
4932 niskyl=niskyl+1
4933 mskyi_sms(niskyl)=abs(h4(i))*fics(j)*fac
4934 iskyi_sms(niskyl,1)=t2main_sms(k+1,ixss(j))
4935 iskyi_sms(niskyl,2)=t2main_sms(kk+1,ix4(i))
4936 iskyi_sms(niskyl,3)=ispmd+1
4937 ENDDO
4938 ENDDO
4939 ENDIF
4940 ENDIF
4941 END DO
4942 ELSE
4943C-----------------------------------------------
4944C-- Main / node contact
4945C-----------------------------------------------
4946 IF(h1(i)/=zero)THEN
4947 fac = t2fac_sms(ig)
4948 IF (t2main_sms(1,ig)*t2main_sms(1,ix1(i))==16) fac = half*max(t2fac_sms(ig),t2fac_sms(ix1(i)))
4949 DO k=1,t2main_sms(1,ig)
4950 DO kk=1,t2main_sms(1,ix1(i))
4951 niskyl=niskyl+1
4952 mskyi_sms(niskyl)=abs(h1(i))*mas*fac
4953 iskyi_sms(niskyl,1)=t2main_sms(k+1,ig)
4954 iskyi_sms(niskyl,2)=t2main_sms(kk+1,ix1(i))
4955 iskyi_sms(niskyl,3)=ispmd+1
4956 ENDDO
4957 ENDDO
4958 END IF
4959 IF(h2(i)/=zero)THEN
4960 fac = t2fac_sms(ig)
4961 IF (t2main_sms(1,ig)*t2main_sms(1,ix2(i))==16) fac = half*max(t2fac_sms(ig),t2fac_sms(ix2(i)))
4962 DO k=1,t2main_sms(1,ig)
4963 DO kk=1,t2main_sms(1,ix2(i))
4964 niskyl=niskyl+1
4965 mskyi_sms(niskyl)=abs(h2(i))*mas*fac
4966 iskyi_sms(niskyl,1)=t2main_sms(k+1,ig)
4967 iskyi_sms(niskyl,2)=t2main_sms(kk+1,ix2(i))
4968 iskyi_sms(niskyl,3)=ispmd+1
4969 ENDDO
4970 ENDDO
4971 END IF
4972 IF(h3(i)/=zero)THEN
4973 fac = t2fac_sms(ig)
4974 IF (t2main_sms(1,ig)*t2main_sms(1,ix3(i))==16) fac = half*max(t2fac_sms(ig),t2fac_sms(ix3(i)))
4975 DO k=1,t2main_sms(1,ig)
4976 DO kk=1,t2main_sms(1,ix3(i))
4977 niskyl=niskyl+1
4978 mskyi_sms(niskyl)=abs(h3(i))*mas*fac
4979 iskyi_sms(niskyl,1)=t2main_sms(k+1,ig)
4980 iskyi_sms(niskyl,2)=t2main_sms(kk+1,ix3(i))
4981 iskyi_sms(niskyl,3)=ispmd+1
4982 ENDDO
4983 ENDDO
4984 END IF
4985 IF(h4(i)/=zero)THEN
4986 fac = t2fac_sms(ig)
4987 IF (t2main_sms(1,ig)*t2main_sms(1,ix4(i))==16) fac = half*max(t2fac_sms(ig),t2fac_sms(ix4(i)))
4988 DO k=1,t2main_sms(1,ig)
4989 DO kk=1,t2main_sms(1,ix4(i))
4990 niskyl=niskyl+1
4991 mskyi_sms(niskyl)=abs(h4(i))*mas*fac
4992 iskyi_sms(niskyl,1)=t2main_sms(k+1,ig)
4993 iskyi_sms(niskyl,2)=t2main_sms(kk+1,ix4(i))
4994 iskyi_sms(niskyl,3)=ispmd+1
4995 ENDDO
4996 ENDDO
4997 END IF
4998 END IF !(IG > NUMNOD) THNE
4999C
5000C-----------------------------------------------
5001 ELSE
5002C-----------------------------------------------
5003C
5004 nn = -ig
5005 IF (isedge_fi(nin)%P(nn)==1)THEN
5006C-----------------------------------------------
5007C-- E2E remote contact
5008C-----------------------------------------------
5009 fici = mas
5010 CALL i24forc_fic(3 ,irtse_fi(nin)%P(1,1) ,nsne ,is2se_fi(nin)%P(1,1) , is2pt_fi(nin)%P(1) ,
5011 + nn ,1 ,fici ,fics ,ixss )
5012C
5013 DO j = 1, 4
5014 IF (fics(j) > zero) THEN
5015 IF(h1(i)/=zero)THEN
5016 fac = t2fac_sms_fi(nin)%P(ixss(j))
5017 IF (t2main_sms_fi(nin)%P(1,ixss(j))*t2main_sms(1,ix1(i))==16)
5018 . fac = half*max(t2fac_sms_fi(nin)%P(ixss(j)),t2fac_sms(ix1(i)))
5019 DO k=1,t2main_sms_fi(nin)%P(1,ixss(j))
5020 DO kk=1,t2main_sms(1,ix1(i))
5021 niskyl=niskyl+1
5022 mskyi_sms(niskyl)=abs(h1(i))*fics(j)*fac
5023 iskyi_sms(niskyl,1)=t2main_sms_fi(nin)%P(k+1,ixss(j))
5024 iskyi_sms(niskyl,2)=t2main_sms(kk+1,ix1(i))
5025 iskyi_sms(niskyl,3)=procamsfi(nin)%P(ixss(j))
5026 ENDDO
5027 ENDDO
5028 END IF
5029 IF(h2(i)/=zero)THEN
5030 fac = t2fac_sms_fi(nin)%P(ixss(j))
5031 IF (t2main_sms_fi(nin)%P(1,ixss(j))*t2main_sms(1,ix2(i))==16)
5032 . fac = half*max(t2fac_sms_fi(nin)%P(ixss(j)),t2fac_sms(ix2(i)))
5033 DO k=1,t2main_sms_fi(nin)%P(1,ixss(j))
5034 DO kk=1,t2main_sms(1,ix2(i))
5035 niskyl=niskyl+1
5036 mskyi_sms(niskyl)=abs(h2(i))*fics(j)*fac
5037 iskyi_sms(niskyl,1)=t2main_sms_fi(nin)%P(k+1,ixss(j))
5038 iskyi_sms(niskyl,2)=t2main_sms(kk+1,ix2(i))
5039 iskyi_sms(niskyl,3)=procamsfi(nin)%P(ixss(j))
5040 ENDDO
5041 ENDDO
5042 END IF
5043 IF(h3(i)/=zero)THEN
5044 fac = t2fac_sms_fi(nin)%P(ixss(j))
5045 IF (t2main_sms_fi(nin)%P(1,ixss(j))*t2main_sms(1,ix3(i))==16)
5046 . fac = half*max(t2fac_sms_fi(nin)%P(ixss(j)),t2fac_sms(ix3(i)))
5047 DO k=1,t2main_sms_fi(nin)%P(1,ixss(j))
5048 DO kk=1,t2main_sms(1,ix3(i))
5049 niskyl=niskyl+1
5050 mskyi_sms(niskyl)=abs(h3(i))*fics(j)*fac
5051 iskyi_sms(niskyl,1)=t2main_sms_fi(nin)%P(k+1,ixss(j))
5052 iskyi_sms(niskyl,2)=t2main_sms(kk+1,ix3(i))
5053 iskyi_sms(niskyl,3)=procamsfi(nin)%P(ixss(j))
5054 ENDDO
5055 ENDDO
5056 END IF
5057 IF(h4(i)/=zero)THEN
5058 fac = t2fac_sms_fi(nin)%P(ixss(j))
5059 IF (t2main_sms_fi(nin)%P(1,ixss(j))*t2main_sms(1,ix4(i))==16)
5060 . fac = half*max(t2fac_sms_fi(nin)%P(ixss(j)),t2fac_sms(ix4(i)))
5061 DO k=1,t2main_sms_fi(nin)%P(1,ixss(j))
5062 DO kk=1,t2main_sms(1,ix4(i))
5063 niskyl=niskyl+1
5064 mskyi_sms(niskyl)=abs(h4(i))*fics(j)*fac
5065 iskyi_sms(niskyl,1)=t2main_sms_fi(nin)%P(k+1,ixss(j))
5066 iskyi_sms(niskyl,2)=t2main_sms(kk+1,ix4(i))
5067 iskyi_sms(niskyl,3)=procamsfi(nin)%P(ixss(j))
5068 ENDDO
5069 ENDDO
5070 END IF
5071 ENDIF
5072 ENDDO
5073 ELSE
5074C-----------------------------------------------
5075C-- Main / node remote contact
5076C-----------------------------------------------
5077 IF(h1(i)/=zero)THEN
5078 fac = t2fac_sms_fi(nin)%P(nn)
5079 IF (t2main_sms_fi(nin)%P(1,nn)*t2main_sms(1,ix1(i))==16)
5080 . fac = half*max(t2fac_sms_fi(nin)%P(nn),t2fac_sms(ix1(i)))
5081 DO k=1,t2main_sms_fi(nin)%P(1,nn)
5082 DO kk=1,t2main_sms(1,ix1(i))
5083 niskyl=niskyl+1
5084 mskyi_sms(niskyl)=abs(h1(i))*mas*fac
5085 iskyi_sms(niskyl,1)=t2main_sms_fi(nin)%P(k+1,nn)
5086 iskyi_sms(niskyl,2)=t2main_sms(kk+1,ix1(i))
5087 iskyi_sms(niskyl,3)=procamsfi(nin)%P(nn)
5088 ENDDO
5089 ENDDO
5090 END IF
5091 IF(h2(i)/=zero)THEN
5092 fac = t2fac_sms_fi(nin)%P(nn)
5093 IF (t2main_sms_fi(nin)%P(1,nn)*t2main_sms(1,ix2(i))==16)
5094 . fac = half*max(t2fac_sms_fi(nin)%P(nn),t2fac_sms(ix2(i)))
5095 DO k=1,t2main_sms_fi(nin)%P(1,nn)
5096 DO kk=1,t2main_sms(1,ix2(i))
5097 niskyl=niskyl+1
5098 mskyi_sms(niskyl)=abs(h2(i))*mas*fac
5099 iskyi_sms(niskyl,1)=t2main_sms_fi(nin)%P(k+1,nn)
5100 iskyi_sms(niskyl,2)=t2main_sms(kk+1,ix2(i))
5101 iskyi_sms(niskyl,3)=procamsfi(nin)%P(nn)
5102 ENDDO
5103 ENDDO
5104 END IF
5105 IF(h3(i)/=zero)THEN
5106 fac = t2fac_sms_fi(nin)%P(nn)
5107 IF (t2main_sms_fi(nin)%P(1,nn)*t2main_sms(1,ix3(i))==16)
5108 . fac = half*max(t2fac_sms_fi(nin)%P(nn),t2fac_sms(ix3(i)))
5109 DO k=1,t2main_sms_fi(nin)%P(1,nn)
5110 DO kk=1,t2main_sms(1,ix3(i))
5111 niskyl=niskyl+1
5112 mskyi_sms(niskyl)=abs(h3(i))*mas*fac
5113 iskyi_sms(niskyl,1)=t2main_sms_fi(nin)%P(k+1,nn)
5114 iskyi_sms(niskyl,2)=t2main_sms(kk+1,ix3(i))
5115 iskyi_sms(niskyl,3)=procamsfi(nin)%P(nn)
5116 ENDDO
5117 ENDDO
5118 END IF
5119 IF(h4(i)/=zero)THEN
5120 fac = t2fac_sms_fi(nin)%P(nn)
5121 IF (t2main_sms_fi(nin)%P(1,nn)*t2main_sms(1,ix4(i))==16)
5122 . fac = half*max(t2fac_sms_fi(nin)%P(nn),t2fac_sms(ix4(i)))
5123 DO k=1,t2main_sms_fi(nin)%P(1,nn)
5124 DO kk=1,t2main_sms(1,ix4(i))
5125 niskyl=niskyl+1
5126 mskyi_sms(niskyl)=abs(h4(i))*mas*fac
5127 iskyi_sms(niskyl,1)=t2main_sms_fi(nin)%P(k+1,nn)
5128 iskyi_sms(niskyl,2)=t2main_sms(kk+1,ix4(i))
5129 iskyi_sms(niskyl,3)=procamsfi(nin)%P(nn)
5130 ENDDO
5131 ENDDO
5132 END IF
5133 ENDIF
5134 END IF
5135 ENDDO
5136C
5137 RETURN
tri7box::t2fac_sms_fi
type(real_pointer), dimension(:), allocatable t2fac_sms_fi
Definition tri7box.F:557
tri7box::procamsfi
type(int_pointer), dimension(:), allocatable procamsfi
Definition tri7box.F:440
tri7box::t2main_sms_fi
type(int_pointer2), dimension(:), allocatable t2main_sms_fi
Definition tri7box.F:558
nl
character *2 function nl()
Definition message.F:2354