czcoork3.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/.
!||====================================================================
!||    czcoork3       ../engine/source/elements/shell/coquez/czcoork3.F
!||--- called by ------------------------------------------------------
!||    czke3          ../engine/source/elements/shell/coquez/czke3.F
!||--- calls      -----------------------------------------------------
!||    clskew3        ../engine/source/elements/sh3n/coquedk/cdkcoor3.F
!||    cortdir3       ../engine/source/elements/shell/coque/cortdir3.F
!||--- uses       -----------------------------------------------------
!||    elbufdef_mod   ../common_source/modules/mat_elem/elbufdef_mod.F90
!||    element_mod    ../common_source/modules/elements/element_mod.F90
!||====================================================================
        SUBROUTINE czcoork3(JFT  ,JLT  ,X     ,IXC  ,PM    ,
     1                      OFFG ,AREA ,AREA_I,VQN  ,VQ    ,
     2                      X13  ,X24  ,Y13   ,Y24  ,MX13  ,
     3                      MX23 ,MX34 ,MY13  ,MY23 ,MY34  ,
     4                      Z1   ,GEO  ,
     5                      ELBUF_STR,SMSTR,NLAY,IREP  ,NPT   ,
     6                      ISMSTR,DIR_A ,DIR_B ,
     7                      PID  ,MAT   ,NGL  ,NPLAT,IPLAT ,
     8                      CORELV,OFF  ,THK  ,NEL)
C-----------------------------------------------
C   M o d u l e s
C----------------------------------------------- 
      USE elbufdef_mod
      use element_mod , only : nixc
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8
#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"
#include      "param_c.inc"
#include      "impl1_c.inc"
#include      "comlock.inc"
#include      "units_c.inc"
#include      "scr17_c.inc"
C-----------------------------------------------
C   D U M M Y   A R G U M E N T S
C-----------------------------------------------
      INTEGER IXC(NIXC,*),JFT,JLT,IREP,NPT,ISMSTR,PID(*),MAT(*),NGL(*),NEL
      INTEGER NPLAT,NLAY,IPLAT(*)
      my_real 
     .   PM(NPROPM,*),GEO(NPROPG,*), X(3,*), 
     .   MX23(*),MY13(*),MY23(*),MY34(*),
     .   X13(*),X24(*),Y13(*),Y24(*),MX13(*),
     .   VQ(MVSIZ,3,3),AREA(*),Z1(*),MX34(*),VQN(MVSIZ,3,4),AREA_I(*)
C     .   DI(6,*),DB(3,4,*)
      my_real
     .   dir_a(nel,*),dir_b(nel,*),offg(*),off(*),
     .   corelv(mvsiz,2,4),thk(*)
      double precision
     .  smstr(*)
      TYPE(elbuf_struct_) :: ELBUF_STR
C-----------------------------------------------
C   L O C A L   V A R I A B L E S
C-----------------------------------------------
      INTEGER NNOD,I,J,K,L,M,II(6),SPLAT,JPLAT(MVSIZ),MAT_1
      PARAMETER (NNOD = 4)
      my_real 
     .   lxyz0(3),deta1(mvsiz),rx(mvsiz), ry(mvsiz), rz(mvsiz),
     .   sx(mvsiz),sy(mvsiz),r11(mvsiz),r12(mvsiz),r13(mvsiz),
     .   r21(mvsiz),r22(mvsiz),r23(mvsiz),r31(mvsiz),r32(mvsiz),
     .   r33(mvsiz),xl2(mvsiz),xl3(mvsiz),xl4(mvsiz),yl2(mvsiz),
     .   yl3(mvsiz),yl4(mvsiz),ssz(mvsiz),
     .   l13(mvsiz),l24(mvsiz),ll(mvsiz),
     .   tol,x13_2(mvsiz),y13_2(mvsiz),x24_2(mvsiz),y24_2(mvsiz),
     .   z2(mvsiz),a_4,sz,sz1,sz2,sl,s1,
     .    
     .   j0,j1,j2
      my_real 
     .   xx1,xx2,xx3,xx4,yy1,yy2,yy3,yy4,zz1,zz2,zz3,zz4
