sdet43.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/.
!||====================================================================
!||    sdet43    ../engine/source/elements/solid/sconnect/sdet43.F
!||--- called by ------------------------------------------------------
!||    scoor43   ../engine/source/elements/solid/sconnect/scoor43.F
!||====================================================================
      SUBROUTINE sdet43(
     .         NEL  ,X    ,Y    ,Z    ,AREA ,AREAP,OFF  ,SOLID_ID,
     .         E1X  ,E1Y  ,E1Z  ,E2X  ,E2Y  ,E2Z  ,E3X  ,E3Y  ,E3Z  )
C-----------------------------------------------
C   I m p l i c i t   T y p e s
C-----------------------------------------------
#include      "implicit_f.inc"
CC-----------------------------------------------
C   C o m m o n   B l o c k s
C-----------------------------------------------
#include "units_c.inc"
#include "comlock.inc"
#include "mvsiz_p.inc"
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      INTEGER NEL
      INTEGER SOLID_ID(*)
C     REAL
      my_real
     .   area(nel),areap(mvsiz,4),off(nel),
     .   x(nel,8) ,y(nel,8),z(nel,8),rxx(nel),ryy(nel),
     .   e1x(*),e2x(*),e3x(*),e1y(*),e2y(*),e3y(*),e1z(*),e2z(*),e3z(*)
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER I,NINDX
      INTEGER INDX(NEL)
C     REAL
      my_real
     .   pg, xx1,xx2,xx3,xx4,yy1,yy2,yy3,yy4,zz1,zz2,zz3,zz4,
     .   l1,l2,xl1,xl2,xl3,xl4,yl2,yl3,yl4,
     .   mx13,my13,mx23,mx34,my23,my34,j0,j1,j2
      my_real 
     .   vcx(4),vcy(4),
     .   p1x(nel), p2x(nel), p3x(nel), p4x(nel),
     .   p1y(nel), p2y(nel), p3y(nel), p4y(nel),
     .   p1z(nel), p2z(nel), p3z(nel), p4z(nel)
      DATA   pg/0.5773502691896257645091488/
C=======================================================================
      nindx = 0  
C     Calcule AREA / P Gauss
      DO i=1,nel 
        p1x(i)=(x(i,1)+x(i,5))*half
        p1y(i)=(y(i,1)+y(i,5))*half
        p1z(i)=(z(i,1)+z(i,5))*half
        p2x(i)=(x(i,2)+x(i,6))*half
        p2y(i)=(y(i,2)+y(i,6))*half
        p2z(i)=(z(i,2)+z(i,6))*half
        p3x(i)=(x(i,3)+x(i,7))*half
        p3y(i)=(y(i,3)+y(i,7))*half
        p3z(i)=(z(i,3)+z(i,7))*half
        p4x(i)=(x(i,4)+x(i,8))*half
        p4y(i)=(y(i,4)+y(i,8))*half
        p4z(i)=(z(i,4)+z(i,8))*half
      ENDDO
C
      DO i=1,nel 
        xx1 = p2x(i)-p1x(i)
        yy1 = p2y(i)-p1y(i)
        zz1 = p2z(i)-p1z(i)
        xl2 = e1x(i)*xx1 + e1y(i)*yy1 + e1z(i)*zz1
        yl2 = e2x(i)*xx1 + e2y(i)*yy1 + e2z(i)*zz1
c
        xx3 = p3x(i)-p1x(i)
        yy3 = p3y(i)-p1y(i)
        zz3 = p3z(i)-p1z(i)
        xl3 = e1x(i)*xx3 + e1y(i)*yy3 + e1z(i)*zz3
        yl3 = e2x(i)*xx3 + e2y(i)*yy3 + e2z(i)*zz3
c
        xx4 = p4x(i)-p1x(i)
        yy4 = p4y(i)-p1y(i)
        zz4 = p4z(i)-p1z(i)
        xl4 = e1x(i)*xx4 + e1y(i)*yy4 + e1z(i)*zz4
        yl4 = e2x(i)*xx4 + e2y(i)*yy4 + e2z(i)*zz4
c
        l1=fourth*(xl2+xl3+xl4)
        l2=fourth*(yl2+yl3+yl4)
c
        vcx(1)=-l1
        vcx(2)= xl2-l1
        vcx(3)= xl3-l1
        vcx(4)= xl4-l1
        vcy(1)=-l2
        vcy(2)=yl2-l2
        vcy(3)=yl3-l2
        vcy(4)=yl4-l2
c
        mx13=(vcx(1)+vcx(3))*half
        my13=(vcy(1)+vcy(3))*half
        mx23=(vcx(2)+vcx(3))*half
        mx34=(vcx(3)+vcx(4))*half
        my23=(vcy(2)+vcy(3))*half
        my34=(vcy(3)+vcy(4))*half   
c
        j1 = (mx23*my13-mx13*my23)*pg
        j2 =-(mx13*my34-mx34*my13)*pg
        j0 = area(i)*fourth
C-----------
        areap(i,1) = j0+j2-j1
        areap(i,2) = j0+j2+j1
        areap(i,3) = j0-j2+j1
        areap(i,4) = j0-j2-j1
        IF (off(i) == one .and. ( area(i) < zero .or.
     .      areap(i,1) < zero .or. areap(i,2) < zero .or.
     .      areap(i,3) < zero .or. areap(i,4) < zero)) THEN
          nindx = nindx+1
          indx(nindx) = i
          off(i) = zero
        ENDIF
      ENDDO
C-----------
      IF (nindx > 0) THEN
        DO i=1,nindx
#include "lockon.inc"
          WRITE(iout ,1000) solid_id(indx(i))
          WRITE(istdo,1100) solid_id(indx(i))
#include "lockoff.inc"
        END DO
      ENDIF         
C-----------
 1000 FORMAT(5x,'NEGATIVE AREA: DELETE CONNECTION SOLID ELEMENT ',i10)
 1100 FORMAT(5x,'NEGATIVE AREA: DELETE CONNECTION SOLID ELEMENT ',i10)
C-----------
      RETURN
      END