i2for28_cin.F

Go to the documentation of this file.

Copyright>        OpenRadioss
Copyright>        Copyright (C) 1986-2025 Altair Engineering Inc.
Copyright>
Copyright>        This program is free software: you can redistribute it and/or modify
Copyright>        it under the terms of the GNU Affero General Public License as published by
Copyright>        the Free Software Foundation, either version 3 of the License, or
Copyright>        (at your option) any later version.
Copyright>
Copyright>        This program is distributed in the hope that it will be useful,
Copyright>        but WITHOUT ANY WARRANTY; without even the implied warranty of
Copyright>        MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
Copyright>        GNU Affero General Public License for more details.
Copyright>
Copyright>        You should have received a copy of the GNU Affero General Public License
Copyright>        along with this program.  If not, see <https://www.gnu.org/licenses/>.
Copyright>
Copyright>
Copyright>        Commercial Alternative: Altair Radioss Software
Copyright>
Copyright>        As an alternative to this open-source version, Altair also offers Altair Radioss
Copyright>        software under a commercial license.  Contact Altair to discuss further if the
Copyright>        commercial version may interest you: https://www.altair.com/radioss/.
!||====================================================================
!||    i2for28_cin   ../engine/source/interfaces/interf/i2for28_cin.F
!||--- called by ------------------------------------------------------
!||    i2for28       ../engine/source/interfaces/interf/i2for28.F
!||--- calls      -----------------------------------------------------
!||    i2forces      ../engine/source/interfaces/interf/i2forces.F
!||--- uses       -----------------------------------------------------
!||    h3d_mod       ../engine/share/modules/h3d_mod.F
!||====================================================================
      SUBROUTINE i2for28_cin(NSN     ,NMN     ,A      ,IRECT  ,DPARA   ,
     2                       MSR     ,NSV     ,IRTL   ,MS     ,WEIGHT  ,
     3                       AR      ,IN      ,X      ,STIFN  ,STIFR   ,
     4                       FSAV    ,DMAST   ,ADM    ,MMASS  ,IDEL2   ,
     5                       SMASS   ,SINER   ,CRST   ,FNCONT ,INDXC   ,
     6                       MINER   ,H3D_DATA,FNCONTP,FTCONTP )
C-----------------------------------------------
C   M o d u l e s
C-----------------------------------------------
      USE h3d_mod
C-----------------------------------------------
C   I m p l i c i t   T y p e s
C-----------------------------------------------
#include      "implicit_f.inc"
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      INTEGER NSN, NMN,
     .   IRECT(4,*), MSR(*), NSV(*), IRTL(*), WEIGHT(*), IDEL2, INDXC(*)
C     REAL
      my_real
     .   X(3,*),A(3,*),AR(3,*), MMASS(*), CRST(2,*),
     .   dpara(7,*), ms(*), in(*),stifn(*),stifr(*),dmast,adm(*),
     .   smass(*), siner(*),fsav(*), fncont(3,*), miner(*),fncontp(3,*),
     .   ftcontp(3,*) 
      TYPE (H3D_DATABASE) :: H3D_DATA
C-----------------------------------------------
C   C o m m o n   B l o c k s
C-----------------------------------------------
#include      "com01_c.inc"
#include      "scr14_c.inc"
#include      "scr16_c.inc"
#include      "impl1_c.inc"
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER NIR, I, J, K, J1,J2,J3,J4, II, L, JJ
C     REAL
      my_real
     .   H(4),XM(4),YM(4),ZM(4),
     .   XMSJ, SS, ST, XMSI, FS(3),SP,SM,TP,TM,FACM,
     .   mx,my,mz,det,fx0,fy0,fz0,ins,
     .   x0,x1,x2,x3,x4,xs,y0,y1,y2,y3,y4,ys,z0,z1,z2,z3,z4,zs,
     .   x12,x22,x32,x42,y12,y22,y32,y42,z12,z22,z32,z42,
     .   xx,yy,zz,xxx,yyy,zzz,xy,yz,zx,xy2,yz2,zx2,s,t,    
     .   a1,a2,a3,b1,b2,b3,c1,c2,c3,mr,mrx,mry,mrz,inx,iny,inz,stf,fact,
     .   fx(4),fy(4),fz(4)
