Skip to content

Commit

Permalink
Merge pull request #62 from hopr-framework/feature.mesh.shift
Browse files Browse the repository at this point in the history
Apply Global Shift to Mesh Coordinates
  • Loading branch information
kopperp authored Jul 12, 2024
2 parents ca49ab0 + 8ad6507 commit f7c6005
Show file tree
Hide file tree
Showing 3 changed files with 184 additions and 13 deletions.
52 changes: 40 additions & 12 deletions src/mesh/mesh.f90
Original file line number Diff line number Diff line change
Expand Up @@ -204,6 +204,11 @@ SUBROUTINE InitMesh()
doScale = (ABS(MeshScale-1.).GT.PP_RealTolerance)
SpaceQuandt = GETREAL('SpaceQuandt','0.1')

! ATTENTION: Scaling is always applied first, then shifting! So `MeshShift` should always be given in scaled coordinates!
postShift = GETLOGICAL('postShiftMesh','.TRUE.') ! apply shifting either after readin or before output
MeshShift = GETREALARRAY('meshShift',3,'0.,0.,0.') ! shifting vector applied to node coordinates during read in
doShift = (NORM2(MeshShift).GT.PP_RealTolerance)

! Curved
IF(useCurveds) THEN
! If curved mesh is read in (e.g. CGNS/GMSH/HDF5) curved boundaries can be rebuilt using the linear mesh
Expand Down Expand Up @@ -574,8 +579,10 @@ SUBROUTINE fillMesh()
CALL abort(__STAMP__,'Not known how to construct mesh')
END SELECT

! apply meshscale after readin
IF(doScale.AND..NOT.postScale) CALL ApplyMeshScale(FirstElem)
! apply meshscale and meshshift after readin
IF( (doShift.AND..NOT.postShift) .OR. (doScale.AND..NOT.postScale) ) THEN
CALL ApplyMeshScale(FirstElem, doMeshScale=(doScale.AND..NOT.postScale), doMeshShift=(doShift.AND..NOT.postShift))
END IF

! get element trafo and ensure right-handed coordinate system
Elem=>FirstElem
Expand Down Expand Up @@ -689,7 +696,9 @@ SUBROUTINE fillMesh()
CALL abort(__STAMP__,'Splitted meshes can only be read in star cd format or cgns!')
END SELECT
CALL Connect2DMesh(FirstSplitElem)
IF(doScale.AND..NOT.postScale) CALL ApplyMeshScale(FirstSplitElem)
IF( (doShift.AND..NOT.postShift) .OR. (doScale.AND..NOT.postScale) ) THEN
CALL ApplyMeshScale(FirstSplitElem, doMeshScale=(doScale.AND..NOT.postScale), doMeshShift=(doShift.AND..NOT.postShift))
END IF
CALL splitToSpline()
CALL curvedEdgesToSurf(keepExistingCurveds=.FALSE.)
CASE(4) !!!only in combination with icem cgns meshes!!!!
Expand Down Expand Up @@ -781,6 +790,9 @@ SUBROUTINE fillMesh()

! apply meshscale before output (default)
IF(doScale.AND.postScale) CALL ApplyMeshScale(FirstElem)
IF( (doShift.AND.postShift) .OR. (doScale.AND.postScale) ) THEN
CALL ApplyMeshScale(FirstElem, doMeshScale=(doScale.AND.postScale), doMeshShift=(doShift.AND.postShift))
END IF
IF(DebugVisu) THEN
CALL netVisu() ! visualize mesh
CALL BCVisu() ! visualize BC
Expand Down Expand Up @@ -1090,18 +1102,21 @@ SUBROUTINE CheckNodeConnectivity()
END DO
END SUBROUTINE CheckNodeConnectivity