C-----------------------------------------------
      DO i=1,6
        ii(i) = nel*(i-1)
      ENDDO
C
      tol=em8
      mat_1 = ixc(1,jft)
      DO i=jft,jlt
        mat(i) = mat_1
        pid(i) = ixc(6,i)
        ngl(i) = ixc(7,i)
      ENDDO
C----------------------------
C     LOCAL SYSTEM
C----------------------------
      DO i=jft,jlt
        rx(i)=x(1,ixc(3,i))+x(1,ixc(4,i))-x(1,ixc(2,i))-x(1,ixc(5,i))
        sx(i)=x(1,ixc(4,i))+x(1,ixc(5,i))-x(1,ixc(2,i))-x(1,ixc(3,i))
        ry(i)=x(2,ixc(3,i))+x(2,ixc(4,i))-x(2,ixc(2,i))-x(2,ixc(5,i))
        sy(i)=x(2,ixc(4,i))+x(2,ixc(5,i))-x(2,ixc(2,i))-x(2,ixc(3,i))
        rz(i)=x(3,ixc(3,i))+x(3,ixc(4,i))-x(3,ixc(2,i))-x(3,ixc(5,i))
       ssz(i)=x(3,ixc(4,i))+x(3,ixc(5,i))-x(3,ixc(2,i))-x(3,ixc(3,i))
      ENDDO 
      k = 0
      CALL clskew3(jft,jlt,k,
     .   rx, ry, rz, 
     .   sx, sy, ssz, 
     .   r11,r12,r13,r21,r22,r23,r31,r32,r33,deta1,offg )
      DO i=jft,jlt
        area(i)=fourth*deta1(i)
        area_i(i)=one/area(i)
        vq(i,1,1)=r11(i)
        vq(i,2,1)=r21(i)
        vq(i,3,1)=r31(i)
        vq(i,1,2)=r12(i)
        vq(i,2,2)=r22(i)
        vq(i,3,2)=r32(i)
        vq(i,1,3)=r13(i)
        vq(i,2,3)=r23(i)
        vq(i,3,3)=r33(i)
      ENDDO 
C--------------------------
C     TRANSMET GLOBAL-->LOCAL
C--------------------------
      DO i=jft,jlt
        j=ixc(2,i)
        k=ixc(3,i)
        l=ixc(4,i)
        m=ixc(5,i)
        lxyz0(1)=fourth*(x(1,l)+x(1,m)+x(1,j)+x(1,k))
        lxyz0(2)=fourth*(x(2,l)+x(2,m)+x(2,j)+x(2,k))
        lxyz0(3)=fourth*(x(3,l)+x(3,m)+x(3,j)+x(3,k))
C
        xx1=x(1,k)-x(1,j)
        yy1=x(2,k)-x(2,j)
        zz1=x(3,k)-x(3,j)
C
        xl2(i)=r11(i)*xx1+r21(i)*yy1+r31(i)*zz1
        yl2(i)=r12(i)*xx1+r22(i)*yy1+r32(i)*zz1
C
        xx2=x(1,j)-lxyz0(1)
        yy2=x(2,j)-lxyz0(2)
        zz2=x(3,j)-lxyz0(3)
        z1(i)=r13(i)*xx2+r23(i)*yy2+r33(i)*zz2
C
        xx3=x(1,l)-x(1,j)
        yy3=x(2,l)-x(2,j)
        zz3=x(3,l)-x(3,j)
        xl3(i)=r11(i)*xx3+r21(i)*yy3+r31(i)*zz3
        yl3(i)=r12(i)*xx3+r22(i)*yy3+r32(i)*zz3