C=======================================================================
C
C------------------------------
C     DUPLICATED FROM I2FOMO3 - SPOTFLAG1 FORMULATION
C------------------------------
C
      nir=4
C
C---  Traitement specifique pour DMAS
C
C---  MMASS(II) initialise a MS(J) a t=0 dans starter
      IF(anim_n(2)+outp_n(2)+h3d_data%N_SCAL_DMAS >0) THEN
        DO ii=1,nmn
          j=msr(ii)
          adm(j) = adm(j)*mmass(ii)
        ENDDO
      ENDIF
C
C------------------------------
C     FORCES ET MOMENTS DES NOEUDS SECONDS 
C     TRANSMIS AUX NOEUDS MAINS SOUS 
C     FORME DE FORCES 
C
C     MASSES ET INERTIES DES NOEUDS SECONDS 
C     TRANSMISES AUX NOEUDS MAINS SOUS 
C     FORME DE MASSES 
C------------------------------
C
      IF(impl_s>0) THEN
      DO ii=1,nsn
       k = indxc(ii)
       IF (k == 0) cycle
       i = nsv(k)
       IF(i>0)THEN
        l=irtl(ii)
C
        s = crst(1,ii)  
        t = crst(2,ii)  
        sp=one + s          
        sm=one - s          
        tp=fourth*(one + t)  
        tm=fourth*(one - t)  
C    
        h(1)=one/nir                                         
        h(2)=one/nir                                          
        h(3)=one/nir                                          
        h(4)=one/nir                                         
C
        nir=4
        DO jj=1,nir
         j=irect(jj,l)
         xm(jj)=x(1,j)
         ym(jj)=x(2,j)
         zm(jj)=x(3,j)
        ENDDO 
        IF(irect(3,l)==irect(4,l)) THEN
         nir=3
         xm(4)=zero
         ym(4)=zero
         zm(4)=zero
        ENDIF
        facm = one / nir
C----------------------------------------------------
C       VITESSE DE ROTATION MOYENNE DU SEGMENT MAIN
C----------------------------------------------------
        x0=facm*(xm(1)+xm(2)+xm(3)+xm(4))
        y0=facm*(ym(1)+ym(2)+ym(3)+ym(4))
        z0=facm*(zm(1)+zm(2)+zm(3)+zm(4))
        DO j=1,nir
         xm(j)=xm(j)-x0
         ym(j)=ym(j)-y0
         zm(j)=zm(j)-z0
        ENDDO 
        xs=x(1,i)-x0
        ys=x(2,i)-y0
        zs=x(3,i)-z0
        xx=0 
        yy=0 
        zz=0 
        xy=0 
        yz=0 
        zx=0
        DO j=1,nir
          xx=xx+ xm(j)*xm(j)
          yy=yy+ ym(j)*ym(j)
          zz=zz+ zm(j)*zm(j)
          xy=xy+ xm(j)*ym(j)
          yz=yz+ ym(j)*zm(j)
          zx=zx+ zm(j)*xm(j)
        ENDDO 
        zzz=xx+yy
        xxx=yy+zz
        yyy=zz+xx 
        xy2=xy*xy
        yz2=yz*yz
        zx2=zx*zx
        det= xxx*yyy*zzz - xxx*yz2 - yyy*zx2 - zzz*xy2 - two*xy*yz*zx
        det=one/det
        b1=zzz*yyy-yz2
        b2=xxx*zzz-zx2
        b3=yyy*xxx-xy2
        c3=zzz*xy+yz*zx
        c1=xxx*yz+zx*xy
        c2=yyy*zx+xy*yz
C
        dpara(1,ii)=det
        dpara(2,ii)=b1
        dpara(3,ii)=b2
        dpara(4,ii)=b3
        dpara(5,ii)=c1
        dpara(6,ii)=c2
        dpara(7,ii)=c3
C
        xmsi=ms(i)*weight(i)
        fs(1)=a(1,i)*weight(i)
        fs(2)=a(2,i)*weight(i)
        fs(3)=a(3,i)*weight(i)
C
        IF (iroddl==1) THEN
          ins=in(i)*weight(i)
          mx=ar(1,i)*weight(i) + ys*fs(3) - zs*fs(2)
          my=ar(2,i)*weight(i) + zs*fs(1) - xs*fs(3)
          mz=ar(3,i)*weight(i) + xs*fs(2) - ys*fs(1)
        ELSE
          ins=zero
          mx=ys*fs(3) - zs*fs(2)
          my=zs*fs(1) - xs*fs(3)
          mz=xs*fs(2) - ys*fs(1)
        ENDIF