SUBROUTINE ApplyMeshScale(StartElem)
SUBROUTINE ApplyMeshScale(StartElem,doMeshScale,doMeshShift)
!===================================================================================================================================
! Apply meshscale to existing mesh (called either directly after readin or before output)
! Apply mesh scaling factor and shift to existing mesh (called either directly after readin or before output)
! ATTENTION: Scaling is always applied first, then shifting! So `MeshShift` should always be given in scaled coordinates!
!===================================================================================================================================
! MODULES
USE MOD_Mesh_Vars, ONLY:tElem,tSide,tEdge
USE MOD_Mesh_Vars, ONLY:N,MeshScale
USE MOD_Mesh_Vars, ONLY:N,MeshScale,MeshShift
! IMPLICIT VARIABLE HANDLING
IMPLICIT NONE
!-----------------------------------------------------------------------------------------------------------------------------------
! INPUT VARIABLES
TYPE(tElem),POINTER,INTENT(IN) :: StartElem ! ?
LOGICAL,OPTIONAL,INTENT(IN) :: doMeshScale ! Flag to indicate whether `MeshScale` should be applied
LOGICAL,OPTIONAL,INTENT(IN) :: doMeshShift ! Flag to indicate whether `MeshShift` should be applied
!-----------------------------------------------------------------------------------------------------------------------------------
! OUTPUT VARIABLES
!-----------------------------------------------------------------------------------------------------------------------------------
Expand All @@ -1110,21 +1125,34 @@ SUBROUTINE ApplyMeshScale(StartElem)
TYPE(tSide),POINTER :: Side ! ?
TYPE(tEdge),POINTER :: Edge ! ?
INTEGER :: iNode,iEdge ! ?
LOGICAL :: doScale,doShift ! local flags for optionals
REAL :: MeshScale_loc ! Contains either MeshScale or Identity operation (*1.)
REAL :: MeshShift_loc(3)! Contains either MeshShift or Identity operation (+0.)
!===================================================================================================================================
! Set MeshScale_loc=MeshScale only if doScale=.TRUE., else use identity operation (*1.)
doScale = .FALSE.
IF (PRESENT(doMeshScale)) doScale = doMeshScale
MeshScale_loc = MERGE(MeshScale,1.,doScale)

! Set MeshShift_loc=MeshShift only if doShift=.TRUE., else use identity operation (+0.)
doShift = .FALSE.
IF (PRESENT(doMeshShift)) doShift = doMeshShift
MeshShift_loc(:) = MERGE(MeshShift(:),0.,doShift)

