s8sderi3.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/.
!||====================================================================
!||    s8sderi3   ../engine/source/elements/solid/solide8s/s8sderi3.F
!||--- called by ------------------------------------------------------
!||    s8sforc3   ../engine/source/elements/solid/solide8s/s8sforc3.F
!||    s8ske3     ../engine/source/elements/solid/solide8s/s8ske3.F
!||====================================================================
      SUBROUTINE s8sderi3(
     1   OFFG,    OFF,     VOLDP,   NGL,
     2   KSI,     ETA,     ZETA,    WI,
     3   X1,      X2,      X3,      X4,
     4   X5,      X6,      X7,      X8,
     5   Y1,      Y2,      Y3,      Y4,
     6   Y5,      Y6,      Y7,      Y8,
     7   Z1,      Z2,      Z3,      Z4,
     8   Z5,      Z6,      Z7,      Z8,
     9   A11,     A12,     A13,     A21,
     A   A22,     A23,     A31,     A32,
     B   A33,     DN_R,    DN_S,    DN_T,
     C   INVJ,    DN_X,    DN_Y,    DN_Z,
     D   VOLN,    NEL)
C-----------------------------------------------
C   I m p l i c i t   T y p e s
C-----------------------------------------------
#include      "implicit_f.inc"
#include      "comlock.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   C o m m o n   B l o c k s
C-----------------------------------------------
#include      "units_c.inc"
#include      "scr17_c.inc"
#include      "impl1_c.inc"
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      INTEGER, INTENT(IN) :: NEL
      !INTEGER ICP,IDEG(*)
C     REAL
      my_real
     .   OFF(*),OFFG(*),KSI,ETA,ZETA,WI,VOLN(*),
     .   A11(MVSIZ), A12(MVSIZ), A13(MVSIZ), 
     .   A21(MVSIZ), A22(MVSIZ), A23(MVSIZ), 
     .   A31(MVSIZ), A32(MVSIZ), A33(MVSIZ) 
      DOUBLE PRECISION 
     .   VOLDP(*)
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER NGL(*), I, J ,ICOR,ep
C     REAL
C                                                                     12
      my_real
     .   x1(mvsiz), x2(mvsiz), x3(mvsiz), x4(mvsiz), 
     .   x5(mvsiz), x6(mvsiz), x7(mvsiz), x8(mvsiz),
     .   y1(mvsiz), y2(mvsiz), y3(mvsiz), y4(mvsiz), 
     .   y5(mvsiz), y6(mvsiz), y7(mvsiz), y8(mvsiz),
     .   z1(mvsiz), z2(mvsiz), z3(mvsiz), z4(mvsiz), 
     .   z5(mvsiz), z6(mvsiz), z7(mvsiz), z8(mvsiz),vol(mvsiz)
     
      DOUBLE PRECISION RI(8),SI(8),TI(8),DN_R(8),DN_S(8),DN_T(8),
     .   DX_R,DX_S,DX_T,DY_R,DY_S,DY_T,DZ_R,DZ_S,DZ_T,DETM1,
     .   detj(mvsiz),invj(9,mvsiz),
     .   dn_x(mvsiz,8),dn_y(mvsiz,8),dn_z(mvsiz,8),c(6,6,mvsiz),detdp
C-----------------------------------------------
      
      DATA ri /-1.0, -1.0,  1.0,  1.0,-1.0, -1.0, 1.0, 1.0/
      DATA si /-1.0, -1.0, -1.0, -1.0, 1.0,  1.0, 1.0, 1.0/
      DATA ti /-1.0,  1.0,  1.0, -1.0,-1.0,  1.0, 1.0,-1.0/
      
!sb  ksi->t   eta->r   zeta->s
      DO i=1,8
        dn_r(i) = ri(i)*(one+ksi*ti(i))*(one+zeta*si(i))/eight
        dn_s(i) = si(i)*(one+ksi*ti(i))*(one+eta*ri(i))/eight
        dn_t(i) = ti(i)*(one+eta*ri(i))*(one+zeta*si(i))/eight
      ENDDO
