diff --git a/src/mesh/mesh.f90 b/src/mesh/mesh.f90 index ab8e451..0c26981 100644 --- a/src/mesh/mesh.f90 +++ b/src/mesh/mesh.f90 @@ -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 @@ -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 @@ -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!!!! @@ -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 @@ -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 !----------------------------------------------------------------------------------------------------------------------------------- @@ -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 @@ -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 @@ -1153,7 +1181,7 @@ 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 @@ -1161,7 +1189,7 @@ SUBROUTINE ApplyMeshScale(StartElem) 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 diff --git a/src/mesh/mesh_vars.f90 b/src/mesh/mesh_vars.f90 index 29ff2ec..8e8d2be 100644 --- a/src/mesh/mesh_vars.f90 +++ b/src/mesh/mesh_vars.f90 @@ -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 diff --git a/tutorials/1-06-curved-postdeform/parameter_cylinderinchannel.ini b/tutorials/1-06-curved-postdeform/parameter_cylinderinchannel.ini new file mode 100644 index 0000000..239bac4 --- /dev/null +++ b/tutorials/1-06-curved-postdeform/parameter_cylinderinchannel.ini @@ -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