iqel03.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/.
!||====================================================================
!||    iqel03   ../engine/source/ale/ale3d/iqel03.F
!||--- called by ------------------------------------------------------
!||    intal2   ../engine/source/ale/inter/intal2.F
!||====================================================================
      SUBROUTINE iqel03(X,IRECT,LMSR,MSR,NSV,ILOC,IRTL,NSN,NSEG,CRST,NOR)
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,INTENT(IN) :: NSN
      INTEGER,INTENT(IN) :: IRECT(4,*), LMSR(*), MSR(*), NSV(*), ILOC(*), IRTL(*),NSEG(*)
      my_real,INTENT(IN) :: x(3,numnod), crst(2,*)
      my_real,INTENT(INOUT) :: nor(3,*)
C-----------------------------------------------
C   C o m m o n   B l o c k s
C-----------------------------------------------
#include      "scr08_a_c.inc"
#include      "com04_c.inc"
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER II, I, J, K, L, JJ, NN, LK, NUM, NPT, K1, K2, K3, NODE1,NODE2, KK
      my_real n1, n2, n3, xx112, xx212, xx312, xx114, xx214, xx314, xx123,
     .        xx223, xx323, xx134, xx234, xx334, ss, tt, tp, tm, fs1, fs2,
     .        fs3, sp, sm, ft1, ft2, ft3, sck, tck, aj4, aj5, aj6, aj7, aj8,
     .        aj9, v1, v2, v3, xmg, xm1, xm2,
     .        xs1,ys1,zs1