!     jacobien
      DO i=1,nel
        dx_r = x1(i)*dn_r(1)+x2(i)*dn_r(2)+x3(i)*dn_r(3)+x4(i)*dn_r(4)
     .        +x5(i)*dn_r(5)+x6(i)*dn_r(6)+x7(i)*dn_r(7)+x8(i)*dn_r(8)
        dx_s = x1(i)*dn_s(1)+x2(i)*dn_s(2)+x3(i)*dn_s(3)+x4(i)*dn_s(4)
     .        +x5(i)*dn_s(5)+x6(i)*dn_s(6)+x7(i)*dn_s(7)+x8(i)*dn_s(8)
        dx_t = x1(i)*dn_t(1)+x2(i)*dn_t(2)+x3(i)*dn_t(3)+x4(i)*dn_t(4)
     .        +x5(i)*dn_t(5)+x6(i)*dn_t(6)+x7(i)*dn_t(7)+x8(i)*dn_t(8)
 
        dy_r = y1(i)*dn_r(1)+y2(i)*dn_r(2)+y3(i)*dn_r(3)+y4(i)*dn_r(4)
     .        +y5(i)*dn_r(5)+y6(i)*dn_r(6)+y7(i)*dn_r(7)+y8(i)*dn_r(8)
        dy_s = y1(i)*dn_s(1)+y2(i)*dn_s(2)+y3(i)*dn_s(3)+y4(i)*dn_s(4)
     .        +y5(i)*dn_s(5)+y6(i)*dn_s(6)+y7(i)*dn_s(7)+y8(i)*dn_s(8)
        dy_t = y1(i)*dn_t(1)+y2(i)*dn_t(2)+y3(i)*dn_t(3)+y4(i)*dn_t(4)
     .        +y5(i)*dn_t(5)+y6(i)*dn_t(6)+y7(i)*dn_t(7)+y8(i)*dn_t(8)
 
        dz_r = z1(i)*dn_r(1)+z2(i)*dn_r(2)+z3(i)*dn_r(3)+z4(i)*dn_r(4)
     .        +z5(i)*dn_r(5)+z6(i)*dn_r(6)+z7(i)*dn_r(7)+z8(i)*dn_r(8)
        dz_s = z1(i)*dn_s(1)+z2(i)*dn_s(2)+z3(i)*dn_s(3)+z4(i)*dn_s(4)
     .        +z5(i)*dn_s(5)+z6(i)*dn_s(6)+z7(i)*dn_s(7)+z8(i)*dn_s(8)
        dz_t = z1(i)*dn_t(1)+z2(i)*dn_t(2)+z3(i)*dn_t(3)+z4(i)*dn_t(4)
     .        +z5(i)*dn_t(5)+z6(i)*dn_t(6)+z7(i)*dn_t(7)+z8(i)*dn_t(8)
      
        detdp  = dx_r*(dy_s*dz_t-dz_s*dy_t)
     .           -dx_s*(dy_r*dz_t-dy_t*dz_r)
     .           +dx_t*(dy_r*dz_s-dy_s*dz_r)
        detj(i) = detdp
        voldp(i) = wi*detdp
        
        IF (detj(i) > zero)THEN
          detm1 = one/detj(i)
 
          invj(1,i) = (dy_s*dz_t-dz_s*dy_t)*detm1
          invj(2,i) = (dz_r*dy_t-dy_r*dz_t)*detm1
          invj(3,i) = (dy_r*dz_s-dy_s*dz_r)*detm1
          invj(4,i) = (dx_t*dz_s-dx_s*dz_t)*detm1
          invj(5,i) = (dx_r*dz_t-dx_t*dz_r)*detm1
          invj(6,i) = (dx_s*dz_r-dx_r*dz_s)*detm1
          invj(7,i) = (dx_s*dy_t-dx_t*dy_s)*detm1
          invj(8,i) = (dx_t*dy_r-dx_r*dy_t)*detm1
          invj(9,i) = (dx_r*dy_s-dx_s*dy_r)*detm1
        ELSE
          
        ENDIF
        
        a11(i) = dx_t !R
        a12(i) = dy_t !R
        a13(i) = dz_t !R
        a21(i) = dx_r !s
        a22(i) = dy_r !S
        a23(i) = dz_r !S
        a31(i) = dx_s !T
        a32(i) = dy_s !T
        a33(i) = dz_s !T
      ENDDO
      
      DO j=1,8
        DO i=1,nel
          dn_x(i,j) = dn_r(j)*invj(1,i)+dn_s(j)*invj(2,i)+dn_t(j)*invj(3,i)
          dn_y(i,j) = dn_r(j)*invj(4,i)+dn_s(j)*invj(5,i)+dn_t(j)*invj(6,i)
          dn_z(i,j) = dn_r(j)*invj(7,i)+dn_s(j)*invj(8,i)+dn_t(j)*invj(9,i)
        ENDDO
      ENDDO
 
      icor = 0
      DO i=1,nel
        off(i)=offg(i)
        IF(off(i)==zero)THEN
         detj(i)=one
         IF (voldp(i)<=zero) voldp(i)=one
        ELSEIF (voldp(i)<=zero ) THEN
         voldp(i)= em20
         off(i) =zero
          icor=1
        ENDIF
      ENDDO
C
      IF (icor>0.AND.impl_s>0) THEN
        DO i=1,nel
          IF(voldp(i)<=zero.AND.off(i)/=zero)THEN
            voldp(i)= em20
            off(i) =zero 
           IF (imp_chk>0) THEN
#include "lockon.inc"
            WRITE(iout ,2001) ngl(i)
#include "lockoff.inc"
            idel7nok = 1
            imp_ir = imp_ir + 1
           ELSEIF (imconv==1) THEN
#include "lockon.inc"
            WRITE(istdo,2000) ngl(i)
            WRITE(iout ,2000) ngl(i)
#include "lockoff.inc"
            idel7nok = 1   
           ENDIF 
          ENDIF
        ENDDO
      END IF
      
      DO i=1,nel
       voln(i) = voldp(i)
      ENDDO
c
 
C
 2000 FORMAT(/' ZERO OR NEGATIVE SUB-VOLUME : DELETE 3D-ELEMENT NB',
     .          i10/)
 2001 FORMAT(/' ZERO OR NEGATIVE SOLID SUB-VOLUME : ELEMENT NB:',
     .          i10/)
      RETURN
      END