! make all nodes unique
Elem=>StartElem
DO WHILE(ASSOCIATED(Elem))
! elem nodes
DO iNode=1,Elem%nNodes
IF(Elem%node(iNode)%np%tmp.NE.-11)THEN
Elem%node(iNode)%np%x=Elem%node(iNode)%np%x*MeshScale
Elem%node(iNode)%np%x=(Elem%node(iNode)%np%x*MeshScale_loc)+MeshShift_loc
Elem%node(iNode)%np%tmp=-11
END IF
END DO
! curved elem nodes
DO iNode=1,Elem%nCurvedNodes
IF(Elem%curvedNode(iNode)%np%tmp.NE.-11)THEN
Elem%curvedNode(iNode)%np%x=Elem%curvedNode(iNode)%np%x*MeshScale
Elem%curvedNode(iNode)%np%x=(Elem%curvedNode(iNode)%np%x*MeshScale_loc)+MeshShift_loc
Elem%curvedNode(iNode)%np%tmp=-11
END IF
END DO
Expand All @@ -1134,14 +1162,14 @@ SUBROUTINE ApplyMeshScale(StartElem)
! side nodes
DO iNode=1,Side%nNodes
IF(side%node(iNode)%np%tmp.NE.-11)THEN
side%node(iNode)%np%x=side%node(iNode)%np%x*MeshScale
side%node(iNode)%np%x=(side%node(iNode)%np%x*MeshScale_loc)+MeshShift_loc
side%node(iNode)%np%tmp=-11
END IF
END DO
! curved side nodes
DO iNode=1,Side%nCurvedNodes
IF(side%curvedNode(iNode)%np%tmp.NE.-11)THEN
side%curvedNode(iNode)%np%x=side%curvedNode(iNode)%np%x*MeshScale
side%curvedNode(iNode)%np%x=(side%curvedNode(iNode)%np%x*MeshScale_loc)+MeshShift_loc
side%curvedNode(iNode)%np%tmp=-11
END IF
END DO
Expand All @@ -1153,15 +1181,15 @@ SUBROUTINE ApplyMeshScale(StartElem)
DO iNode=1,2
IF(.NOT.ASSOCIATED(edge%node(iNode)%np)) CYCLE
IF(edge%Node(iNode)%np%tmp.NE.-11)THEN
edge%node(iNode)%np%x=edge%node(iNode)%np%x*MeshScale
edge%node(iNode)%np%x=(edge%node(iNode)%np%x*MeshScale_loc)+MeshShift_loc
edge%node(iNode)%np%tmp=-11
END IF
END DO
! curved edge nodes
IF(ASSOCIATED(edge%curvedNode))THEN
DO iNode=1,N+1
IF(edge%curvedNode(iNode)%np%tmp.NE.-11)THEN
edge%curvedNode(iNode)%np%x=edge%curvedNode(iNode)%np%x*MeshScale
edge%curvedNode(iNode)%np%x=(edge%curvedNode(iNode)%np%x*MeshScale_loc)+MeshShift_loc
edge%curvedNode(iNode)%np%tmp=-11
END IF
END DO
Expand Down
4 changes: 3 additions & 1 deletion src/mesh/mesh_vars.f90
Original file line number Diff line number Diff line change
Expand Up @@ -220,9 +220,11 @@ MODULE MOD_Mesh_Vars
TYPE(tElem) ,POINTER :: firstElem ! pointer to first element in order to start a loop
TYPE(tElem) ,POINTER :: firstSplitElem ! pointer to first element in splitted elem list for curveds
LOGICAL :: doScale ! scaling factor gt realtolerance
LOGICAL :: preScale ! apply scaling after readin or before output
LOGICAL :: postScale ! apply scaling after readin or before output
REAL :: MeshScale ! scaling factor applied to Node Coordinates during read in
LOGICAL :: doShift ! scaling factor gt realtolerance
LOGICAL :: postShift ! apply scaling after readin or before output
REAL :: MeshShift(3) ! shifting vector applied to Node Coordinates during read in
REAL :: SpaceQuandt ! Characteristic length in the mesh. Used as tolerance
REAL :: minDX ! smallest edge length
REAL :: maxDX(3) ! Used for search mesh
Expand Down
141 changes: 141 additions & 0 deletions tutorials/1-06-curved-postdeform/parameter_cylinderinchannel.ini
Original file line number Diff line number Diff line change
@@ -0,0 +1,141 @@
DEFVAR=(INT): ia = 004 ! no. elems in angular direction top/bottom
DEFVAR=(INT): ial= 002 ! no. elems in angular direction left (inflow)
DEFVAR=(INT): iar= 005 ! no. elems in angular direction right (wake)
DEFVAR=(INT): ir = 006 ! no. elems in radial direction
DEFVAR=(INT): iz = 001 ! no. elems in z
DEFVAR=(INT): iw = 012 ! no. elems in wake
DEFVAR=(REAL): fr = 1.2 ! radial stretching factor in ring
DEFVAR=(REAL): fw = 1.1 ! streamwise stretching factor in wake
DEFVAR=(REAL): rm = 1.0 ! middle square dim
DEFVAR=(REAL): r0 = 5.93969696196699920 ! = sqrt(2.)*ymax/(PostDeform_R0)
DEFVAR=(REAL): r1 = 5.65685424949238019 ! = sqrt(2.)*ymin/(PostDeform_R0)
DEFVAR=(REAL): lx = 56.5685424949238020 ! = sqrt(2.)*(xtotal-2.)/(PostDeform_R0)
DEFVAR=(REAL): lz = 1. ! length of domain in z
!================================================================================================================================= !
! 2D cylinder in channel following (all measures factor 10 larger) 2D-2 case by:
! - "Benchmark Computations of Laminar Flow Around a Cylinder", Schäfer, Turek, 1996.
! (http://www.featflow.de/en/benchmarks/cfdbenchmarking/flow/dfg_benchmark2_re100.html)
!
! ymin ymin xmax - 2*ymin
! |<-------->|<--------->|<------------------------------------->|
! -- .--------------------------------------------------------------.
! ^ | ^ | |
! ymax | | | ---- | |
! | | **|** ^ | |
! v | ** | ** R0| | |
! -- | ** 0------> | | wake |
! ^ | ** ** v | |
! ymin | | ***** ---- | |
! v | | |
! -- '--------------------------------------------------------------'
! |<------------------------------------------------------------>|
! xtotal
!
! Explanation:
! ----------------
!
! The cylinder is defined in [-1,1]^1 with diameter=1, and is created using a radius of PostDeform_R0=0.5. The r0 and r1 values,
! i.e. the maximum/minimum y-values of the mesh (r0 and r1), are calculated as follows to account for the PostDeform:
! r0 = sqrt(2.)*ymax/(PostDeform_R0) ,with ymax = 2.1
! r1 = sqrt(2.)*ymin/(PostDeform_R0) ,with ymin = 2.0
!
! Similarly, the maximum x coordinate of the domain "lx" (i.e. from 0 to end of domain) can be computed as:
! lx = sqrt(2.)*(xtotal-ymin)/(PostDeform_R0) ,with xtotal = 22
!
! Note that the cylinder is thus embedded in a rectangle that is longer in positive y-axis side as detailed in the reference.
! The mesh consits of 5 zones. 4 to build the first part around the cylinder:
! - left: [-r1,-r1] to [-rm,-rm]
! - right: [ rm,-rm] to [ r1, r0]
! - upper: [-rm, rm] to [ r1, r0]
! - lower: [-r1,-r1] to [ rm,-rm]
!
! And lastly the "wake" part that is a single structured block that extends the mesh in streamwise direction to x-coordinate "lx".
! - wake: [r1,-r1,0] to [lx,r0,0]
!
! The final step is to shift the final mesh by vector:
! MeshShift = (/ ymin, -(ymax-ymin)/2, 0 /)
! such that the domain is symmetric around the y-axis and the domain is within [0,22.] x [-2.05,2.05] x [0,1.]
! In order to obtain domain from paper with y\in[0,H], use instead:
! MeshShift = (/ ymin, ymin, 0/)
!
!================================================================================================================================= !
! OUTPUT
!================================================================================================================================= !
ProjectName = CylinderInChannel ! Name of output files
Debugvisu = T ! Visualize mesh and boundary conditions (tecplot ascii)
checkElemJacobians=T
NVisu=12

!================================================================================================================================= !
! MESH
!================================================================================================================================= !
Mode =1 ! Mode for Cartesian boxes
nZones =5 ! number of boxes
! left
Corner =(/-r1,-r1,0. ,,-rm,-rm,0. ,,-rm,rm,0. ,, -r1,r0,0.,, -r1,-r1,lz ,,-rm,-rm,lz ,,-rm,rm,lz ,, -r1,r0,lz /)
nElems =(/ir,ial,iz/) ! number of elements in each direction
BCIndex =(/1,0,7,0,5,6/) ! Indices of Boundary Conditions
factor =(/-fr,1.,1./) ! stretching
elemtype =108 ! element type (108: Hexahedral)
! right
Corner =(/rm,-rm,0. ,,r1,-r1,0. ,,r1,r0,0. ,, rm,rm,0.,, rm,-rm,lz ,,r1,-r1,lz ,,r1,r0,lz ,, rm,rm,lz /)
nElems =(/ir,iar,iz/) ! number of elements in each direction
!BCIndex =(/1,0,3,0,7,6/) ! Indices of Boundary Conditions
BCIndex =(/1,0,0,0,7,6/) ! Indices of Boundary Conditions
elemtype =108 ! element type (108: Hexahedral)
factor =(/fr,1.,1./) ! stretching
! upper
Corner =(/-rm,rm,0. ,,rm,rm,0. ,,r1,r0,0. ,, -r1,r0,0.,, -rm,rm,lz ,,rm,rm,lz ,,r1,r0,lz ,, -r1,r0,lz /)
nElems =(/ia,ir,iz/) ! number of elements in each direction
BCIndex =(/1,7,0,4,0,6/) ! Indices of Boundary Conditions
elemtype =108 ! element type (108: Hexahedral)
factor =(/1.,fr,1./) ! stretching
! lower
Corner =(/-r1,-r1,0. ,,r1,-r1,0. ,,rm,-rm,0. ,, -rm,-rm,0.,, -r1,-r1,lz ,,r1,-r1,lz ,,rm,-rm,lz ,, -rm,-rm,lz /)
nElems =(/ia,ir,iz/) ! number of elements in each direction
BCIndex =(/1,2,0,7,0,6/) ! Indices of Boundary Conditions
elemtype =108 ! element type (108: Hexahedral)
factor =(/1.,-fr,1./) ! stretching
! wake
Corner =(/r1,-r1,0. ,,lx,-r1,0. ,,lx,r0,0. ,, r1,r0,0.,, r1,-r1,lz ,,lx,-r1,lz ,,lx,r0,lz ,, r1,r0,lz /)
nElems =(/iw,iar,iz/) ! number of elements in each direction
BCIndex =(/1,2,3,4,0,6/) ! Indices of Boundary Conditions
factor =(/fw,1.,1./) ! stretching
elemtype =108 ! element type (108: Hexahedral)

useCurveds =T
BoundaryOrder=6

!================================================================================================================================= !
! BOUNDARY CONDITIONS
!================================================================================================================================= !
BoundaryName=BC_zminus ! BC index 1 (from position in parameterfile)
BoundaryType=(/1,0,0,1/) ! (/ Type, curveIndex, State, alpha /)
BoundaryName=BC_wall_lower ! BC index 2
BoundaryType=(/3,0,0,0/)
BoundaryName=BC_xplus ! BC index 3
BoundaryType=(/24,0,1,0/)
BoundaryName=BC_wall_upper ! BC index 4
BoundaryType=(/3,0,0,0/)
BoundaryName=BC_xminus ! BC index 5
BoundaryType=(/2,0,0,0/)
BoundaryName=BC_zplus ! BC index 6
BoundaryType=(/1,0,0,-1/)
BoundaryName=BC_wall_cylinder ! BC index 7
BoundaryType=(/3,0,0,0/)
vv=(/0.,0.,lz/)

MeshPostDeform=3 ! deforms [-1,1]^2 to a cylinder with radius Postdeform_R0 and goes back to box [-4,4]^2
PostDeform_R0=0.5

MeshShift = (/2.,-0.05,0./) ! Shift final mesh such that it is within [0,22.] x [-2.05,2.05] x [0,1.]

!================================================================================================================================= !
! Z Correction
!================================================================================================================================= !
OrientZ = T
dozcorrection = T
zPeriodic = T
zLength = lz
zstart = 0.
nElemsZ = iz

0 comments on commit f7c6005

Please sign in to comment.