C-----------------------------------------------
C   S o u r c e   L i n e s
C-----------------------------------------------
      DO ii=1,nsn
        i=nsv(ii)
        j=iloc(ii)
        k=msr(j)
        l=irtl(ii)
        DO jj=1,4
          nn = msr(irect(jj,l))
          ix(jj) = nn
          xx1(jj) = x(1,nn)
          xx2(jj) = x(2,nn)
          xx3(jj) = x(3,nn)
        ENDDO !next JJ
        xs1=x(1,i)
        ys1=x(2,i)
        zs1=x(3,i)
        xx112=xx1(1)-xx1(2)
        xx212=xx2(1)-xx2(2)
        xx312=xx3(1)-xx3(2)
        xx114=xx1(1)-xx1(4)
        xx214=xx2(1)-xx2(4)
        xx314=xx3(1)-xx3(4)
        xx123=xx1(2)-xx1(3)
        xx223=xx2(2)-xx2(3)
        xx323=xx3(2)-xx3(3)
        xx134=xx1(3)-xx1(4)
        xx234=xx2(3)-xx2(4)
        xx334=xx3(3)-xx3(4)
        ss=crst(1,ii)
        tt=crst(2,ii)
        tp=fourth*(one + tt)
        tm=fourth*(one - tt)
        fs1=tp*xx134-tm*xx112
        fs2=tp*xx234-tm*xx212
        fs3=tp*xx334-tm*xx312
        ! Specific Case
        IF(ix(3) == ix(4) .AND. tm == zero) THEN
          fs1=-xx112
          fs2=-xx212
          fs3=-xx312
        ENDIF
 
        sp=one + ss
        sm=one - ss
        ft1=-sm*xx114-sp*xx123
        ft2=-sm*xx214-sp*xx223
        ft3=-sm*xx314-sp*xx323
 
        n1=fs2*ft3-fs3*ft2
        n2=fs3*ft1-fs1*ft3
        n3=fs1*ft2-fs2*ft1
 
        sck=abs(ss) - one
        tck=abs(tt) - one
        
        lk=l
        num=nseg(j+1)-nseg(j)
        IF(num <= 4 .AND. (abs(sck) <= fiveem2.OR.abs(tck) <= fiveem2) ) THEN  
          npt = nseg(j)-1
          n1  = zero
          n2  = zero
          n3  = zero    
          IF(abs(sck) > fiveem2 .OR. abs(tck) > fiveem2) THEN
            !C----------------------------------------------------
            !C     TREATMENT SPECIFIC TO CORNERS
            !C---------------------------------------------------
            IF(k == ix(1)) THEN
              k1=1
              k2=2
              k3=4
            ELSEIF(k == ix(2)) THEN
              k1=2
              k2=3
              k3=1
            ELSEIF(k == ix(3)) THEN
              k1=3
              k2=4
              k3=2
              IF(ix(3) == ix(4)) k2=1
            ELSE
              k1=4
              k2=1
              k3=3
            ENDIF
            aj4=xx1(k2)-xx1(k1)
            aj5=xx2(k2)-xx2(k1)
            aj6=xx3(k2)-xx3(k1)
            aj7=xx1(k3)-xx1(k1)
            aj8=xx2(k3)-xx2(k1)
            aj9=xx3(k3)-xx3(k1)
            n1=n1+aj5*aj9-aj6*aj8
            n2=n2+aj6*aj7-aj4*aj9
            n3=n3+aj4*aj8-aj5*aj7
            v1=xs1-xx1(k1)
            v2=ys1-xx2(k1)
            v3=zs1-xx3(k1)
            xmg=sqrt(aj4**2+aj5**2+aj6**2)
            xm1=(v1*aj4+v2*aj5+v3*aj6)/xmg
            xmg=sqrt(aj7**2+aj8**2+aj9**2)
            xm2=(v1*aj7+v2*aj8+v3*aj9)/xmg
            node1=ix(k1)
            node2=ix(k2)
            IF(xm2 > xm1) node2=ix(k3)
            DO jj=1,num
              l=lmsr(npt+jj)
              IF(l /= lk) THEN
                DO kk=1,4
                  nn=msr(irect(kk,l))
                  ix(kk)=nn
                  xx1(kk)=x(1,nn)
                  xx2(kk)=x(2,nn)
                  xx3(kk)=x(3,nn)
                ENDDO !next KK
                IF(k == ix(1))THEN
                  k1=1
                  k2=2
                  k3=4
                ELSEIF(k == ix(2))THEN
                  k1=2
                  k2=3
                  k3=1
                ELSEIF(k == ix(3))THEN
                  k1=3
                  k2=4
                  k3=2
                  IF(ix(3) == ix(4)) k2=1
                ELSE
                  k1=4
                  k2=1
                  k3=3
                ENDIF
                IF(node2 == ix(k2) .OR. node2 == ix(k3)) THEN
                  aj4 = xx1(k2)-xx1(k1)
                  aj5 = xx2(k2)-xx2(k1)
                  aj6 = xx3(k2)-xx3(k1)
                  aj7 = xx1(k3)-xx1(k1)
                  aj8 = xx2(k3)-xx2(k1)
                  aj9 = xx3(k3)-xx3(k1)
                  n1  = n1+aj5*aj9-aj6*aj8
                  n2  = n2+aj6*aj7-aj4*aj9
                  n3  = n3+aj4*aj8-aj5*aj7
                ENDIF
              endif!IF(L /= LK)
            ENDDO !next JJ
          ELSE
            !C------------------------------------------------------
            !C     TREATMENT SPECIFIC TO CORNERS
            !C------------------------------------------------------
            DO  jj=1,num 
              l=lmsr(npt+jj)
              DO  kk=1,4
                 nn=msr(irect(kk,l))
                 ix(kk)  = nn
                 xx1(kk) = x(1,nn)
                 xx2(kk) = x(2,nn)
                 xx3(kk) = x(3,nn)
              ENDDO !next KK
              IF(k == ix(1)) THEN
                k1=1
                k2=2
                k3=4
              ELSEIF(k == ix(2))THEN
                k1=2
                k2=3
                k3=1
              ELSEIF(k == ix(3))THEN
                k1=3
                k2=4
                k3=2
                IF(ix(3) == ix(4)) k2=1
              ELSE
                k1=4
                k2=1
                k3=3
              ENDIF
              aj4=xx1(k2)-xx1(k1)
              aj5=xx2(k2)-xx2(k1)
              aj6=xx3(k2)-xx3(k1)
              aj7=xx1(k3)-xx1(k1)
              aj8=xx2(k3)-xx2(k1)
              aj9=xx3(k3)-xx3(k1)
              n1=n1+aj5*aj9-aj6*aj8
              n2=n2+aj6*aj7-aj4*aj9
              n3=n3+aj4*aj8-aj5*aj7
            ENDDO !next JJ
          ENDIF
        ENDIF
        xmg       = sqrt(n1*n1+n2*n2+n3*n3)
        n1        = n1/xmg
        n2        = n2/xmg
        n3        = n3/xmg
        nor(1,ii) = n1
        nor(2,ii) = n2
        nor(3,ii) = n3
      ENDDO !next II
C-----------------------------------------------      
      RETURN
      END