C
        a1=det*(mx*b1+my*c3+mz*c2)
        a2=det*(my*b2+mz*c1+mx*c3)
        a3=det*(mz*b3+mx*c2+my*c1)
C
        fx0=fs(1)*facm
        fy0=fs(2)*facm
        fz0=fs(3)*facm
C------------------------------------------------------
C     FORCES TRANSMISES AUX NOEUDS MAINS
C------------------------------------------------------
        fx(1:4) = zero
        fy(1:4) = zero
        fz(1:4) = zero
        DO jj=1,nir
          j=irect(jj,l)
          fx(jj)=fx0 + a2*zm(jj) - a3*ym(jj)
          fy(jj)=fy0 + a3*xm(jj) - a1*zm(jj)
          fz(jj)=fz0 + a1*ym(jj) - a2*xm(jj)
          a(1,j)=a(1,j) + fx(jj)
          a(2,j)=a(2,j) + fy(jj)
          a(3,j)=a(3,j) + fz(jj)
        ENDDO 
C------------------------------------------------------
C     INERTIES => MASSES
C------------------------------------------------------
C
        inx=ins + xmsi*(xs*xs+ys*ys+zs*zs)
        mrx = (b1+c3+c2)
        mry = (b2+c1+c3)
        mrz = (b3+c2+c1)
        mr=det*inx*max(mrx,mry,mrz)
C
C------------------------------------------------------
C     MASSES TRANSMISES AUX NOEUDS MAINS
C------------------------------------------------------
        xmsi=max(facm*xmsi,mr)
        dmast = dmast + nir*xmsi - ms(i)
        IF (anim_n(2)+outp_n(2)+h3d_data%N_SCAL_DMAS >0) THEN
         DO jj=1,nir
          j=irect(jj,l)
          adm(j) = adm(j) + xmsi - facm*ms(i)
         ENDDO 
        ENDIF
C
        IF (iroddl == 1) THEN
          stf = (facm*stifn(i)
     .         + det*max(mrx,mry,mrz)*(stifr(i)+stifn(i)*(xs*xs+ys*ys+zs*zs)))*weight(i)
        ELSE
          stf = (facm*stifn(i)
     .         + det*max(mrx,mry,mrz)*(stifn(i)*(xs*xs+ys*ys+zs*zs)))*weight(i)
        ENDIF
C
        DO jj=1,nir
         j=irect(jj,l)
         ms(j)=ms(j)+xmsi
         stifn(j)=stifn(j) + stf
        ENDDO 
C
        IF(idel2/=0.AND.ms(i)/=0.)smass(ii)=ms(i)
        ms(i)=zero
        stifn(i)=em20
C  
        IF (iroddl==1) THEN
          IF(idel2/=0.AND.in(i)/=0.)siner(ii)=in(i)
          in(i)=zero
          stifr(i)=em20
        ENDIF
C
       ENDIF
C----
       IF(i>0)THEN
C---     output of tied contact forces       
         CALL i2forces(x       ,fs     ,fx     ,fy      ,fz     ,
     .                irect(1,l),nir  ,fsav   ,fncont  ,fncontp,
     .                ftcontp ,weight ,h3d_data,i      ,h)              
       ENDIF                            
C----
      ENDDO
      
      ELSE
C
      DO ii=1,nsn
       k = indxc(ii)
       IF (k == 0) cycle
       i = nsv(k)
C
       IF(i>0)THEN
        l=irtl(ii)
C
        s = crst(1,ii)  
        t = crst(2,ii)  
        sp=one + s          
        sm=one - s          
        tp=fourth*(one + t)  
        tm=fourth*(one - t)  
C
        h(1)=one/nir                                         
        h(2)=one/nir                                          
        h(3)=one/nir                                          
        h(4)=one/nir                                           