C
        xx4=x(1,m)-x(1,j)
        yy4=x(2,m)-x(2,j)
        zz4=x(3,m)-x(3,j)
        xl4(i)=r11(i)*xx4+r21(i)*yy4+r31(i)*zz4
        yl4(i)=r12(i)*xx4+r22(i)*yy4+r32(i)*zz4
      ENDDO
C----------------------------
C     SMALL STRAIN
C----------------------------
      IF(ismstr==1.OR.ismstr==2)THEN
        DO i=jft,jlt
          IF(abs(offg(i))==two)THEN
            xl2(i)=smstr(ii(1)+i)
            yl2(i)=smstr(ii(2)+i)
            xl3(i)=smstr(ii(3)+i)
            yl3(i)=smstr(ii(4)+i)
            xl4(i)=smstr(ii(5)+i)
            yl4(i)=smstr(ii(6)+i)
            z1(i)=zero
            area(i)=half*
     .              ((xl2(i)-xl4(i))*yl3(i)-xl3(i)*(yl2(i)-yl4(i)))
            area_i(i)=one/max(em20,area(i))
          ELSE
            smstr(ii(1)+i)=xl2(i)
            smstr(ii(2)+i)=yl2(i)
            smstr(ii(3)+i)=xl3(i)
            smstr(ii(4)+i)=yl3(i)
            smstr(ii(5)+i)=xl4(i)
            smstr(ii(6)+i)=yl4(i)
         ENDIF
        ENDDO
      ENDIF
      IF(ismstr==1)THEN
        DO i=jft,jlt
            IF(offg(i)==one)offg(i)=two
        ENDDO
      ENDIF
C----------------------------
C     Orthotropy (later)
C----------------------------
      IF (irep > 0) THEN
       CALL cortdir3(elbuf_str,dir_a,dir_b ,jft    ,jlt    ,
     .               nlay   ,irep   ,rx    ,ry     ,rz     , 
     .               sx     ,sy     ,ssz   ,r11    ,r21    ,
     .               r31    ,r12    ,r22   ,r32    ,nel    )
      ENDIF
 
      DO i=jft,jlt
        lxyz0(1)=fourth*(xl2(i)+xl3(i)+xl4(i))
        lxyz0(2)=fourth*(yl2(i)+yl3(i)+yl4(i))
        corelv(i,1,1)=-lxyz0(1)
        corelv(i,1,2)=xl2(i)-lxyz0(1)
        corelv(i,1,3)=xl3(i)-lxyz0(1)
        corelv(i,1,4)=xl4(i)-lxyz0(1)
        corelv(i,2,1)=-lxyz0(2)
        corelv(i,2,2)=yl2(i)-lxyz0(2)
        corelv(i,2,3)=yl3(i)-lxyz0(2)
        corelv(i,2,4)=yl4(i)-lxyz0(2)
        x13(i)=(corelv(i,1,1)-corelv(i,1,3))*half
        x24(i)=(corelv(i,1,2)-corelv(i,1,4))*half
        y13(i)=(corelv(i,2,1)-corelv(i,2,3))*half
        y24(i)=(corelv(i,2,2)-corelv(i,2,4))*half
C
        mx13(i)=(corelv(i,1,1)+corelv(i,1,3))*half
        mx23(i)=(corelv(i,1,2)+corelv(i,1,3))*half
        mx34(i)=(corelv(i,1,3)+corelv(i,1,4))*half
        my13(i)=(corelv(i,2,1)+corelv(i,2,3))*half
        my23(i)=(corelv(i,2,2)+corelv(i,2,3))*half
        my34(i)=(corelv(i,2,3)+corelv(i,2,4))*half
