detcord.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/.
!||====================================================================
!||    detcord                ../starter/source/initial_conditions/detonation/detcord.F
!||--- called by ------------------------------------------------------
!||    m5in2                  ../starter/source/initial_conditions/detonation/m5in2.F
!||    m5in2t                 ../starter/source/initial_conditions/detonation/m5in2t.F
!||    m5in3                  ../starter/source/initial_conditions/detonation/m5in3.F
!||--- calls      -----------------------------------------------------
!||--- uses       -----------------------------------------------------
!||    detonators_mod         ../starter/share/modules1/detonators_mod.F
!||====================================================================
      SUBROUTINE detcord(DETONATOR_CORD,X,XC,YC,ZC,VDET,VDET2,ALT,BT,TB,HAS_DETONATOR,IOPT)
C-----------------------------------------------
C   D e s c r i p t i o n
C-----------------------------------------------
      !IOPT=0  : def=3
      !IOPT=1  : piecewise linear  - multiple segments (experimental / osbslete)
      !IOPT=2  : instantaneous - multiple segments (experimental / obsolete)
      !IOPT=3  : Centripetal-Catmull-Rom SPLINE interpolation + projection along neutral fiber
C-----------------------------------------------
C   M o d u l e s
C-----------------------------------------------
      USE detonators_mod
C-----------------------------------------------
C   I m p l i c i t   T y p e s
C-----------------------------------------------
#include      "implicit_f.inc"
C-----------------------------------------------
C   C o m m o n   B l o c k s
C-----------------------------------------------
#include      "mvsiz_p.inc"
#include      "vect01_c.inc"
#include      "scr11_c.inc"
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      INTEGER :: IOPT
      my_real :: x(3,*),xc(mvsiz),yc(mvsiz),zc(mvsiz),bt(mvsiz),vdet,vdet2,alt,tb(mvsiz)
      TYPE(detonator_cord_struct_)::DETONATOR_CORD
      LOGICAL HAS_DETONATOR(MVSIZ)
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER :: J, I1,I2, II,K, iTdet,  i, NPTS, NSEG,IDIST(MVSIZ),NDETCORD,NNOD
      my_real :: ddmx, dtos, xlp2, ylp2, xlp1     , zlp2 ,ylp1,
     .           zlp1, xl0 , yl0 , zl0 , xl1      , yl1  ,zl1 , xl2 ,yl2 ,zl2,
     .           ps1 , ps2 , dl1 , dl2 , dl(mvsiz), s1   ,s2  , s3,  tdetc, tp1, tp2, th, th1,th2,
     .           tc1,tc2,tc0,tc,p1p2,dh,alpha,knots(4),zp(3),
     .           zh1(3),dist1,t1,zh2(3),dist2,t2,dist(mvsiz),xx,yy,zz,d,local_pt(4,3),c(3),t,dd,
     .           len,len1,len2
 
      type spline_path
         INTEGER :: NUM_POINTS
         my_real,ALLOCATABLE,DIMENSION(:,:) :: control_point
         my_real,ALLOCATABLE,DIMENSION(:)   :: length
         my_real,ALLOCATABLE,DIMENSION(:)   :: cumulative_length
         my_real,ALLOCATABLE,DIMENSION(:,:) :: knots
      end type spline_path
 
      type(spline_path),TARGET :: USER_SPLINE_PATH
 
      my_real, POINTER, DIMENSION(:,:) :: ptr
 
