From 0760af85159bc94407b8fef701beebb3c10e6278 Mon Sep 17 00:00:00 2001
From: Quentin Rodier <quentin.rodier@meteo.fr>
Date: Wed, 23 Jun 2021 11:10:19 +0200
Subject: [PATCH] T.Nagel 23/06/2021: IBM Object file reading to compute the
Level Set function for any surface
---
src/MNH/ibm_generls.f90 | 562 ++++++++++++++++++++++++++++++++++++++++
src/MNH/ibm_prep_ls.f90 | 313 ++++++++++++++++++----
2 files changed, 824 insertions(+), 51 deletions(-)
create mode 100644 src/MNH/ibm_generls.f90
diff --git a/src/MNH/ibm_generls.f90 b/src/MNH/ibm_generls.f90
new file mode 100644
index 000000000..24592d216
--- /dev/null
+++ b/src/MNH/ibm_generls.f90
@@ -0,0 +1,562 @@
+!MNH_LIC Copyright 1994-2018 CNRS, Meteo-France and Universite Paul Sabatier
+!MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
+!MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
+!MNH_LIC for details. version 1.
+!
+! #######################
+MODULE MODI_IBM_GENERLS
+ ! #######################
+ !
+ INTERFACE
+ !
+ SUBROUTINE IBM_GENERLS(PIBM_FACES,PNORM_FACES,PV1,PV2,PV3,PX_MIN,PY_MIN,PX_MAX,PY_MAX,PPHI)
+ !
+ REAL, DIMENSION(:,:,:) ,INTENT(IN) :: PIBM_FACES
+ REAL, DIMENSION(:,:) ,INTENT(IN) :: PNORM_FACES,PV1,PV2,PV3
+ REAL, DIMENSION(:,:,:,:) ,INTENT(INOUT) :: PPHI
+ REAL ,INTENT(IN) :: PX_MIN,PY_MIN,PX_MAX,PY_MAX
+ !
+ END SUBROUTINE IBM_GENERLS
+ !
+ END INTERFACE
+ !
+END MODULE MODI_IBM_GENERLS
+!
+! #####################################
+SUBROUTINE IBM_GENERLS(PIBM_FACES,PNORM_FACES,PV1,PV2,PV3,PX_MIN,PY_MIN,PX_MAX,PY_MAX,PPHI)
+ ! #####################################
+ !
+ !
+ !**** IBM_GENERLS computes the Level Set function for any surface
+ !
+ ! PURPOSE
+ ! -------
+ !**** The purpose of this routine is to estimate the level set
+ ! containing XYZ minimalisation interface locations
+
+ ! METHOD
+ ! ------
+ !**** Iterative system and minimization of the interface distance
+ !
+ ! EXTERNAL
+ ! --------
+ ! SUBROUTINE ?
+ !
+ ! IMPLICIT ARGUMENTS
+ ! ------------------
+ ! MODD_?
+ !
+ ! REFERENCE
+ ! ---------
+ ! The method is based on '3D Distance from a Point to a Triangle'
+ ! a technical report from Mark W. Jones, University of Wales Swansea
+ !
+ ! AUTHORS
+ ! ------
+ ! Tim Nagel, Valéry Masson & Robert Schoetter
+ !
+ ! MODIFICATIONS
+ ! -------------
+ ! Original 01/06/2021
+ !
+ !------------------------------------------------------------------------------
+ !
+ !**** 0. DECLARATIONS
+ ! ---------------
+ !
+ ! module
+ USE MODE_POS
+ USE MODE_ll
+ USE MODE_IO
+ USE MODD_ARGSLIST_ll, ONLY : LIST_ll
+ !
+ ! declaration
+ USE MODD_IBM_PARAM_n
+ USE MODD_IBM_LSF
+ USE MODD_DIM_n, ONLY: NIMAX,NJMAX,NKMAX
+ USE MODD_PARAMETERS, ONLY: JPVEXT,JPHEXT,XUNDEF
+ USE MODD_GRID_n, ONLY: XXHAT,XYHAT,XZHAT,XZZ
+ USE MODD_METRICS_n, ONLY: XDXX,XDYY,XDZZ
+ USE MODD_VAR_ll, ONLY: IP
+ USE MODD_CST, ONLY: XMNH_EPSILON
+ !
+ ! interface
+ USE MODI_SHUMAN
+ USE MODI_IBM_INTERPOS
+ USE MODI_IBM_DETECT
+ USE MODI_INI_CST
+ !
+ IMPLICIT NONE
+ !
+ ! 0.1 declarations of arguments
+ !
+ REAL, DIMENSION(:,:,:) ,INTENT(IN) :: PIBM_FACES !faces coordinates
+ REAL, DIMENSION(:,:) ,INTENT(IN) :: PNORM_FACES !normal
+ REAL, DIMENSION(:,:) ,INTENT(IN) :: PV1,PV2,PV3
+ REAL, DIMENSION(:,:,:,:) ,INTENT(INOUT) :: PPHI ! LS functions
+ REAL ,INTENT(IN) :: PX_MIN,PY_MIN,PX_MAX,PY_MAX
+ !
+ !------------------------------------------------------------------------------
+ !
+ ! 0.2 declaration of local variables
+ !
+ INTEGER :: JI,JJ,JK,JN,JM,JI2,JJ2,JK2 ! loop index
+ INTEGER :: JI_MIN,JI_MAX,JJ_MIN,JJ_MAX,JK_MIN,JK_MAX,IIU,IJU,IKU ! loop boundaries
+ REAL :: Z_DIST_TEST1,Z_DIST_TEST2 ! saving distances
+ REAL :: Z_DIST_TEST3,Z_DIST_TEST4,ZDIST_REF0
+ INTEGER :: INUMB_FACES ! number of faces
+ REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZXHATM,ZYHATM,ZZHATM,ZDP0PP0PAST
+ CHARACTER(LEN=1) :: YPOS
+ REAL, DIMENSION(3) :: ZP1P0,ZP1P2,ZP0PP0,ZP1PP0,ZP2PP0,ZP3PP0,ZPP0P1,ZPP0P2,ZPP0P3
+ REAL, DIMENSION(3) :: ZPP0PPP0,ZPPP0P1,ZPPP0P2,ZP2P1,ZP2P0,ZP2P3,ZP3P2,ZP3P1
+ REAL, DIMENSION(3) :: ZPP0,ZFT1,ZFT2,ZFT3,ZFT1B,ZFT2B,ZFT3B,ZR,ZPPP0,ZP3P0,ZP0P1
+ REAL, DIMENSION(3) :: ZPPP0P3,ZP1P3,ZPCP0,ZR0
+ REAL, DIMENSION(:), ALLOCATABLE :: ZSTEMP,ZRDIR,ZVECTDISTPLUS,ZVECTDISTMOINS,ZVECTDIST!,ZFACE
+ REAL, DIMENSION(:,:), ALLOCATABLE :: ZC
+ REAL :: ZF1,ZF2,ZF3,ZF1B,ZF2B,ZF3B,ZDPP0PPP0
+ REAL :: ZT,ZSIGN,ZS,ZDIST,ZDP0PP0,ZNNORM,ZRN,ZPHI_OLD
+ TYPE(LIST_ll), POINTER :: TZFIELDS_ll ! list of fields to exchange
+ INTEGER :: IINFO_ll,IMI ! return code of parallel routine
+ INTEGER :: IIE,IIB,IJB,IJE,IKE,IKB,ZBPLUS
+ LOGICAL :: GABOVE_ROOF,LFACE,LDZ
+ LOGICAL, DIMENSION(:), ALLOCATABLE :: ZFACE
+ INTEGER :: ZCOUNT,ZIDX,ZII,ZCHANGE,ZCHANGE1
+ REAL :: ZDIFF,ZMIN_DIFF,ZDX
+ !
+ !------------------------------------------------------------------------------
+ !
+ ! 0.3 allocation
+ !
+ NULLIFY(TZFIELDS_ll)
+ IIU = SIZE(PPHI,1)
+ IJU = SIZE(PPHI,2)
+ IKU = SIZE(PPHI,3)
+ IIB=1+JPHEXT
+ IIE=IIU-JPHEXT
+ IJB=1+JPHEXT
+ IJE=IJU-JPHEXT
+ IKB=1+JPVEXT
+ IKE=IKU-JPVEXT
+ !
+ JK_MIN = 1 + JPVEXT
+ JK_MAX = IKU - JPVEXT
+ !
+ CALL GET_INDICE_ll (JI_MIN,JJ_MIN,JI_MAX,JJ_MAX)
+ !
+ ALLOCATE(ZXHATM(IIU,IJU,IKU))
+ ALLOCATE(ZYHATM(IIU,IJU,IKU))
+ ALLOCATE(ZZHATM(IIU,IJU,IKU))
+ !
+ !-------------------------------------------------------------------------------
+ !
+ !**** 1. PRELIMINARIES
+ ! ----------------
+ !
+ INUMB_FACES = SIZE(PIBM_FACES,1)
+ ALLOCATE(ZC(INUMB_FACES,3))
+ ALLOCATE(ZSTEMP(1))
+ ALLOCATE(ZRDIR(1))
+ PPHI = -XUNDEF
+ ALLOCATE(ZDP0PP0PAST(IIU,IJU,IKU))
+ ZDP0PP0PAST = 0.
+ ALLOCATE(ZVECTDIST(10000))
+ ALLOCATE(ZVECTDISTPLUS(10000))
+ ALLOCATE(ZVECTDISTMOINS(10000))
+ ALLOCATE(ZFACE(10000))
+ ZFACE=.FALSE.
+ !
+ !-------------------------------------------------------------------------------
+ !
+ !**** 2. EXECUTIONS
+ ! -------------
+ !
+ JM=1
+ YPOS = 'P'
+ !
+ CALL IBM_INTERPOS(ZXHATM,ZYHATM,ZZHATM,YPOS)
+ ZDX = ZXHATM(JI_MIN+1,JJ_MIN,JK_MIN)-ZXHATM(JI_MIN,JJ_MIN,JK_MIN)
+ !
+ DO JK = JK_MIN,JK_MAX
+ DO JJ = JJ_MIN,JJ_MAX
+ DO JI = JI_MIN,JI_MAX
+ ZCOUNT = 1
+ ZVECTDIST = -999.
+ DO JN = 1,INUMB_FACES
+ LFACE=.FALSE.
+ !***Calcul of the face center
+ ZC(JN,1)=(PIBM_FACES(JN,1,1)+PIBM_FACES(JN,2,1)+PIBM_FACES(JN,3,1))/3.
+ ZC(JN,2)=(PIBM_FACES(JN,1,2)+PIBM_FACES(JN,2,2)+PIBM_FACES(JN,3,2))/3.
+ ZC(JN,3)=(PIBM_FACES(JN,1,3)+PIBM_FACES(JN,2,3)+PIBM_FACES(JN,3,3))/3.
+ !***Norm normalization
+ ZNNORM = SQRT(PNORM_FACES(JN,1)**2+PNORM_FACES(JN,2)**2+PNORM_FACES(JN,3)**2)
+ !***Vector between the face center and the current grid point
+ ZPCP0(1) = ZXHATM(JI,JJ,JK)-ZC(JN,1)
+ ZPCP0(2) = ZYHATM(JI,JJ,JK)-ZC(JN,2)
+ ZPCP0(3) = ZZHATM(JI,JJ,JK)-ZC(JN,3)
+ ZSIGN = ZPCP0(1)*PNORM_FACES(JN,1)+ &
+ ZPCP0(2)*PNORM_FACES(JN,2)+ &
+ ZPCP0(3)*PNORM_FACES(JN,3)
+ !***Various vectors
+ ZP1P0(1) = ZXHATM(JI,JJ,JK)-PIBM_FACES(JN,1,1)
+ ZP1P0(2) = ZYHATM(JI,JJ,JK)-PIBM_FACES(JN,1,2)
+ ZP1P0(3) = ZZHATM(JI,JJ,JK)-PIBM_FACES(JN,1,3)
+ ZP3P0(1) = ZXHATM(JI,JJ,JK)-PIBM_FACES(JN,3,1)
+ ZP3P0(2) = ZYHATM(JI,JJ,JK)-PIBM_FACES(JN,3,2)
+ ZP3P0(3) = ZZHATM(JI,JJ,JK)-PIBM_FACES(JN,3,3)
+ ZP0P1(1) = PIBM_FACES(JN,1,1)-ZXHATM(JI,JJ,JK)
+ ZP0P1(2) = PIBM_FACES(JN,1,2)-ZYHATM(JI,JJ,JK)
+ ZP0P1(3) = PIBM_FACES(JN,1,3)-ZZHATM(JI,JJ,JK)
+ ZP2P0(1) = ZXHATM(JI,JJ,JK)-PIBM_FACES(JN,2,1)
+ ZP2P0(2) = ZYHATM(JI,JJ,JK)-PIBM_FACES(JN,2,2)
+ ZP2P0(3) = ZZHATM(JI,JJ,JK)-PIBM_FACES(JN,2,3)
+ !***Equation (3) of Jones (1995)
+ IF(ZP1P0(1)==0.AND.ZP1P0(2)==0.AND.ZP1P0(3)==0) THEN
+ WRITE(*,*) 'ZP1P0(1,2,3)',ZP1P0(1),ZP1P0(2),ZP1P0(3)
+ ZDP0PP0 = 0.
+ ELSE
+ ZDP0PP0 = SQRT(ZP0P1(1)**2+ZP0P1(2)**2+ZP0P1(3)**2)* &
+ ((ZP1P0(1)*PNORM_FACES(JN,1)+ZP1P0(2)*PNORM_FACES(JN,2)+&
+ ZP1P0(3)*PNORM_FACES(JN,3))/( &
+ SQRT((ZP1P0(1))**2+(ZP1P0(2))**2+(ZP1P0(3))**2)*ZNNORM))
+ END IF
+ !***Equation (4) of Jones (1995)
+ ZP0PP0(1) = -ZDP0PP0*(PNORM_FACES(JN,1)/ZNNORM)
+ ZP0PP0(2) = -ZDP0PP0*(PNORM_FACES(JN,2)/ZNNORM)
+ ZP0PP0(3) = -ZDP0PP0*(PNORM_FACES(JN,3)/ZNNORM)
+ !***Equation (5) of Jones (1995)
+ ZPP0(1) = ZXHATM(JI,JJ,JK)+ZP0PP0(1)
+ ZPP0(2) = ZYHATM(JI,JJ,JK)+ZP0PP0(2)
+ ZPP0(3) = ZZHATM(JI,JJ,JK)+ZP0PP0(3)
+ !
+ ZP1PP0(1)=ZPP0(1)-PIBM_FACES(JN,1,1)
+ ZP1PP0(2)=ZPP0(2)-PIBM_FACES(JN,1,2)
+ ZP1PP0(3)=ZPP0(3)-PIBM_FACES(JN,1,3)
+ !
+ ZP2PP0(1)=ZPP0(1)-PIBM_FACES(JN,2,1)
+ ZP2PP0(2)=ZPP0(2)-PIBM_FACES(JN,2,2)
+ ZP2PP0(3)=ZPP0(3)-PIBM_FACES(JN,2,3)
+ !
+ ZP3PP0(1)=ZPP0(1)-PIBM_FACES(JN,3,1)
+ ZP3PP0(2)=ZPP0(2)-PIBM_FACES(JN,3,2)
+ ZP3PP0(3)=ZPP0(3)-PIBM_FACES(JN,3,3)
+ !
+ ZPP0P1(1)=PIBM_FACES(JN,1,1)-ZPP0(1)
+ ZPP0P1(2)=PIBM_FACES(JN,1,2)-ZPP0(2)
+ ZPP0P1(3)=PIBM_FACES(JN,1,3)-ZPP0(3)
+ !
+ ZPP0P2(1)=PIBM_FACES(JN,2,1)-ZPP0(1)
+ ZPP0P2(2)=PIBM_FACES(JN,2,2)-ZPP0(2)
+ ZPP0P2(3)=PIBM_FACES(JN,2,3)-ZPP0(3)
+ !
+ ZPP0P3(1)=PIBM_FACES(JN,3,1)-ZPP0(1)
+ ZPP0P3(2)=PIBM_FACES(JN,3,2)-ZPP0(2)
+ ZPP0P3(3)=PIBM_FACES(JN,3,3)-ZPP0(3)
+ !
+ !***Calculation of f1,f2,f3 (Jones (1995))
+ ZFT1= CROSSPRODUCT(PV1(JN,:),ZP1PP0)
+ ZFT2= CROSSPRODUCT(PV2(JN,:),ZP2PP0)
+ ZFT3= CROSSPRODUCT(PV3(JN,:),ZP3PP0)
+
+ ZF1 =ZFT1(1)*PNORM_FACES(JN,1)+ &
+ ZFT1(2)*PNORM_FACES(JN,2)+ &
+ ZFT1(3)*PNORM_FACES(JN,3)
+
+ ZF2 =ZFT2(1)*PNORM_FACES(JN,1)+ &
+ ZFT2(2)*PNORM_FACES(JN,2)+ &
+ ZFT2(3)*PNORM_FACES(JN,3)
+
+ ZF3 =ZFT3(1)*PNORM_FACES(JN,1)+ &
+ ZFT3(2)*PNORM_FACES(JN,2)+ &
+ ZFT3(3)*PNORM_FACES(JN,3)
+ !***Point anticlockwise of V1 and clockwise of V2
+ IF (ZF1.GE.0.AND.ZF2.LE.0) THEN
+ ZFT1B = CROSSPRODUCT(ZPP0P1,ZPP0P2)
+ ZF1B = ZFT1B(1)*PNORM_FACES(JN,1)+ &
+ ZFT1B(2)*PNORM_FACES(JN,2)+ &
+ ZFT1B(3)*PNORM_FACES(JN,3)
+ IF (ZF1B<0) THEN
+ ZP1P2(:) = PIBM_FACES(JN,2,:)-PIBM_FACES(JN,1,:)
+ ZR = CROSSPRODUCT(CROSSPRODUCT(ZPP0P2,ZPP0P1),ZP1P2)
+ ZRN = SQRT(ZR(1)**2+ZR(2)**2+ZR(3)**2)
+ !***Eq. (10) of Jones(1995)
+ ZDPP0PPP0 = SQRT(ZPP0P1(1)**2+ZPP0P1(2)**2+ZPP0P1(3)**2)* &
+ ((ZPP0P1(1)*ZR(1)+ZPP0P1(2)*ZR(2)+ZPP0P1(3)*ZR(3))/( &
+ SQRT(ZPP0P1(1)**2+ZPP0P1(2)**2+ZPP0P1(3)**2)*ZRN))! &
+ ZPP0PPP0 = ZDPP0PPP0*(ZR/ZRN)
+ ZPPP0 = ZPP0+ZPP0PPP0
+ ZPPP0P1 = PIBM_FACES(JN,1,:)-ZPPP0
+ ZP2P1 = PIBM_FACES(JN,1,:)-PIBM_FACES(JN,2,:)
+ ZRDIR = SIGN(1.,SCALPRODUCT(ZPPP0P1,ZP2P1))
+ ZT = SQRT(ZPPP0P1(1)**2+ZPPP0P1(2)**2+ZPPP0P1(3)**2)/ &
+ SQRT(ZP2P1(1)**2+ZP2P1(2)**2+ZP2P1(3)**2)*ZRDIR(1)
+ IF (ZT.GE.0.AND.ZT.LE.1) THEN
+ ZDIST =SQRT(ZDPP0PPP0**2+ZDP0PP0**2)
+ ELSEIF (ZT<0.) THEN
+ ZDIST = SQRT(ZP1P0(1)**2+ZP1P0(2)**2+ZP1P0(3)**2)
+ ELSEIF (ZT>1.) THEN
+ ZDIST = SQRT(ZP2P0(1)**2+ZP2P0(2)**2+ZP2P0(3)**2)
+ ELSE
+ WRITE(*,*) '*****************************'
+ WRITE(*,*) '********ERROR in ZT **********'
+ WRITE(*,*) '******* calculation **********'
+ WRITE(*,*) '**** (stopped execution) ****'
+ WRITE(*,*) '*****************************'
+ call Print_msg( NVERB_FATAL, 'GEN', 'IBM_PREP_LS', 'Error in ZT calculation' )
+ ENDIF
+ ELSE
+ ZDIST = ZDP0PP0
+ LFACE = .TRUE.
+ ENDIF
+ !***Point anticlockwise of V2 and clockwise of V3
+ ELSEIF (ZF2.GE.0.AND.ZF3.LE.0) THEN
+ ZFT2B = CROSSPRODUCT(ZPP0P2,ZPP0P3)
+ ZF2B = ZFT2B(1)*PNORM_FACES(JN,1)+ &
+ ZFT2B(2)*PNORM_FACES(JN,2)+ &
+ ZFT2B(3)*PNORM_FACES(JN,3)
+ IF (ZF2B<0) THEN
+ ZP2P3(:) = PIBM_FACES(JN,3,:)-PIBM_FACES(JN,2,:)
+ ZR = CROSSPRODUCT(CROSSPRODUCT(ZPP0P3,ZPP0P2),ZP2P3)
+ ZRN = SQRT(ZR(1)**2+ZR(2)**2+ZR(3)**2)
+ ZDPP0PPP0 = SQRT(ZPP0P2(1)**2+ZPP0P2(2)**2+ZPP0P2(3)**2)* &
+ ((ZPP0P2(1)*ZR(1)+ZPP0P2(2)*ZR(2)+ZPP0P2(3)*ZR(3))/( &
+ SQRT(ZPP0P2(1)**2+ZPP0P2(2)**2+ZPP0P2(3)**2)*ZRN))! &
+ ZPP0PPP0 = ZDPP0PPP0*(ZR/ZRN)
+ ZPPP0 = ZPP0+ZPP0PPP0
+ ZPPP0P2 = PIBM_FACES(JN,2,:)-ZPPP0
+ ZP3P2 = PIBM_FACES(JN,2,:)-PIBM_FACES(JN,3,:)
+ ZRDIR = SIGN(1.,SCALPRODUCT(ZPPP0P2,ZP3P2))
+ ZT = SQRT(ZPPP0P2(1)**2+ZPPP0P2(2)**2+ZPPP0P2(3)**2)/ &
+ SQRT(ZP3P2(1)**2+ZP3P2(2)**2+ZP3P2(3)**2)*ZRDIR(1)
+ IF (ZT.GE.0.AND.ZT.LE.1) THEN
+ ZDIST = SQRT(ZDPP0PPP0**2+ZDP0PP0**2)
+ ELSEIF (ZT<0.) THEN
+ ZDIST = SQRT(ZP2P0(1)**2+ZP2P0(2)**2+ZP2P0(3)**2)
+ ELSEIF (ZT>1.) THEN
+ ZDIST = SQRT(ZP3P0(1)**2+ZP3P0(2)**2+ZP3P0(3)**2)
+ ELSE
+ WRITE(*,*) '*****************************'
+ WRITE(*,*) '********ERROR in ZT **********'
+ WRITE(*,*) '******* calculation **********'
+ WRITE(*,*) '**** (stopped execution) ****'
+ WRITE(*,*) '*****************************'
+ call Print_msg( NVERB_FATAL, 'GEN', 'IBM_PREP_LS', 'Error in ZT calculation' )
+ ENDIF
+ ELSE
+ ZDIST = ZDP0PP0
+ LFACE = .TRUE.
+ ENDIF
+ !***Point anticlockwise of V3 and clockwise of V1
+ ELSEIF (ZF3.GE.0.AND.ZF1.LE.0) THEN
+ ZFT3B = CROSSPRODUCT(ZPP0P3,ZPP0P1)
+ ZF3B = ZFT3B(1)*PNORM_FACES(JN,1)+ &
+ ZFT3B(2)*PNORM_FACES(JN,2)+ &
+ ZFT3B(3)*PNORM_FACES(JN,3)
+ IF (ZF3B<0) THEN
+ ZP3P1(:) = PIBM_FACES(JN,1,:)-PIBM_FACES(JN,3,:)
+ ZR = CROSSPRODUCT(CROSSPRODUCT(ZPP0P1,ZPP0P3),ZP3P1)
+ ZRN = SQRT(ZR(1)**2+ZR(2)**2+ZR(3)**2)
+ ZDPP0PPP0 = SQRT(ZPP0P3(1)**2+ZPP0P3(2)**2+ZPP0P3(3)**2)* &
+ ((ZPP0P3(1)*ZR(1)+ZPP0P3(2)*ZR(2)+ZPP0P3(3)*ZR(3))/( &
+ SQRT((ZPP0P3(1))**2+(ZPP0P3(2))**2+(ZPP0P3(3))**2)*ZRN))! &
+ ZPP0PPP0 = ZDPP0PPP0*(ZR/ZRN)
+ ZPPP0 = ZPP0+ZPP0PPP0
+ ZPPP0P3 = PIBM_FACES(JN,3,:)-ZPPP0
+ ZP1P3 = PIBM_FACES(JN,3,:)-PIBM_FACES(JN,1,:)
+ ZRDIR = SIGN(1.,SCALPRODUCT(ZPPP0P3,ZP1P3))
+ ZT = SQRT(ZPPP0P3(1)**2+ZPPP0P3(2)**2+ZPPP0P3(3)**2)/ &
+ SQRT(ZP1P3(1)**2+ZP1P3(2)**2+ZP1P3(3)**2)*ZRDIR(1)
+ IF (ZT.GE.0.AND.ZT.LE.1) THEN
+ ZDIST = SQRT(ZDPP0PPP0**2+ZDP0PP0**2)
+ ELSEIF (ZT<0.) THEN
+ ZDIST = SQRT(ZP3P0(1)**2+ZP3P0(2)**2+ZP3P0(3)**2)
+ ELSEIF (ZT>1.) THEN
+ ZDIST = SQRT(ZP1P0(1)**2+ZP1P0(2)**2+ZP1P0(3)**2)
+ ELSE
+ WRITE(*,*) '*****************************'
+ WRITE(*,*) '********ERROR in ZT **********'
+ WRITE(*,*) '******* calculation **********'
+ WRITE(*,*) '**** (stopped execution) ****'
+ WRITE(*,*) '*****************************'
+ call Print_msg( NVERB_FATAL, 'GEN', 'IBM_PREP_LS', 'Error in ZT calculation' )
+ ENDIF
+ ELSE
+ ZDIST = ZDP0PP0
+ LFACE = .TRUE.
+ ENDIF
+ ELSE
+ WRITE(*,*) '*****************************'
+ WRITE(*,*) '********ERROR in ZF **********'
+ WRITE(*,*) '******* instruction **********'
+ WRITE(*,*) '**** (stopped execution) ****'
+ WRITE(*,*) '*****************************'
+ call Print_msg( NVERB_FATAL, 'GEN', 'IBM_PREP_LS', 'Error in ZF instruction' )
+ ENDIF
+ ZDIST = SIGN(ZDIST,-ZSIGN)
+ ZDIST = ANINT(ZDIST*10.E5) / 10.E5
+ PPHI(JI,JJ,JK,JM) = ANINT(PPHI(JI,JJ,JK,JM)*10.E5) / 10.E5
+ IF (ABS(ZDIST).LE.ABS(PPHI(JI,JJ,JK,JM))) THEN
+ ZPHI_OLD = PPHI(JI,JJ,JK,JM)
+ IF (ABS(ZDIST)==ABS(PPHI(JI,JJ,JK,JM))) THEN
+ IF (ABS(ZDP0PP0).GT.ABS(ZDP0PP0PAST(JI,JJ,JK))) THEN
+ PPHI(JI,JJ,JK,JM) = ZDIST
+ ZDP0PP0PAST(JI,JJ,JK) = ZDP0PP0
+ ENDIF
+ ELSE
+ PPHI(JI,JJ,JK,JM) = ZDIST
+ ENDIF
+ IF (ABS(ZDIST).LT.ABS(ZPHI_OLD)) THEN
+ ZDP0PP0PAST(JI,JJ,JK) = ZDP0PP0
+ ENDIF
+ ENDIF
+ IF (ABS(PPHI(JI,JJ,JK,JM)).GT.(SQRT(3.)*4.)) THEN
+ PPHI(JI,JJ,JK,JM) = -999.
+ ENDIF
+ IF (ABS(ZDIST).LT.(SQRT(3.)*4.)) THEN
+ ZVECTDIST(ZCOUNT)=ZDIST
+ ZFACE(ZCOUNT)=LFACE
+ ZCOUNT = ZCOUNT +1
+ ENDIF
+ ENDDO
+ ZVECTDISTPLUS=ZVECTDIST
+ ZVECTDISTMOINS=ZVECTDIST
+ WHERE (ZVECTDIST.GT.0)
+ ZVECTDISTMOINS=-999.
+ ENDWHERE
+ WHERE (ZVECTDIST.LT.0)
+ ZVECTDISTPLUS=999.
+ ENDWHERE
+ IF (ANY(ZVECTDIST.GT.0.).AND.(ABS(ABS(MINVAL(ZVECTDISTPLUS))-ABS(MAXVAL(ZVECTDISTMOINS))).LT.10.E-6)) THEN
+ ZMIN_DIFF = 1.
+ ZIDX = 0
+ DO ZII = 1, SIZE(ZVECTDIST)
+ ZDIFF = ABS(ZVECTDIST(ZII)-MINVAL(ZVECTDISTPLUS))
+ IF ( ZDIFF < ZMIN_DIFF) THEN
+ ZIDX = ZII
+ ZMIN_DIFF = ZDIFF
+ ENDIF
+ ENDDO
+ IF (ZFACE(ZIDX)) THEN
+ PPHI(JI,JJ,JK,JM) = MINVAL(ZVECTDISTPLUS)
+ ENDIF
+ ENDIF
+ ENDDO
+ ENDDO
+ ENDDO
+
+DO JJ=JJ_MIN,JJ_MAX
+DO JI=JI_MIN,JI_MAX
+GABOVE_ROOF=.FALSE.
+DO JK=IKB, IKE
+ ! check if point is flagged as not calculated
+ IF (PPHI(JI,JJ,JK,JM)==-999.) THEN
+ ! check if point is already above a point that encountered a point near the
+ ! surface (that can be outside or inside a building)
+ ! check if that point was inside (if outside, the value of the levelset
+ ! stays at -999.)
+ IF (GABOVE_ROOF .AND. PPHI(JI,JJ,JK-1,JM) > XIBM_EPSI) THEN
+ PPHI(JI,JJ,JK,JM) = 999.
+ CYCLE
+ END IF
+ ! check if the point of the column have not encoutered a near-building
+ ! surface point with a physical value of the level set
+ IF (.NOT. GABOVE_ROOF) THEN
+ ! if the point above has a physical value for the level set, then the
+ ! status inside (999) or outside (-999) is given to all points below,
+ ! depending if this point above (that needs not to be the point at the top
+ ! of the model!) is inside or outside
+ ! checks if the point above has a physical value for the levelset
+ IF (JK<IKE .AND. ABS (PPHI(JI,JJ,JK+1,JM)) < 900.) THEN
+ ! if the point above is inside, all points below are set inside
+ IF (PPHI(JI,JJ,JK+1,JM)>XIBM_EPSI) PPHI(JI,JJ,IKB:JK,JM) = 999.
+ ! indicate for further processing of points above the current point
+ ! that we have encountered a physical value of the level set, near the
+ ! surface building
+ GABOVE_ROOF = .TRUE.
+ END IF
+ CYCLE
+ ENDIF
+ END IF
+ ! if we have never encoutered a roof or point near a building form above,
+ ! then, we are outside, and nothing is changed (value -999 kept)
+ END DO
+ PPHI(JI,JJ,IKB-1,JM) = PPHI(JI,JJ,IKB,JM)
+ PPHI(JI,JJ,IKE+1,JM) = PPHI(JI,JJ,IKE,JM)
+END DO
+END DO
+
+
+JN=1
+PPHI(:,:,IKB-1,JN)=2*PPHI(:,:,IKB,JN)-PPHI(:,:,IKB+1,JN)
+PPHI(:,:,IKE+1,JN)=2*PPHI(:,:,IKE,JN)-PPHI(:,:,IKE-1,JN)
+PPHI(IIB-1,:,:,JN) = PPHI( IIB ,:,:,JN)
+PPHI(IIE+1,:,:,JN) = PPHI( IIE ,:,:,JN)
+PPHI(:,IJB-1,:,JN) = PPHI(:, IJB ,:,JN)
+PPHI(:,IJE+1,:,JN) = PPHI(:, IJE ,:,JN)
+
+PPHI(:,:,:,2)=MXM(PPHI(:,:,:,1))
+PPHI(:,:,:,3)=MYM(PPHI(:,:,:,1))
+PPHI(:,:,:,4)=MZM(PPHI(:,:,:,1))
+
+NULLIFY(TZFIELDS_ll)
+DO JN=2,4
+ PPHI(:,:,IKB-1,JN)=2*PPHI(:,:,IKB,JN)-PPHI(:,:,IKB+1,JN)
+ PPHI(:,:,IKE+1,JN)=2*PPHI(:,:,IKE,JN)-PPHI(:,:,IKE-1,JN)
+ PPHI(IIB-1,:,:,JN) = PPHI( IIB ,:,:,JN)
+ PPHI(IIE+1,:,:,JN) = PPHI( IIE ,:,:,JN)
+ PPHI(:,IJB-1,:,JN) = PPHI(:, IJB ,:,JN)
+ PPHI(:,IJE+1,:,JN) = PPHI(:, IJE ,:,JN)
+ENDDO
+
+PPHI(:,:,:,5)=MYM(PPHI(:,:,:,2))
+PPHI(:,:,:,6)=MXM(PPHI(:,:,:,4))
+PPHI(:,:,:,7)=MYM(PPHI(:,:,:,4))
+NULLIFY(TZFIELDS_ll)
+DO JN=5,7
+ PPHI(:,:,IKB-1,JN)=2*PPHI(:,:,IKB,JN)-PPHI(:,:,IKB+1,JN)
+ PPHI(:,:,IKE+1,JN)=2*PPHI(:,:,IKE,JN)-PPHI(:,:,IKE-1,JN)
+ PPHI(IIB-1,:,:,JN) = PPHI( IIB ,:,:,JN)
+ PPHI(IIE+1,:,:,JN) = PPHI( IIE ,:,:,JN)
+ PPHI(:,IJB-1,:,JN) = PPHI(:, IJB ,:,JN)
+ PPHI(:,IJE+1,:,JN) = PPHI(:, IJE ,:,JN)
+ENDDO
+WHERE (ABS(PPHI(:,:,:,:)).LT.XIBM_EPSI) PPHI(:,:,:,:)=2.*XIBM_EPSI
+
+
+ !COMPLETE PPHI ON THE HALO OF EACH SUBDOMAINS
+ DO JN=1,7
+ CALL ADD3DFIELD_ll(TZFIELDS_ll,PPHI(:,:,:,JN),'IBM_GENERLS::PPHI')
+ ENDDO
+ CALL UPDATE_HALO_ll(TZFIELDS_ll,IINFO_ll)
+ CALL CLEANLIST_ll(TZFIELDS_ll)
+
+ !
+ !-------------------------------------------------------------------------------
+ !
+ !**** X. DEALLOCATIONS/CLOSES
+ ! -----------------------
+ !
+ !DEALLOCATE(ZDP0PP0,ZDIST,ZC,ZSTEMP)
+ DEALLOCATE(ZC,ZSTEMP)
+ DEALLOCATE(ZXHATM,ZYHATM,ZZHATM)
+ !
+ RETURN
+ !
+CONTAINS
+ !
+ FUNCTION CROSSPRODUCT(PA,PB) RESULT(CROSS)
+ !
+ REAL, DIMENSION(3) :: CROSS
+ REAL, DIMENSION(3), INTENT(IN) :: PA, PB
+ CROSS(1) = PA(2) * PB(3) - PA(3) * PB(2)
+ CROSS(2) = PA(3) * PB(1) - PA(1) * PB(3)
+ CROSS(3) = PA(1) * PB(2) - PA(2) * PB(1)
+ END FUNCTION CROSSPRODUCT
+
+ FUNCTION SCALPRODUCT(PA,PB) RESULT(SCAL)
+ !
+ REAL :: SCAL
+ REAL, DIMENSION(3), INTENT(IN) :: PA, PB
+ SCAL = PA(1)*PB(1)+PA(2)*PB(2)+PA(3)*PB(3)
+ END FUNCTION SCALPRODUCT
+
+END SUBROUTINE IBM_GENERLS
diff --git a/src/MNH/ibm_prep_ls.f90 b/src/MNH/ibm_prep_ls.f90
index 22ba92cc5..063269439 100644
--- a/src/MNH/ibm_prep_ls.f90
+++ b/src/MNH/ibm_prep_ls.f90
@@ -37,10 +37,27 @@ SUBROUTINE IBM_PREP_LS(OIBM,HIBM_TYPE,PPHI)
!
! METHOD
! ------
- !**** Three main steps
- ! - read input ASCII files
- ! - Types of topography:
- ! IDEA : idealized obstacles (x,y coordinates)
+ !**** Types of topography:
+ ! 1)GENE : generalized obstacles (x,y coordinates)
+ ! => read the informations (triangles constituting the
+ ! faces of obstacles) from an .obj file.
+ ! The .obj file must have a particular organization:
+ ! a) A line with 'usemtl' indicates the 2 materials
+ ! of each side of the interface. Only the faces with
+ ! their external face in contact with the outside air
+ ! are read (mat2=air)
+ ! b) A line starting with 'v' indicates the location
+ ! (x,y,z coordinates) of a face vortex.
+ ! c) A line starting with 'f' indicates the vortices
+ ! constituting the face.
+ ! usemtl mat1:mat2
+ ! v xv1 yv1 zv1
+ ! v xv2 yv2 zv2
+ ! v xv3 yv3 zv3
+ ! v xv4 yv4 zv4
+ ! f 1 2 3
+ ! f 1 3 4
+ ! 2)IDEA : idealized obstacles (x,y coordinates)
!
! EXTERNAL
! --------
@@ -52,14 +69,17 @@ SUBROUTINE IBM_PREP_LS(OIBM,HIBM_TYPE,PPHI)
!
! REFERENCE
! ---------
+ ! For the generalized case, the method is based on '3D Distance from a
+ ! Point to a Triangle' a technical report from Mark W. Jones, University
+ ! of Wales Swansea [Jones (1995)].
!
- ! AUTHOR
+ ! AUTHORS
! ------
- ! Franck Auguste (CERFACS-AE)
+ ! Franck Auguste (CERFACS-AE), Tim Nagel (Météo-France)
!
! MODIFICATIONS
! -------------
- ! Original 01/01/2019
+ ! Original 01/06/2021
!
!------------------------------------------------------------------------------
!
@@ -86,6 +106,7 @@ SUBROUTINE IBM_PREP_LS(OIBM,HIBM_TYPE,PPHI)
! interface
USE MODI_SHUMAN
USE MODI_GDIV
+ USE MODI_IBM_GENERLS
USE MODI_IBM_IDEALRP
USE MODI_IBM_IDEALEE
!
@@ -93,6 +114,7 @@ SUBROUTINE IBM_PREP_LS(OIBM,HIBM_TYPE,PPHI)
USE MODD_CST
USE MODD_GRID_n
USE MODE_GRIDPROJ
+ USE MODE_MSG
!
IMPLICIT NONE
!
@@ -100,55 +122,30 @@ SUBROUTINE IBM_PREP_LS(OIBM,HIBM_TYPE,PPHI)
!
! 0.1 declarations of arguments
!
- LOGICAL ,INTENT(IN) :: OIBM ! flag for immersed boundary method
- CHARACTER(LEN=4) ,INTENT(IN) :: HIBM_TYPE ! switch generalized/idealised object
- REAL, DIMENSION(:,:,:,:) ,INTENT(INOUT) :: PPHI ! LS functions
+ LOGICAL ,INTENT(IN) :: OIBM ! flag for immersed boundary method
+ CHARACTER(LEN=4) ,INTENT(IN) :: HIBM_TYPE ! switch generalized/idealized object
+ REAL, DIMENSION(:,:,:,:) ,INTENT(INOUT) :: PPHI ! LS functions
!
!------------------------------------------------------------------------------
!
! 0.2 declaration of local variables
!
- INTEGER :: JN,JM,JNM,JL,JMM,JI,JJ,JK,JF,JV ! loop index
- INTEGER :: IIU,IJU,KII,KJJ
+ INTEGER :: JN,JM,JNM,JL,JMM,JI,JJ,JK,JF,JV ! loop index
+ INTEGER :: IIU,IJU
REAL :: ZX_MIN,ZX_MAX,ZY_MIN,ZY_MAX,DX_LOW,DY_LOW,DX_HIGH,DY_HIGH
- INTEGER :: JI2,JJ2,JK2,JI3,JJ3,JK3
- INTEGER :: JIM1,JIP1,JIM2,JIP2,JIM4,JIP4
- INTEGER :: JJM1,JJP1,JJM2,JJP2,JJM4,JJP4
- INTEGER :: JI2_MIN,JI2_MAX,JJ2_MIN,JJ2_MAX
- INTEGER :: IIB,IIE,IJB,IJE,IKB,IKE,ILOOP,JLOOP,KLOOP
- INTEGER :: IGRIB,IIBM_LEVEL,KIBM_LEVEL,IIBM_MIDDLE,KIBM_LEVEL2
- INTEGER :: KIII,KJJJ,KIIM1,KIIP1,KJJM1,KJJP1
- INTEGER :: ILUIBMIDEA,IRESPIBMGENE,ILUIBMGENE,IRESPIBMIDEA ! integers for open/read files
- INTEGER :: IIBM_NUMB_NODE_SURF ! number of surface points (generalized case)
- INTEGER :: IIBM_NUMB_TYPE_SURF ! number of surface type (idealized case)
- INTEGER :: IIBM_TYPE_SURF ! type of surfaces
- INTEGER :: IIBM_NUMB_SURF ! number of surfaces in each type
- REAL :: ZIBM_X1,ZIBM_X2,ZIBM_Y1,ZIBM_Y2,ZIBM_Z1,ZIBM_Z2 ! location of surface points for one object
+ INTEGER :: ILUIBMIDEA,IRESPIBMGENE,ILUIBMGENE,IRESPIBMIDEA ! integers for open/read files
+ INTEGER :: IIBM_NUMB_NODE_SURF ! number of surface points (generalized case)
+ INTEGER :: IIBM_NUMB_TYPE_SURF ! number of surface type (idealized case)
+ INTEGER :: IIBM_TYPE_SURF ! type of surfaces
+ INTEGER :: IIBM_NUMB_SURF ! number of surfaces in each type
+ REAL :: ZIBM_X1,ZIBM_X2,ZIBM_Y1,ZIBM_Y2,ZIBM_Z1,ZIBM_Z2 ! location of surface points for one object
REAL :: ZIBM_TYPE_SURF
- REAL, DIMENSION(:,:), ALLOCATABLE :: ZIBM_XYZ1,ZIBM_XYZ2 ! location of surface points for all object
- REAL, DIMENSION(:,:), ALLOCATABLE :: ZV1,ZV1_2,ZV2,ZV2_2,ZV3,ZV3_2
- REAL, DIMENSION(:,:), ALLOCATABLE :: NORM_FACES,NORM_FACES2
+ REAL, DIMENSION(:,:), ALLOCATABLE :: ZIBM_XYZ1,ZIBM_XYZ2 ! location of surface points for all object
+ REAL, DIMENSION(:,:), ALLOCATABLE :: ZV1,ZV1_2,ZV2,ZV2_2,ZV3,ZV3_2 ! face vectors
+ REAL, DIMENSION(:,:), ALLOCATABLE :: NORM_FACES,NORM_FACES2 ! norm of the faces
REAL, DIMENSION(:,:,:), ALLOCATABLE :: ZIBM_FACES,ZIBM_FACES2,ZIBM_FACES2b ! extremities of triangle faces for all object
- REAL :: XXX,YYY,ZZZ
- INTEGER :: IRESPIBMREAL,ILUIBMREAL ! reading/writing ASCII files
- REAL, DIMENSION(:,:,:) , ALLOCATABLE :: ZSURF,ZINDI,ZTMP2,ZTMP3 ! SSF and ISF functions + temporary arrays
- REAL, DIMENSION(:,:,:,:), ALLOCATABLE :: ZTMP,ZPHI ! temporary arrays
- REAL :: ZLAT,ZLON,ZHEI,ZIBM_HEI ! lat/lon/z coordinates
- INTEGER :: KNUM,KBOR,KIBM_BATMIN,KIBM_BATMAX ! index of buildings (BATI, REAL)
- CHARACTER(LEN=10) :: YHCOUNT1 ! reading/writing ASCII files
- CHARACTER(LEN=24) :: YHCOUNT2
- LOGICAL :: GHCOUNT3
- REAL :: ZXX,ZYY,ZXM2,ZYM2,ZII,ZJJ ! temporary values
- REAL :: ZMAX1,ZMAX2,ZMAX3
- REAL :: IIMAX,IJMAX,IKMAX
- REAL :: ZXX1,ZYY1,ZZZ1,ZXX2,ZYY2,ZZZ2
- REAL :: ZTES1,ZTES2,ZTES3,ZTES4,ZDIS,ZHIGH,ZHORI
- REAL :: ZLATMIN,ZLATMAX,ZLONMIN,ZLONMAX,ZXM2MIN,ZXM2MAX, ZYM2MIN, ZYM2MAX
- REAL :: SIGN1,SIGN2,SIGN3,SIGN4,ZHEI2
- REAL :: ZX1,ZY1,ZX2,ZY2,ZIND
TYPE(LIST_ll), POINTER :: TZFIELDS_ll
INTEGER :: IINFO_ll
- LOGICAL :: LCAEP,LAZF
CHARACTER(LEN=12) :: HFILEGENE, HFILEIDEA
CHARACTER(LEN=100) :: YSTRING,YSTRING2
INTEGER :: NS1,NS2,NS3,NS4,NS5,NS6
@@ -158,12 +155,24 @@ SUBROUTINE IBM_PREP_LS(OIBM,HIBM_TYPE,PPHI)
!------------------------------------------------------------------------------
!
! 0.3 Allocation
+ ILUIBMGENE = 42
ILUIBMIDEA = 43
+ HFILEGENE = "ibm_gene.obj"
HFILEIDEA = "ibm_idea.nam"
!
IIU = SIZE(XXHAT)
IJU = SIZE(XYHAT)
!
+ DX_LOW = XXHAT(2)-XXHAT(1)
+ DX_HIGH = XXHAT(IIU)-XXHAT(IIU-1)
+ DY_LOW = XYHAT(2)-XYHAT(1)
+ DY_HIGH = XYHAT(IJU)-XYHAT(IJU-1)
+ !
+ !Collect the face up to 10 gridsize out of the current processor
+ ZX_MIN = XXHAT(1)-10.*DX_LOW
+ ZX_MAX = XXHAT(IIU)+(11.)*DX_HIGH
+ ZY_MIN = XYHAT(1)-10.*DY_LOW
+ ZY_MAX = XYHAT(IJU)+(11.)*DY_HIGH
!
!------------------------------------------------------------------------------
!
@@ -171,6 +180,9 @@ SUBROUTINE IBM_PREP_LS(OIBM,HIBM_TYPE,PPHI)
! ----------------
!
! Read input files in order to compute interface location
+ ! - 'ibm_gene.obj' for generalized case
+ ! => read the informations (triangles constituting the
+ ! faces of obstacles) from an .obj file
! - 'm_ideal.nam' for idealized case
! (NUMB_NODE_SURF is the number of objects)
! (NUMB_TYPE_SURF is the number of surface types:
@@ -179,7 +191,181 @@ SUBROUTINE IBM_PREP_LS(OIBM,HIBM_TYPE,PPHI)
! ( NUMB_SURF is the objects number in each type)
!
!
- IF ((HIBM_TYPE=='IDEA')) THEN
+ !--------------------------------
+ !--------Generalized case--------
+ !--------------------------------
+ IF (HIBM_TYPE=='GENE') THEN
+ !
+ !Allocate the tables containing the vortices, the faces locations,
+ !the norms, which are needed to calculate the LSF
+ ALLOCATE(ZIBM_XYZ1(5400000,3))
+ ALLOCATE(ZIBM_FACES(3150000,3,3))
+ ALLOCATE(ZIBM_FACES2b(3150000,3,3))
+ ALLOCATE(NORM_FACES(3150000,3))
+ ALLOCATE(ZV1(3150000,3))
+ ALLOCATE(ZV2(3150000,3))
+ ALLOCATE(ZV3(3150000,3))
+ !
+ OPEN(ILUIBMGENE , FILE= HFILEGENE , IOSTAT=IRESPIBMGENE , STATUS='OLD')
+ !
+ JV=1
+ JF=0
+ JCOUNT=0
+ !
+ !Only the faces that are in contact with the air (external faces of
+ !the obstacles) are read.
+ DO
+ IF (IRESPIBMGENE/=0) EXIT
+ READ(UNIT=ILUIBMGENE,FMT='(A100)',IOSTAT=IRESPIBMGENE) YSTRING
+ NS1=LEN(TRIM(YSTRING))
+ IF (TRIM(YSTRING(1:7))=='usemtl') THEN
+ IF (TRIM(YSTRING(NS1-3:NS1))==':air') THEN
+ JN=1
+ ELSE
+ JN=0
+ ENDIF
+ ENDIF
+ IF (TRIM(YSTRING(1:2))=='v') THEN
+ NS2=INDEX(TRIM(YSTRING)," ",back=.true.)
+ NS3=LEN(TRIM(YSTRING(:NS2)))
+ NS4=INDEX(TRIM(YSTRING(:NS3))," ",back=.true.)
+ NS5=LEN(TRIM(YSTRING(:NS4)))
+ NS6=INDEX(TRIM(YSTRING(:NS5))," ",back=.true.)
+ READ(YSTRING(NS6:NS5) , *) ZIBM_XYZ1(JV,1)
+ READ(YSTRING(NS4:NS3) , *) ZIBM_XYZ1(JV,2)
+ READ(YSTRING(NS2:NS1) , *) ZIBM_XYZ1(JV,3)
+ !FIXME temporary spatial modification
+ ZIBM_XYZ1(JV,1) = ZIBM_XYZ1(JV,1) +200.
+ ZIBM_XYZ1(JV,2) = ZIBM_XYZ1(JV,2) +200.
+ JV=JV+1
+ ENDIF
+ IF (JN==1.AND.TRIM(YSTRING(1:2))=='f') THEN
+ NS2=INDEX(TRIM(YSTRING)," ",back=.true.)
+ NS3=LEN(TRIM(YSTRING(:NS2)))
+ NS4=INDEX(TRIM(YSTRING(:NS3))," ",back=.true.)
+ NS5=LEN(TRIM(YSTRING(:NS4)))
+ NS6=INDEX(TRIM(YSTRING(:NS5))," ",back=.true.)
+ READ(YSTRING(NS6:NS5) , *) ZN1
+ READ(YSTRING(NS4:NS3) , *) ZN2
+ READ(YSTRING(NS2:NS1) , *) ZN3
+ ! If the face extremities are far outside of the processor they are not read
+ IF (ZIBM_XYZ1(ZN1,1)<ZX_MIN.AND.ZIBM_XYZ1(ZN2,1)<ZX_MIN.AND.ZIBM_XYZ1(ZN3,1)<ZX_MIN) CYCLE
+ IF (ZIBM_XYZ1(ZN1,1)>ZX_MAX.AND.ZIBM_XYZ1(ZN2,1)>ZX_MAX.AND.ZIBM_XYZ1(ZN3,1)>ZX_MAX) CYCLE
+ IF (ZIBM_XYZ1(ZN1,2)<ZY_MIN.AND.ZIBM_XYZ1(ZN2,2)<ZY_MIN.AND.ZIBM_XYZ1(ZN3,2)<ZY_MIN) CYCLE
+ IF (ZIBM_XYZ1(ZN1,2)>ZY_MAX.AND.ZIBM_XYZ1(ZN2,2)>ZY_MAX.AND.ZIBM_XYZ1(ZN3,2)>ZY_MAX) CYCLE
+ JF=JF+1
+ !
+ ZIBM_FACES(JF,1,1) = ZIBM_XYZ1(ZN1,1)
+ ZIBM_FACES(JF,1,2) = ZIBM_XYZ1(ZN1,2)
+ ZIBM_FACES(JF,1,3) = ZIBM_XYZ1(ZN1,3)
+ ZIBM_FACES(JF,2,1) = ZIBM_XYZ1(ZN2,1)
+ ZIBM_FACES(JF,2,2) = ZIBM_XYZ1(ZN2,2)
+ ZIBM_FACES(JF,2,3) = ZIBM_XYZ1(ZN2,3)
+ ZIBM_FACES(JF,3,1) = ZIBM_XYZ1(ZN3,1)
+ ZIBM_FACES(JF,3,2) = ZIBM_XYZ1(ZN3,2)
+ ZIBM_FACES(JF,3,3) = ZIBM_XYZ1(ZN3,3)
+ !
+ ZNA(1) = ZIBM_FACES(JF,1,1)-ZIBM_FACES(JF,2,1)
+ ZNA(2) = ZIBM_FACES(JF,1,2)-ZIBM_FACES(JF,2,2)
+ ZNA(3) = ZIBM_FACES(JF,1,3)-ZIBM_FACES(JF,2,3)
+ ZNB(1) = ZIBM_FACES(JF,1,1)-ZIBM_FACES(JF,3,1)
+ ZNB(2) = ZIBM_FACES(JF,1,2)-ZIBM_FACES(JF,3,2)
+ ZNB(3) = ZIBM_FACES(JF,1,3)-ZIBM_FACES(JF,3,3)
+ !Elimination of '1D' faces
+ IF (ZNA(1)==0..AND.ZNA(2)==0..AND.ZNA(3)==0.) CYCLE
+ IF (ZNB(1)==0..AND.ZNB(2)==0..AND.ZNB(3)==0.) CYCLE
+ IF (ZNA(2)==0..AND.ZNA(3)==0..AND.ZNB(2)==0..AND.ZNB(3)==0.) CYCLE
+ IF (ZNA(1)==ZNB(1).AND.ZNA(2)==ZNB(2).AND.ZNA(3)==ZNB(3)) CYCLE
+ JCOUNT=JCOUNT+1
+ NORM_FACES(JCOUNT,:)= CROSSPRODUCT(ZNA,ZNB)
+ ZIBM_FACES2b(JCOUNT,:,:)=ZIBM_FACES(JF,:,:)
+ !
+ !Equation (6) of Jones (1995)
+ ZV1(JCOUNT,1) = (ZIBM_FACES(JF,1,1)-ZIBM_FACES(JF,2,1))/ SQRT((ZIBM_FACES(JF,1,1)-ZIBM_FACES(JF,2,1))**2 + &
+ (ZIBM_FACES(JF,1,2)-ZIBM_FACES(JF,2,2))**2 +(ZIBM_FACES(JF,1,3)-ZIBM_FACES(JF,2,3))**2)+ &
+ (ZIBM_FACES(JF,1,1)-ZIBM_FACES(JF,3,1))/ SQRT((ZIBM_FACES(JF,1,1)-ZIBM_FACES(JF,3,1))**2 + &
+ (ZIBM_FACES(JF,1,2)-ZIBM_FACES(JF,3,2))**2 +(ZIBM_FACES(JF,1,3)-ZIBM_FACES(JF,3,3))**2)
+ !
+ ZV1(JCOUNT,2) = (ZIBM_FACES(JF,1,2)-ZIBM_FACES(JF,2,2))/ SQRT((ZIBM_FACES(JF,1,1)-ZIBM_FACES(JF,2,1))**2 + &
+ (ZIBM_FACES(JF,1,2)-ZIBM_FACES(JF,2,2))**2 +(ZIBM_FACES(JF,1,3)-ZIBM_FACES(JF,2,3))**2)+ &
+ (ZIBM_FACES(JF,1,2)-ZIBM_FACES(JF,3,2))/ SQRT((ZIBM_FACES(JF,1,1)-ZIBM_FACES(JF,3,1))**2 + &
+ (ZIBM_FACES(JF,1,2)-ZIBM_FACES(JF,3,2))**2 +(ZIBM_FACES(JF,1,3)-ZIBM_FACES(JF,3,3))**2)
+ !
+ ZV1(JCOUNT,3) = (ZIBM_FACES(JF,1,3)-ZIBM_FACES(JF,2,3))/ SQRT((ZIBM_FACES(JF,1,1)-ZIBM_FACES(JF,2,1))**2 + &
+ (ZIBM_FACES(JF,1,2)-ZIBM_FACES(JF,2,2))**2 +(ZIBM_FACES(JF,1,3)-ZIBM_FACES(JF,2,3))**2)+ &
+ (ZIBM_FACES(JF,1,3)-ZIBM_FACES(JF,3,3))/ SQRT((ZIBM_FACES(JF,1,1)-ZIBM_FACES(JF,3,1))**2 + &
+ (ZIBM_FACES(JF,1,2)-ZIBM_FACES(JF,3,2))**2 +(ZIBM_FACES(JF,1,3)-ZIBM_FACES(JF,3,3))**2)
+ !
+ ZV2(JCOUNT,1) = (ZIBM_FACES(JF,2,1)-ZIBM_FACES(JF,3,1))/ SQRT((ZIBM_FACES(JF,2,1)-ZIBM_FACES(JF,3,1))**2 + &
+ (ZIBM_FACES(JF,2,2)-ZIBM_FACES(JF,3,2))**2 +(ZIBM_FACES(JF,2,3)-ZIBM_FACES(JF,3,3))**2)+ &
+ (ZIBM_FACES(JF,2,1)-ZIBM_FACES(JF,1,1))/ SQRT((ZIBM_FACES(JF,2,1)-ZIBM_FACES(JF,1,1))**2 + &
+ (ZIBM_FACES(JF,2,2)-ZIBM_FACES(JF,1,2))**2 +(ZIBM_FACES(JF,2,3)-ZIBM_FACES(JF,1,3))**2)
+ !
+ ZV2(JCOUNT,2) = (ZIBM_FACES(JF,2,2)-ZIBM_FACES(JF,3,2))/ SQRT((ZIBM_FACES(JF,2,1)-ZIBM_FACES(JF,3,1))**2 + &
+ (ZIBM_FACES(JF,2,2)-ZIBM_FACES(JF,3,2))**2 +(ZIBM_FACES(JF,2,3)-ZIBM_FACES(JF,3,3))**2)+ &
+ (ZIBM_FACES(JF,2,2)-ZIBM_FACES(JF,1,2))/ SQRT((ZIBM_FACES(JF,2,1)-ZIBM_FACES(JF,1,1))**2 + &
+ (ZIBM_FACES(JF,2,2)-ZIBM_FACES(JF,1,2))**2 +(ZIBM_FACES(JF,2,3)-ZIBM_FACES(JF,1,3))**2)
+ !
+ ZV2(JCOUNT,3) = (ZIBM_FACES(JF,2,3)-ZIBM_FACES(JF,3,3))/ SQRT((ZIBM_FACES(JF,2,1)-ZIBM_FACES(JF,3,1))**2 + &
+ (ZIBM_FACES(JF,2,2)-ZIBM_FACES(JF,3,2))**2 +(ZIBM_FACES(JF,2,3)-ZIBM_FACES(JF,3,3))**2)+ &
+ (ZIBM_FACES(JF,2,3)-ZIBM_FACES(JF,1,3))/ SQRT((ZIBM_FACES(JF,2,1)-ZIBM_FACES(JF,1,1))**2 + &
+ (ZIBM_FACES(JF,2,2)-ZIBM_FACES(JF,1,2))**2 +(ZIBM_FACES(JF,2,3)-ZIBM_FACES(JF,1,3))**2)
+ !
+ ZV3(JCOUNT,1) = (ZIBM_FACES(JF,3,1)-ZIBM_FACES(JF,1,1))/ SQRT((ZIBM_FACES(JF,3,1)-ZIBM_FACES(JF,1,1))**2 + &
+ (ZIBM_FACES(JF,3,2)-ZIBM_FACES(JF,1,2))**2 +(ZIBM_FACES(JF,3,3)-ZIBM_FACES(JF,1,3))**2)+ &
+ (ZIBM_FACES(JF,3,1)-ZIBM_FACES(JF,2,1))/ SQRT((ZIBM_FACES(JF,3,1)-ZIBM_FACES(JF,2,1))**2 + &
+ (ZIBM_FACES(JF,3,2)-ZIBM_FACES(JF,2,2))**2 +(ZIBM_FACES(JF,3,3)-ZIBM_FACES(JF,2,3))**2)
+ !
+ ZV3(JCOUNT,2) = (ZIBM_FACES(JF,3,2)-ZIBM_FACES(JF,1,2))/ SQRT((ZIBM_FACES(JF,3,1)-ZIBM_FACES(JF,1,1))**2 + &
+ (ZIBM_FACES(JF,3,2)-ZIBM_FACES(JF,1,2))**2 +(ZIBM_FACES(JF,3,3)-ZIBM_FACES(JF,1,3))**2)+ &
+ (ZIBM_FACES(JF,3,2)-ZIBM_FACES(JF,2,2))/ SQRT((ZIBM_FACES(JF,3,1)-ZIBM_FACES(JF,2,1))**2 + &
+ (ZIBM_FACES(JF,3,2)-ZIBM_FACES(JF,2,2))**2 +(ZIBM_FACES(JF,3,3)-ZIBM_FACES(JF,2,3))**2)
+ !
+ ZV3(JCOUNT,3) = (ZIBM_FACES(JF,3,3)-ZIBM_FACES(JF,1,3))/ SQRT((ZIBM_FACES(JF,3,1)-ZIBM_FACES(JF,1,1))**2 + &
+ (ZIBM_FACES(JF,3,2)-ZIBM_FACES(JF,1,2))**2 +(ZIBM_FACES(JF,3,3)-ZIBM_FACES(JF,1,3))**2)+ &
+ (ZIBM_FACES(JF,3,3)-ZIBM_FACES(JF,2,3))/ SQRT((ZIBM_FACES(JF,3,1)-ZIBM_FACES(JF,2,1))**2 + &
+ (ZIBM_FACES(JF,3,2)-ZIBM_FACES(JF,2,2))**2 +(ZIBM_FACES(JF,3,3)-ZIBM_FACES(JF,2,3))**2)
+ !
+ ENDIF
+ !
+ IF (JN==1.AND.TRIM(YSTRING(1:2))=='vn') THEN
+ WRITE(*,*) '*****************************'
+ WRITE(*,*) '***** vn found in .obj ******'
+ WRITE(*,*) '***** Unable to read it *****'
+ WRITE(*,*) '**** (stopped execution) ****'
+ WRITE(*,*) '*****************************'
+ call Print_msg( NVERB_FATAL, 'GEN', 'IBM_PREP_LS', 'Unable to read vn found in .obj' )
+ ENDIF
+ !
+ IF (JN==1.AND.TRIM(YSTRING(1:2))=='vt') THEN
+ WRITE(*,*) '*****************************'
+ WRITE(*,*) '***** vt found in .obj ******'
+ WRITE(*,*) '***** Unable to read it *****'
+ WRITE(*,*) '**** (stopped execution) ****'
+ WRITE(*,*) '*****************************'
+ call Print_msg( NVERB_FATAL, 'GEN', 'IBM_PREP_LS', 'Unable to read vt found in .obj' )
+ ENDIF
+ !
+ END DO
+ !
+ ALLOCATE(ZIBM_FACES2(JCOUNT,3,3))
+ ALLOCATE(NORM_FACES2(JCOUNT,3))
+ ALLOCATE(ZV1_2(JCOUNT,3))
+ ALLOCATE(ZV2_2(JCOUNT,3))
+ ALLOCATE(ZV3_2(JCOUNT,3))
+ !
+ NORM_FACES2 = NORM_FACES(:JCOUNT,:)
+ ZV1_2 = ZV1(:JCOUNT,:)
+ ZV2_2 = ZV2(:JCOUNT,:)
+ ZV3_2 = ZV3(:JCOUNT,:)
+ ZIBM_FACES2 = ZIBM_FACES2b(:JCOUNT,:,:)
+ !
+ ENDIF
+ !
+ !----------------------------
+ !---------Idealized case-----
+ !----------------------------
+ IF (HIBM_TYPE=='IDEA') THEN
!
OPEN(ILUIBMIDEA , FILE= HFILEIDEA , IOSTAT=IRESPIBMIDEA , FORM='FORMATTED' , &
STATUS='OLD', ACCESS='SEQUENTIAL', ACTION='READ')
@@ -216,16 +402,41 @@ SUBROUTINE IBM_PREP_LS(OIBM,HIBM_TYPE,PPHI)
! -------------
!
! Computations of volumic fraction (VF) and Level Set function (LS) for all kinds of initialization
- ! idealized shape => construction of VF/LS function using analytical
+ ! generalized shape => construction of LS function (z<z_interface <=> phi>0)
+ ! the method is based on '3D Distance from a Point to a Triangle' [Jones (1995)].
+ ! => conversion of LS function to VF function (Sussman, JCP (1994)
+ ! idealized shape => construction of VF/LS function using analytical
! locations of interface (ellipsoidal/parallelepipedic shapes)
!
- IF ((HIBM_TYPE=='IDEA')) then
+ IF (HIBM_TYPE=='GENE') THEN
+ CALL IBM_GENERLS(ZIBM_FACES2,NORM_FACES2,ZV1_2,ZV2_2,ZV3_2,ZX_MIN,ZY_MIN,ZX_MAX,ZY_MAX,PPHI)
+ ENDIF
+ !
+ IF (HIBM_TYPE=='IDEA') then
DO JN=1,JNM
- !
IF (abs(ZIBM_XYZ2(JN,7)-1.).lt.XIBM_EPSI) CALL IBM_IDEALRP(JN,ZIBM_XYZ2,PPHI)
IF (abs(ZIBM_XYZ2(JN,7)-2.).lt.XIBM_EPSI) CALL IBM_IDEALEE(JN,ZIBM_XYZ2,PPHI)
ENDDO
- !
ENDIF
!
+CONTAINS
+ !
+ FUNCTION CROSSPRODUCT(PA,PB) RESULT(CROSS)
+ !
+ REAL, DIMENSION(3) :: CROSS
+ REAL :: VAL
+ REAL, DIMENSION(3), INTENT(IN) :: PA, PB
+ !
+ CROSS(1) = PA(2) * PB(3) - PA(3) * PB(2)
+ CROSS(2) = PA(3) * PB(1) - PA(1) * PB(3)
+ CROSS(3) = PA(1) * PB(2) - PA(2) * PB(1)
+ !
+ VAL = (CROSS(1)**2+CROSS(2)**2+CROSS(3)**2)**(0.5)
+ !
+ CROSS(1) = CROSS(1)/VAL
+ CROSS(2) = CROSS(2)/VAL
+ CROSS(3) = CROSS(3)/VAL
+ !
+ END FUNCTION CROSSPRODUCT
+ !
END SUBROUTINE IBM_PREP_LS
--
GitLab