Skip to content
Snippets Groups Projects
Commit 0760af85 authored by RODIER Quentin's avatar RODIER Quentin
Browse files

T.Nagel 23/06/2021: IBM Object file reading to compute the Level Set function for any surface

parent af743026
No related branches found
No related tags found
No related merge requests found
Expand all Hide whitespace changes
Inline Side-by-side
src/MNH/ibm_generls.f90 0 → 100644
+ 562
− 0
View file @ 0760af85
This diff is collapsed.
src/MNH/ibm_prep_ls.f90
+ 262
− 51
View file @ 0760af85
......@@ -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
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment

Documentation | Mentions légales | CGU | Accessibilité : non conforme