sfor_n2s3.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/.
!||====================================================================
!||    sfor_n2s3       ../engine/source/elements/solid/solide4/sfor_n2s3.F
!||--- called by ------------------------------------------------------
!||    sfor_3n2s3      ../engine/source/elements/solid/solide/sfor_4n2s4.F90
!||    sfor_n2stria    ../engine/source/elements/solid/solide4/sfor_n2stria.F
!||    sfor_n2stria2   ../engine/source/elements/solid/solide10/sfor_n2stria2.F
!||====================================================================
      SUBROUTINE sfor_n2s3( XI,     YI,     ZI, FORC_N,     
     .                      X1,     Y1,     Z1, FOR_T1, 
     .                      X2,     Y2,     Z2, FOR_T2,
     .                      X3,     Y3,     Z3, FOR_T3,
     .                    STIF,   IFC1,  IFCTL, PENMIN,
     .                  PENREF, FKTMAX,    LL , FQMAX ,
     .                     VX1,     VX2,   VX3,  VXI  ,
     .                     VY1,     VY2,   VY3,  VYI  ,
     .                     VZ1,     VZ2,   VZ3,  VZI  ,
     .                     NEL,E_DISTOR,   DT1)
C-----------------------------------------------
C   I m p l i c i t   T y p e s
C-----------------------------------------------
#include      "implicit_f.inc"
C-----------------------------------------------
C   G l o b a l   P a r a m e t e r s
C-----------------------------------------------
#include      "mvsiz_p.inc"
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      INTEGER, INTENT (IN)  :: NEL
      INTEGER, INTENT (OUT) :: IFCTL
      INTEGER, DIMENSION(MVSIZ), INTENT (IN) :: IFC1
      my_real, DIMENSION(MVSIZ), INTENT (IN) ::  X1, X2, X3,
     .                                           Y1, Y2, Y3,
     .                                           Z1, Z2, Z3,
     .                                           XI, YI, ZI,
     .                     VX1,     VX2,   VX3,  VXI  ,
     .                     VY1,     VY2,   VY3,  VYI  ,
     .                     vz1,     vz2,   vz3,  vzi  
      my_real, INTENT (IN) :: fqmax,dt1
      my_real, DIMENSION(MVSIZ), INTENT (INOUT) :: fktmax
      my_real, DIMENSION(MVSIZ), INTENT (IN) :: stif,ll,penmin, penref
      my_real, DIMENSION(MVSIZ,3), INTENT (INOUT) :: forc_n, for_t1, for_t2, for_t3
      my_real, DIMENSION(NEL),  INTENT (INOUT) :: e_distor
C-----------------------------------------------
C   C o m m o n   B l o c k s
C-----------------------------------------------
C
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER I,J,IFCTL1,ie
C                                                                     12
      my_real nx(mvsiz),ny(mvsiz),nz(mvsiz),fn(mvsiz),
     .        xsc(mvsiz),ysc(mvsiz),zsc(mvsiz),pene(mvsiz),
     .        xa(mvsiz),ya(mvsiz),za(mvsiz),area(mvsiz),
     .        xb(mvsiz),yb(mvsiz),zb(mvsiz),
     .        xc(mvsiz),yc(mvsiz),zc(mvsiz),
     .        la(mvsiz),lb(mvsiz),lc(mvsiz),stifkt(mvsiz),
     .        rx, ry, rz, sx, sy, sz,vxc,vyc,vzc,
     .        x42,y42, z42, x31, y31, z31,fx,fy,fz,
     .        sax,say,saz,sbx,sby,sbz,scx,scy,scz,
     .        trx,try,trz,tsx,tsy,tsz,ttx,tty,ttz,
     .        tr2,ts2,tt2,aaa,bbb,vr,vs,vt,nnx,nny,nnz,
     .        xab,xbc,xca,yab,ybc,yca,zab,zbc,zca,
     .        xia,  xib,  xic, yia,  yib,  yic,
     .        zia,  zib,  zic, h0,norm,s2,fac,s2min,lj,
     .        f_q,f_c,kts,zerom,tx,ty,tz,pendr,dx,dy,dz,dn
C----------------------------
        ifctl = 0
        ifctl1 = 0
C-------diff in velocity as 1er sorting
         DO i=1,nel
           IF (ifc1(i)==0) cycle
           ifctl1=1
         END DO
C
       IF (ifctl1==1) THEN
         pene(1:nel) =zero
         zerom = -two*em03
         DO i=1,nel
           IF (ifc1(i)==0) cycle
           rx =x2(i)-x1(i)
           ry =y2(i)-y1(i)
           rz =z2(i)-z1(i)
           sx =x3(i)-x1(i)
           sy =y3(i)-y1(i)
           sz =z3(i)-z1(i)
           nx(i)=ry*sz - rz*sy
           ny(i)=rz*sx - rx*sz
           nz(i)=rx*sy - ry*sx
           area(i) = nx(i)*nx(i)+ny(i)*ny(i)+nz(i)*nz(i)
           norm=one/max(em20,sqrt(area(i)))
           nx(i)=nx(i)*norm
           ny(i)=ny(i)*norm
           nz(i)=nz(i)*norm
           bbb = (x3(i)-xi(i))*nx(i) +
     .           (y3(i)-yi(i))*ny(i) +
     .           (z3(i)-zi(i))*nz(i) -penmin(i)
           pene(i) = max(zero,-bbb)  
           IF (area(i)<penmin(i)*ll(i)) pene(i)=min(pene(i),em01*penmin(i))
         ENDDO
