ingrbric_dx.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/.
!||====================================================================
!||    ingrbric_dx            ../starter/source/interfaces/interf1/ingrbric_dx.f
!||--- called by ------------------------------------------------------
!||    lecint                 ../starter/source/interfaces/interf1/lecint.F
!||--- calls      -----------------------------------------------------
!||    ancmsg                 ../starter/source/output/message/message.F
!||--- uses       -----------------------------------------------------
!||    message_mod            ../starter/share/message_module/message_mod.F
!||====================================================================
      SUBROUTINE ingrbric_dx(NBRIC    , IBUFSSG, GLOBAL_GAP     , IXS        , X      ,  
     .                       NOINT    , TITR   , IS_GAP_COMPUTED, PM         , IPM    ,
     .                       IDDLEVEL , ISTIFF , AUTO_RHO       , AUTO_LENGTH,
     .                       MULTI_FVM)
C-----------------------------------------------
C   D e s c r i p t i o n
C-----------------------------------------------
C This subroutine is computing mesh size from the brick group which is
C related to interface type 18. A gap value is then set consequently.
C
C A check is also introduced about aspect ratio 
C because 
C  - in this case it is not obvious to determine if computed gap is the expected one
C  - it is recommended to use uniform mesh size with colocated scheme : ie law151
C-----------------------------------------------
C   M o d u l e s
C-----------------------------------------------
      USE message_mod
      USE multi_fvm_mod
      USE names_and_titles_mod , ONLY : nchartitle
      use element_mod , only : nixs
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      "com04_c.inc"
#include      "param_c.inc"
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      INTEGER,INTENT(IN) :: NBRIC, NOINT, IDDLEVEL,ISTIFF
      INTEGER,INTENT(IN) :: IBUFSSG(*), IXS(NIXS,NUMELS), IPM(NPROPMI,NUMMAT)
      my_real,INTENT(INOUT) :: GLOBAL_GAP
      my_real,INTENT(IN) :: x(3,numnod)
      CHARACTER(LEN=NCHARTITLE) ,INTENT(IN) :: TITR
      LOGICAL, INTENT(INOUT) :: IS_GAP_COMPUTED
      my_real , INTENT(IN) :: pm(npropm,nummat)
      my_real, INTENT(INOUT) :: auto_rho, auto_length
      TYPE(multi_fvm_struct), INTENT(IN) :: MULTI_FVM      
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER :: I, J, J2, IEDG, CONNECT1(12), CONNECT2(12),IE,IMAT, ENUM, ILAW
      my_real :: MIN_X,MIN_Y,MIN_Z
      my_real :: MAX_X,MAX_Y,MAX_Z  
      my_real :: xx,yy,zz
      my_real :: xx2,yy2,zz2
      my_real :: dx,dy,dz
      my_real :: diag, diag_max , len_edge(12), lmax, lmin, ratio2
      my_real :: rho_max, rho0
      LOGICAL :: CHECK_ASPECT 
      CHARACTER(LEN=NCHARTITLE) :: MSGTITL
      CHARACTER*10 :: CHAR_ID
C-----------------------------------------------
C   S o u r c e   L i n e s
C-----------------------------------------------
     
      is_gap_computed = .false.
     
      !-----------------------------------------
      ! COMPUTE GLOBAL GAP
      !-----------------------------------------
      IF(global_gap == zero)THEN
        !global gap required to estimate a stiffness value (automatic stiffness value when ISTIFF==2)
        !global gap also required when a constant gap is not input IGAP=0
        diag_max = em20
        DO i=1,nbric
          ie = ibufssg(i)
          IF(ie==0)EXIT
          max_x = -ep20
          max_y = -ep20
          max_z = -ep20
          min_x = ep20
          min_y = ep20
          min_z = ep20 
          DO j=2,9
            IF(ixs(j,ibufssg(i))==0)EXIT
            xx = x(1,ixs(j,ie))
            yy = x(2,ixs(j,ie)) 
            zz = x(3,ixs(j,ie)) 
            IF(xx < min_x)min_x=xx
            IF(yy < min_y)min_y=yy
            IF(zz < min_z)min_z=zz    
            IF(xx > max_x)max_x=xx
            IF(yy > max_y)max_y=yy
            IF(zz > max_z)max_z=zz 
          ENDDO
          dx = min_x-max_x
          dy = min_y-max_y
          dz = min_z-max_z                    
          diag = sqrt(dx*dx+dy*dy+dz*dz)  
          diag = sqrt(three)*diag
          diag = half*diag            !   sqrt(3)/2 * DIAG = ~ 1.5*MESH_SIZE
          IF(diag > diag_max)diag_max=diag
        END DO
        global_gap = diag_max
        auto_length = sqrt(three)*third*diag_max
        is_gap_computed = .true.
       ENDIF
       
      !-----------------------------------------
      ! DETERMINE GLOBAL DENSITY
      !-----------------------------------------
      rho_max = zero
      IF(istiff == 2)THEN
        rho_max=zero
        DO i=1,nbric
          ie=ibufssg(i)      
          IF(ie == 0)EXIT
          imat=ixs(1,ie)
          rho0=pm(89,imat)
          ilaw=ipm(2,imat)
          IF(ilaw == 51 .OR. ilaw == 151)THEN
            rho_max = max(rho_max,pm(91,imat)) ! use rho_max(1:nsubmat) in case of multi material laws
          ELSE
            rho_max = max(rho_max, rho0) ! monomaterial case
          ENDIF
        END DO
      ENDIF  
      auto_rho = rho_max     
      
      !-----------------------------------------
      ! CHECK ASPECT RATIO
      !-----------------------------------------
      check_aspect=.false.
      IF(iddlevel==1 .AND. multi_fvm%IS_USED)check_aspect=.true.
      IF (check_aspect)THEN
        !edge connectivity
        connect1(1:12)=(/1,1,1,2,2,3,3,4,5,5,6,7/)
        connect2(1:12)=(/2,4,5,3,6,4,7,8,6,8,7,8/) 
        enum=0
        DO i=1,nbric
          DO iedg=1,12
            j=1+connect1(iedg)
            j2=1+connect2(iedg)
            xx  = x(1,ixs(j,ibufssg(i)))
            yy  = x(2,ixs(j,ibufssg(i))) 
            zz  = x(3,ixs(j,ibufssg(i))) 
            xx2 = x(1,ixs(j2,ibufssg(i)))
            yy2 = x(2,ixs(j2,ibufssg(i))) 
            zz2 = x(3,ixs(j2,ibufssg(i)))
            dx = xx-xx2
            dy = yy-yy2
            dz = zz-zz2                    
            len_edge(iedg) = dx*dx + dy*dy + dz*dz 
          ENDDO
          lmin=minval(len_edge)
          lmax=maxval(len_edge)
          ratio2 = lmax/lmin ! ratio of squared values
          IF(ratio2 > 6.25 .AND. ENUM < 10)then
            char_id=' '
            WRITE(char_id,fmt='(I0)')ixs(11,ibufssg(i))
            msgtitl='CHECK ASPECT RATIO CELL ID ='//char_id 
            CALL ancmsg(msgid=1826, msgtype=msgwarning, anmode=aninfo, i1=noint, c1=titr, c2=msgtitl)
            enum=enum+1
          ENDIF
          IF(ENUM == 10)exit
        ENDDO      
      ENDIF
C------------------------------------------------------------
      RETURN
      END