C
        j1=irect(1,l)
        j2=irect(2,l)
        j3=irect(3,l)
        j4=irect(4,l)
        x1=x(1,j1)
        y1=x(2,j1)
        z1=x(3,j1)
        x2=x(1,j2)
        y2=x(2,j2)
        z2=x(3,j2)
        x3=x(1,j3)
        y3=x(2,j3)
        z3=x(3,j3)
        x4=x(1,j4)
        y4=x(2,j4)
        z4=x(3,j4)
        x0=fourth*(x1+x2+x3+x4)
        y0=fourth*(y1+y2+y3+y4)
        z0=fourth*(z1+z2+z3+z4)
        x1=x1-x0
        y1=y1-y0
        z1=z1-z0
        x2=x2-x0
        y2=y2-y0
        z2=z2-z0
        x3=x3-x0
        y3=y3-y0
        z3=z3-z0
        x4=x4-x0
        y4=y4-y0
        z4=z4-z0
        xs=x(1,i)-x0
        ys=x(2,i)-y0
        zs=x(3,i)-z0
C
        x12=x1*x1
        x22=x2*x2
        x32=x3*x3
        x42=x4*x4 
        y12=y1*y1
        y22=y2*y2
        y32=y3*y3
        y42=y4*y4 
        z12=z1*z1 
        z22=z2*z2
        z32=z3*z3 
        z42=z4*z4 
        xx=x12 + x22 + x32 + x42 
        yy=y12 + y22 + y32 + y42 
        zz=z12 + z22 + z32 + z42 
        xy=x1*y1 + x2*y2 + x3*y3 + x4*y4 
        yz=y1*z1 + y2*z2 + y3*z3 + y4*z4 
        zx=z1*x1 + z2*x2 + z3*x3 + z4*x4
        zzz=xx+yy
        xxx=yy+zz
        yyy=zz+xx 
        xy2=xy*xy
        yz2=yz*yz
        zx2=zx*zx
        det= xxx*yyy*zzz - xxx*yz2 - yyy*zx2 - zzz*xy2 - two*xy*yz*zx
        det=one/det
        b1=zzz*yyy-yz2
        b2=xxx*zzz-zx2
        b3=yyy*xxx-xy2
        c3=zzz*xy+yz*zx
        c1=xxx*yz+zx*xy
        c2=yyy*zx+xy*yz
C
        dpara(1,ii)=det
        dpara(2,ii)=b1
        dpara(3,ii)=b2
        dpara(4,ii)=b3
        dpara(5,ii)=c1
        dpara(6,ii)=c2
        dpara(7,ii)=c3
C
        xmsi=ms(i)*weight(i)
        fs(1)=a(1,i)*weight(i)
        fs(2)=a(2,i)*weight(i)
        fs(3)=a(3,i)*weight(i)
        IF (iroddl==1) THEN
          ins=in(i)*weight(i)
          mx=ar(1,i)*weight(i) + ys*fs(3) - zs*fs(2)
          my=ar(2,i)*weight(i) + zs*fs(1) - xs*fs(3)
          mz=ar(3,i)*weight(i) + xs*fs(2) - ys*fs(1)
        ELSE
          ins=zero
          mx=ys*fs(3) - zs*fs(2)
          my=zs*fs(1) - xs*fs(3)
          mz=xs*fs(2) - ys*fs(1)
        ENDIF
C
        a1=det*(mx*b1+my*c3+mz*c2)
        a2=det*(my*b2+mz*c1+mx*c3)
        a3=det*(mz*b3+mx*c2+my*c1)
C
        fx0=fs(1)*fourth
        fy0=fs(2)*fourth
        fz0=fs(3)*fourth
C------------------------------------------------------
C     FORCES TRANSMISES AUX NOEUDS MAINS
C------------------------------------------------------
        fx(1) = fx0 + a2*z1 - a3*y1
        fy(1) = fy0 + a3*x1 - a1*z1
        fz(1) = fz0 + a1*y1 - a2*x1
        fx(2) = fx0 + a2*z2 - a3*y2
        fy(2) = fy0 + a3*x2 - a1*z2
        fz(2) = fz0 + a1*y2 - a2*x2
        fx(3) = fx0 + a2*z3 - a3*y3
        fy(3) = fy0 + a3*x3 - a1*z3
        fz(3) = fz0 + a1*y3 - a2*x3
        fx(4) = fx0 + a2*z4 - a3*y4
        fy(4) = fy0 + a3*x4 - a1*z4
        fz(4) = fz0 + a1*y4 - a2*x4