C-----------------------------------------------
C   S o u r c e   L i n e s
C-----------------------------------------------
      nnod = detonator_cord%NUMNOD
      dtos = ep20
 
      IF(iopt == 2)THEN
      !INSTANTANEOUS IGNITION OF THE WHOLE CORD
        DO j=1,nnod-1
          !first node
          ii   = detonator_cord%NODES(j)
          xlp1 = x(1,ii)
          ylp1 = x(2,ii)
          zlp1 = x(3,ii)
          !second node
          ii   = detonator_cord%NODES(j+1)
          xlp2 = x(1,ii)
          ylp2 = x(2,ii)
          zlp2 = x(3,ii)
          !
          xl0 = (xlp1-xlp2)
          yl0 = (ylp1-ylp2)
          zl0 = (zlp1-zlp2)
          !loop on elems to compute burning time from line [P1,P2]
          DO i=lft,llt
            xl1 = (xc(i)-xlp1)
            yl1 = (yc(i)-ylp1)
            zl1 = (zc(i)-zlp1)
            xl2 = (xc(i)-xlp2)
            yl2 = (yc(i)-ylp2)
            zl2 = (zc(i)-zlp2)
            ps1 = xl1*xl0+yl1*yl0+zl1*zl0
            ps2 = xl2*xl0+yl2*yl0+zl2*zl0
            IF(ps1*ps2 > zero)THEN
              !projection point not on the segment [P1,P2]
              dl1 = sqrt(xl1**2+yl1**2+zl1**2)
              dl2 = sqrt(xl2**2+yl2**2+zl2**2)
              dl(i)=min(dl1,dl2)
            ELSE
              !projection along the segment [P1,P2]
              s1 = yl1*zl0 - zl1*yl0
              s2 = zl1*xl0 - xl1*zl0
              s3 = xl1*yl0 - yl1*xl0
              dl(i)=sqrt((s1**2+s2**2+s3**2)/(xl0**2+yl0**2+zl0**2))
            ENDIF
            bt(i) =alt+dl(i)/vdet
            IF(bt(i) < abs(tb(i))) tb(i)=-bt(i)
          END DO !next I
        END DO !J=1,NNOD
 
      ELSEIF(iopt == 1)THEN
      ! DETONATING CORD HAS ITS OWN DETONATION VELOCITY ALONG THE PATH
        DO j=1,nnod-1
          !first node
          ii   = detonator_cord%NODES(j)
          xlp1 = x(1,ii)
          ylp1 = x(2,ii)
          zlp1 = x(3,ii)
          tp1  = detonator_cord%TDET_PATH(j)
          !second node
          ii   = detonator_cord%NODES(j+1)
          xlp2 = x(1,ii)
          ylp2 = x(2,ii)
          zlp2 = x(3,ii)
          tp2  = detonator_cord%TDET_PATH(j+1)
          !loop on elems to compute burning time from line [P1,P2]
          DO i=lft,llt
            tdetc = zero
            xl0 = (xlp1-xlp2)
            yl0 = (ylp1-ylp2)
            zl0 = (zlp1-zlp2)
            xl1 = (xc(i)-xlp1)
            yl1 = (yc(i)-ylp1)
            zl1 = (zc(i)-zlp1)
            xl2 = (xc(i)-xlp2)
            yl2 = (yc(i)-ylp2)
            zl2 = (zc(i)-zlp2)
            ps1 = xl1*xl0+yl1*yl0+zl1*zl0
            ps2 = xl2*xl0+yl2*yl0+zl2*zl0
            dl1 = sqrt(xl1**2+yl1**2+zl1**2)
            dl2 = sqrt(xl2**2+yl2**2+zl2**2)
            tc1 = tp1 + dl1  /vdet
            tc2 = tp2 + dl2  /vdet
            tc  = min(tc1,tc2)
            IF(ps1*ps2 <= zero)THEN
              s1    = yl1*zl0 - zl1*yl0
              s2    = zl1*xl0 - xl1*zl0
              s3    = xl1*yl0 - yl1*xl0
              !
              p1p2  = (xl0**2+yl0**2+zl0**2)
              dl2   = (s1**2+s2**2+s3**2)/p1p2
              dl(i) = sqrt(dl2)
              p1p2  = sqrt(p1p2)
              dh    = sqrt((xl1**2+yl1**2+zl1**2)-dl2)  !P1H
              !
              th1   = tp1 + dh/vdet2
              th2   = tp2 + (p1p2-dh)/vdet2
              th    = min(th1,th2)
              tc0   = th  + dl(i)/vdet
              tc    = min(tc,tc0)
            ENDIF
            bt(i)     = tc
            IF(bt(i) < abs(tb(i))) tb(i)=-bt(i)
          END DO !next I
        END DO !J=1,NNOD
 
      ELSEIF(iopt == 3)THEN
 
        !CATMULL-ROM PARAMETER
        ! ALPHA = 0.0 : uniform
        ! ALPHA = 0.5 : centripetal
        ! ALPHA = 1.0 : chordal
        alpha = half
 
        IF(vdet2 == zero)vdet2=vdet
 
        !initialize burning times by following the neutral fiber.
        npts = nnod    !number of point composing the cord path
        nseg = nnod-1  !number of segment composing the cord path
        ALLOCATE(user_spline_path%control_point(0:npts+1,3))
        ALLOCATE(user_spline_path%length(npts-1))
        ALLOCATE(user_spline_path%cumulative_length(npts-1))
        ALLOCATE(user_spline_path%knots(npts-1,4) )
        !PT(0    ,1:3)= will be defined by energy minimum
        !PT(NPTS+1,1:3)= will be defined by energy minimum
 
        !===CONTROL POINTS===!
        DO j=1,nnod
          ii   = detonator_cord%NODES(j)
          xlp1 = x(1,ii)
          ylp1 = x(2,ii)
          zlp1 = x(3,ii)
          user_spline_path%control_point(j,1:3)=(/xlp1, ylp1, zlp1/)
        ENDDO
 
        !===END POINTS===!
        ptr=>user_spline_path%control_point(0:npts+1,1:3)
        !ADDING FIRST END POINTS - BENDING ENERGY MINIMUM
        ptr(1+0,1)=half*(five*ptr(1+1,1)-four*ptr(1+2,1)+ptr(1+3,1))
        ptr(1+0,2)=half*(five*ptr(1+1,2)-four*ptr(1+2,2)+ptr(1+3,2))
        ptr(1+0,3)=half*(five*ptr(1+1,3)-four*ptr(1+2,3)+ptr(1+3,3))
        !ADDING LAST END POINTS - BENDING ENERGY MINIMUM
        ptr(1+npts+1,1)=half*(1.*ptr(1+npts-2,1)-four*ptr(1+npts-1,1)+five*ptr(1+npts,1))
        ptr(1+npts+1,2)=half*(1.*ptr(1+npts-2,2)-four*ptr(1+npts-1,2)+five*ptr(1+npts,2))
        ptr(1+npts+1,3)=half*(1.*ptr(1+npts-2,3)-four*ptr(1+npts-1,3)+five*ptr(1+npts,3))
 
        !DEFINE KNOTS
        DO j=1,nseg
          CALL cr_spline_knots(ptr(j  ,1), knots, alpha)
          user_spline_path%knots(j,1:4)=knots(1:4)
        ENDDO
 
       !COMPUTE DISTANCE FROM CENTROID Z TO CONTROL POINT ptr(1+:)
       !  idist(J) is the shortest point in PT from Z(J)
       DO j=lft,llt
         dist(j)  = 1e20
         idist(j) = 0
         DO i=1,npts
           xx = xc(j)-ptr(1+i,1)
           yy = yc(j)-ptr(1+i,2)
           zz = zc(j)-ptr(1+i,3)
           dd = sqrt(xx*xx+yy*yy+zz*zz)
           IF(dd < dist(j))THEN
             dist(j)  = dd
             idist(j) = i
           ENDIF
         ENDDO
       ENDDO
 
       !SPLINE LENGTHS
       DO i=1,nseg
         local_pt(1,1:3)=ptr(1+i-1,1:3)
         local_pt(2,1:3)=ptr(1+i+0,1:3)
         local_pt(3,1:3)=ptr(1+i+1,1:3)
         local_pt(4,1:3)=ptr(1+i+2,1:3)
         t=1.0
         CALL cr_spline_length(local_pt,alpha,t,len)
         user_spline_path%length(i) = len
         IF(i > 1)THEN
           user_spline_path%cumulative_length(i) = user_spline_path%cumulative_length(i-1)  + len
         ELSE
           user_spline_path%cumulative_length(1) =  len
         ENDIF
         !print *, "spline-i=", I
         !print *, "len =", len
       ENDDO
 
       !TEST POINT PROJECTIONS
       DO j=lft,llt
         !nearest point : IDIST(J)
         i=idist(j)
         i=min(i,nseg)
         k=i
         t=half
         local_pt(1,1:3)=ptr(1+i-1,1:3)
         local_pt(2,1:3)=ptr(1+i+0,1:3)
         local_pt(3,1:3)=ptr(1+i+1,1:3)
         local_pt(4,1:3)=ptr(1+i+2,1:3)
         zp(1:3) =(/xc(j),yc(j),zc(j)/)
         CALL cr_spline_point_proj(local_pt,zp(1:3),alpha,zh1,dist1,t1)
         CALL cr_spline_length(local_pt,alpha,t1,len1)
         c(1:3)=zh1(1:3)
         t=t1 ! position [0,1]
         len=len1
         !two adjacent spline on given node IDIST(J), unless for first spline I==1
         IF(i >= 2)THEN
           i=i-1
           local_pt(1,1:3)=ptr(1+i-1,1:3)
           local_pt(2,1:3)=ptr(1+i+0,1:3)
           local_pt(3,1:3)=ptr(1+i+1,1:3)
           local_pt(4,1:3)=ptr(1+i+2,1:3)
           CALL cr_spline_point_proj(local_pt,zp(1:3),alpha,zh2,dist2,t2)
           CALL cr_spline_length(local_pt,alpha,t2,len2)
           IF(dist2 < dist1)THEN
             c(1:3)=zh2(1:3)
             t=t2
             k=i
             len = len2
           ENDIF
         ENDIF
         !print *, "point=", ZP(1:3)
         !print *, "    ->", C(1:3)
         !!CALL CR_SPLINE_LENGTH(LOCAL_PT,ALPHA,T,LEN)
         IF(k > 1)len=len+user_spline_path%cumulative_length(k-1)
         bt(j) = detonator_cord%TDET_PATH(1) + len/vdet2
         IF(bt(j) < abs(tb(j))) tb(j)=-bt(j)
       ENDDO
 
      ENDIF
 
      DO i=lft,llt
        has_detonator(i) = .true.
      ENDDO
 
      dto=dtos
 
C-----------------------------------------------
      RETURN
      END SUBROUTINE detcord