C
        x13_2(i) =x13(i)*x13(i)
        y13_2(i) =y13(i)*y13(i)
        x24_2(i) =x24(i)*x24(i)
        y24_2(i) =y24(i)*y24(i)
        l13(i)=x13_2(i)+y13_2(i)
        l24(i)=x24_2(i)+y24_2(i)
        ll(i)=half*(l13(i)+l24(i))
                 s1=em01*thk(i)*thk(i)
        ll(i)=max(ll(i),s1)
      ENDDO
      IF (imp_chk > 0) THEN
       s1 =.577350269189626
       DO i=jft,jlt
        j1=(mx23(i)*my13(i)-mx13(i)*my23(i))*s1
        j2=-(mx13(i)*my34(i)-mx34(i)*my13(i))*s1
        j0=area(i)*fourth
        xx1=j0+j2-j1
        xx2=j0+j2+j1
        xx3=j0-j2+j1
        xx4=j0-j2-j1
        IF(offg(i)/=zero)THEN
         IF(xx1<=zero.OR.xx2<=zero.OR.
     .      xx3<=zero.OR.xx4<=zero)THEN
#include "lockon.inc"
            WRITE(iout ,2001) ngl(i)
#include "lockoff.inc"
            idel7nok = 1
            imp_ir = imp_ir + 1
         ENDIF 
        ENDIF 
       ENDDO
      ENDIF 
C--------------------------
      DO i=jft,jlt
       z2(i)=z1(i)*z1(i)
       IF (z2(i)<ll(i)*tol.OR.npt==1) THEN
         z1(i)=zero
       ELSE
C--------------------------------------------------
C WARPING SPECIAL TREATMENT
C full projection eliminer drilling rotations and rigid rotations
C--------------------------------------------------------------------------  
         a_4=area(i)*fourth
C----------  node N ----------
         sz1=mx13(i)*y24(i)-my13(i)*x24(i)
         sz2=a_4+sz1
         sz=z2(i)*l24(i)
         sl=one/sqrt(sz+sz2*sz2)
         vqn(i,1,1)=-z1(i)*y24(i)
         vqn(i,2,1)= z1(i)*x24(i)
         vqn(i,3,1)= sz2*sl
         vqn(i,1,3)=-vqn(i,1,1)
         vqn(i,2,3)=-vqn(i,2,1)
         vqn(i,1,1)= vqn(i,1,1)*sl
         vqn(i,2,1)= vqn(i,2,1)*sl
C
         sz2=a_4-sz1
         sl=one/sqrt(sz+sz2*sz2)
         vqn(i,1,3)= vqn(i,1,3)*sl
         vqn(i,2,3)= vqn(i,2,3)*sl
         vqn(i,3,3)= sz2*sl
C
         sz1=mx13(i)*y13(i)-my13(i)*x13(i)
         sz2=a_4+sz1
         sz=z2(i)*l13(i)
         sl=one/sqrt(sz+sz2*sz2)
         vqn(i,1,2)=-z1(i)*y13(i)
         vqn(i,2,2)= z1(i)*x13(i)
         vqn(i,3,2)= sz2*sl
         vqn(i,1,4)=-vqn(i,1,2)
         vqn(i,2,4)=-vqn(i,2,2)
         vqn(i,1,2)= vqn(i,1,2)*sl
         vqn(i,2,2)= vqn(i,2,2)*sl
C
         sz2=a_4-sz1
         sl=one/sqrt(sz+sz2*sz2)
         vqn(i,1,4)= vqn(i,1,4)*sl
         vqn(i,2,4)= vqn(i,2,4)*sl
         vqn(i,3,4)= sz2*sl
       ENDIF
      ENDDO
C
      nplat=jft-1
      splat= 0
      DO i=jft,jlt
       IF (z1(i)==zero) THEN
        nplat=nplat+1
        iplat(nplat)=i
       ELSE
        splat=splat+1
        jplat(splat)=i
       ENDIF
      ENDDO 
      DO i=nplat+1,jlt
       iplat(i)=jplat(i-nplat)
      ENDDO 
C
      DO i=jft,jlt
       off(i)=offg(i)
      ENDDO
C
      RETURN
 2001 FORMAT(/' ZERO OR NEGATIVE SHELL SUB-AREA : ELEMENT NB:',i8/)
      END