c
        a(1,j1)=a(1,j1) + fx(1)
        a(2,j1)=a(2,j1) + fy(1)
        a(3,j1)=a(3,j1) + fz(1)
        a(1,j2)=a(1,j2) + fx(2)
        a(2,j2)=a(2,j2) + fy(2)
        a(3,j2)=a(3,j2) + fz(2)
        a(1,j3)=a(1,j3) + fx(3)
        a(2,j3)=a(2,j3) + fy(3)
        a(3,j3)=a(3,j3) + fz(3)
        a(1,j4)=a(1,j4) + fx(4)
        a(2,j4)=a(2,j4) + fy(4)
        a(3,j4)=a(3,j4) + fz(4)
C------------------------------------------------------
C     INERTIES => MASSES
C------------------------------------------------------
C
        inx=ins + xmsi*(xs*xs+ys*ys+zs*zs)
        mrx = (b1+c3+c2)
        mry = (b2+c1+c3)
        mrz = (b3+c2+c1)
        mr=det*inx*max(mrx,mry,mrz)
C
C------------------------------------------------------
C     MASSES TRANSMISES AUX NOEUDS MAINS
C------------------------------------------------------
C
        fact = one
        IF (iroddl==1) THEN
        IF (miner(j1)>zero.AND.miner(j2)>zero.AND.miner(j3)>zero.AND.miner(j4)>zero) THEN
C--       Inertie transmise sous forme d'inertie
          fact = zero
        ENDIF
        ENDIF
C
        xmsi=fourth*xmsi+mr*fact
C
        dmast = dmast + four*xmsi - ms(i)
        IF (anim_n(2)+outp_n(2)+h3d_data%N_SCAL_DMAS >0) THEN
         adm(j1) = adm(j1) + xmsi - fourth*ms(i)
         adm(j2) = adm(j2) + xmsi - fourth*ms(i)
         adm(j3) = adm(j3) + xmsi - fourth*ms(i)
         adm(j4) = adm(j4) + xmsi - fourth*ms(i)
        ENDIF
C
        ms(j1)=ms(j1)+xmsi
        ms(j2)=ms(j2)+xmsi
        ms(j3)=ms(j3)+xmsi
        ms(j4)=ms(j4)+xmsi
C
        IF (iroddl == 1) THEN
          stf = (fourth*stifn(i)
     .         + det*max(mrx,mry,mrz)*(stifr(i)+stifn(i)*(xs*xs+ys*ys+zs*zs)))*weight(i)
        ELSE
          stf = (fourth*stifn(i)
     .         + det*max(mrx,mry,mrz)*(stifn(i)*(xs*xs+ys*ys+zs*zs)))*weight(i)
        ENDIF
C
        stifn(j1)=stifn(j1) + stf
        stifn(j2)=stifn(j2) + stf
        stifn(j3)=stifn(j3) + stf
        stifn(j4)=stifn(j4) + stf
C
        IF (iroddl==1) THEN
          in(j1)=in(j1)+inx*fourth*(one-fact)
          in(j2)=in(j2)+inx*fourth*(one-fact)
          in(j3)=in(j3)+inx*fourth*(one-fact)
          in(j4)=in(j4)+inx*fourth*(one-fact)
        ENDIF
C
        IF(idel2/=0.AND.ms(i)/=0.)smass(ii)=ms(i)
        ms(i)=zero
        stifn(i)=em20
C  
        IF (iroddl==1) THEN
          IF(idel2/=0.AND.in(i)/=0.)siner(ii)=in(i)
          in(i)=zero
          stifr(i)=em20
        ENDIF
C
       ENDIF
C----
       IF(i>0)THEN
C---     output of tied contact forces
         CALL i2forces(x       ,fs     ,fx     ,fy      ,fz     ,
     .                irect(1,l),nir  ,fsav   ,fncont  ,fncontp,
     .                ftcontp ,weight ,h3d_data,i      ,h)    
       ENDIF                        
C----
      ENDDO
C
      ENDIF
C
C Traitement specifique pour ADM
C
      IF(anim_n(2)+outp_n(2)+h3d_data%N_SCAL_DMAS >0)THEN
#include "vectorize.inc"
        DO ii=1,nmn
          j=msr(ii)
          adm(j) = adm(j)/max(mmass(ii),em20)
        ENDDO
      ENDIF
C
C
      RETURN
      END