ingrbric_dx.F File Reference

#include "implicit_f.inc"
#include "com04_c.inc"
#include "param_c.inc"

Go to the source code of this file.

Functions/Subroutines
subroutine	ingrbric_dx (nbric, ibufssg, global_gap, ixs, x, noint, titr, is_gap_computed, pm, ipm, iddlevel, istiff, auto_rho, auto_length, multi_fvm)

Function/Subroutine Documentation

ingrbric_dx()

subroutine ingrbric_dx	(	integer, intent(in)	nbric,
		integer, dimension(*), intent(in)	ibufssg,
		intent(inout)	global_gap,
		integer, dimension(nixs,numels), intent(in)	ixs,
		dimension(3,numnod), intent(in)	x,
		integer, intent(in)	noint,
		character(len=nchartitle), intent(in)	titr,
		logical, intent(inout)	is_gap_computed,
		dimension(npropm,nummat), intent(in)	pm,
		integer, dimension(npropmi,nummat), intent(in)	ipm,
		integer, intent(in)	iddlevel,
		integer, intent(in)	istiff,
		intent(inout)	auto_rho,
		intent(inout)	auto_length,
		type(multi_fvm_struct), intent(in)	multi_fvm )

Definition at line 32 of file ingrbric_dx.F.

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
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 a i 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 , max_ratio, 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