!     3  A 
!    / \
!   /   \
!  /     \
! 1-------2 C
         DO i=1,nel
           IF(pene(i) == zero) cycle
C-------- if other sub-tria ? still using elementary nodal normal      
            xa(i) = x3(i)   
            ya(i) = y3(i)    
            za(i) = z3(i) 
            xb(i) = x1(i)   
            yb(i) = y1(i)    
            zb(i) = z1(i)   
            xc(i) = x2(i)   
            yc(i) = y2(i)    
            zc(i) = z2(i)   
         ENDDO
         DO i=1,nel
           IF(pene(i) == zero) cycle
           xab = xb(i)-xa(i)
           yab = yb(i)-ya(i)
           zab = zb(i)-za(i)
           xbc = xc(i)-xb(i)
           ybc = yc(i)-yb(i)
           zbc = zc(i)-zb(i)
           xca = xa(i)-xc(i)
           yca = ya(i)-yc(i)
           zca = za(i)-zc(i)
           
           xia = xa(i)-xi(i)
           yia = ya(i)-yi(i)
           zia = za(i)-zi(i)
           xib = xb(i)-xi(i)
           yib = yb(i)-yi(i)
           zib = zb(i)-zi(i)
           xic = xc(i)-xi(i)
           yic = yc(i)-yi(i)
           zic = zc(i)-zi(i)
           sx = - yab*zca + zab*yca
           sy = - zab*xca + xab*zca
           sz = - xab*yca + yab*xca
           s2 = sx*sx+sy*sy+sz*sz
           sax = yib*zic - zib*yic
           say = zib*xic - xib*zic
           saz = xib*yic - yib*xic
           la(i) = (sx*sax+sy*say+sz*saz)/s2
           sbx = yic*zia - zic*yia
           sby = zic*xia - xic*zia
           sbz = xic*yia - yic*xia
           lb(i) = (sx*sbx+sy*sby+sz*sbz)/s2
           lc(i) = one - la(i) - lb(i)
           lj = min(la(i),lb(i),lc(i))
           IF (lj<zerom) pene(i)=min(pene(i),penmin(i))
           IF(la(i)<zero)THEN
             IF(lb(i)<zero)THEN
               la(i) = zero
               lb(i) = zero
               lc(i) = one
             ELSEIF(lc(i)<zero)THEN
               lc(i) = zero
               la(i) = zero
               lb(i) = one
             ELSE
               la(i) = zero
               aaa = lb(i) + lc(i)
               lb(i) = lb(i)/aaa
               lc(i) = lc(i)/aaa
             ENDIF
           ELSEIF(lb(i)<zero)THEN
             IF(lc(i)<zero)THEN
               lb(i) = zero
               lc(i) = zero
               la(i) = one
             ELSE
               lb(i) = zero
               aaa = lc(i) + la(i)
               lc(i) = lc(i)/aaa
               la(i) = la(i)/aaa
             ENDIF
           ELSEIF(lc(i)<zero)THEN
               lc(i) = zero
               aaa = la(i) + lb(i)
               la(i) = la(i)/aaa
               lb(i) = lb(i)/aaa
           ENDIF
         ENDDO
         f_q = ep02
         DO i=1,nel
           IF(pene(i) == zero) cycle
           pendr  = (pene(i)/penref(i))**2
           fac = min(fqmax,f_q*pendr)
           fktmax(i) =max(fktmax(i),(one+three*fac))
           fn(i) = (fac+one)*stif(i)*pene(i)
         ENDDO
         DO i=1,nel
           IF(pene(i) == zero) cycle
            dx = vxi(i) - lb(i)*vx1(i) - lc(i)*vx2(i)- la(i)*vx3(i)
            dy = vyi(i) - lb(i)*vy1(i) - lc(i)*vy2(i)- la(i)*vy3(i)
            dz = vzi(i) - lb(i)*vz1(i) - lc(i)*vz2(i)- la(i)*vz3(i)
            dn = (nx(i)*dx + ny(i)*dy + nz(i)*dz)*dt1
            e_distor(i) = e_distor(i) - fn(i)*dn
         ENDDO
        DO i=1,nel
          IF (pene(i) ==zero) cycle
          fx = nx(i)*fn(i)
          fy = ny(i)*fn(i)
          fz = nz(i)*fn(i)
          forc_n(i,1) = forc_n(i,1) - fx
          forc_n(i,2) = forc_n(i,2) - fy
          forc_n(i,3) = forc_n(i,3) - fz
          for_t1(i,1) = for_t1(i,1) + fx*lb(i)
          for_t1(i,2) = for_t1(i,2) + fy*lb(i)
          for_t1(i,3) = for_t1(i,3) + fz*lb(i)
          for_t2(i,1) = for_t2(i,1) + fx*lc(i)
          for_t2(i,2) = for_t2(i,2) + fy*lc(i)
          for_t2(i,3) = for_t2(i,3) + fz*lc(i)
          for_t3(i,1) = for_t3(i,1) + fx*la(i)
          for_t3(i,2) = for_t3(i,2) + fy*la(i)
          for_t3(i,3) = for_t3(i,3) + fz*la(i)
          ifctl = 1
        ENDDO
       END IF ! (IFCTL1==1) THEN
c-----------
      RETURN
      END