diff --git a/src/basis/basis.f90 b/src/basis/basis.f90 index af66153..b006cde 100644 --- a/src/basis/basis.f90 +++ b/src/basis/basis.f90 @@ -105,7 +105,7 @@ SUBROUTINE fillBasisMapping() USE MOD_Basis_Vars,ONLY:QuadMap,QuadMapInv,VisuQuadMap,VisuQuadMapInv,Vdm_visu_Quad,D_visu_Quad USE MOD_Basis_Vars,ONLY:TetraMap,TetraMapInv !,Vdm_visu_Tetra,D_visu_Tetra USE MOD_Basis_Vars,ONLY:PyraMap,PyraMapInv !,Vdm_visu_Pyra,D_visu_Pyra -USE MOD_Basis_Vars,ONLY:PrismMap,PrismMapInv !,Vdm_visu_Prism,D_visu_Prism +USE MOD_Basis_Vars,ONLY:PrismMap,PrismMap2,PrismMapInv,Vdm_visu_Prism,D_visu_Prism,VisuPrismMap,VisuPrismMapInv USE MOD_Basis_Vars,ONLY:HexaMap,HexaMapInv,Vdm_visu_Hexa,D_visu_Hexa,VisuHexaMap,VisuHexaMapInv USE MOD_Basis_Vars,ONLY:MapSideToVol USE MOD_Basis_Vars,ONLY:EdgeToTria,EdgeToQuad @@ -132,7 +132,8 @@ SUBROUTINE fillBasisMapping() CALL getBasisMappingQuad(N,nNodes,QuadMap,QuadMapInv) CALL getBasisMappingTetra(N,nNodes,TetraMap,TetraMapInv) CALL getBasisMappingPyra(N,nNodes,PyraMap,PyraMapInv) -CALL getBasisMappingPrism(N,nNodes,PrismMap,PrismMapInv) +CALL getBasisMappingPrism(N,nNodes,PrismMap,PrismMap2,PrismMapInv) +CALL getBasisMappingPrism(nVisu,nNodes,VisuPrismMap,bMapInv=VisuPrismMapInv) ! basis mappings used for visualization (often number of visu points != order) + mapping tria->quad(tensorproduct) required CALL getBasisMappingQuad(nVisu,nNodes,VisuTriaMap,VisuTriaMapInv) @@ -144,7 +145,8 @@ SUBROUTINE fillBasisMapping() CALL getTriaBasis(N,nVisu+1,Vdm_visu_Tria,D_visu_Tria) CALL getQuadBasis(N,nVisu+1,Vdm_visu_Quad,D_visu_Quad) !CALL getTetraBasis(N,nVisu+1,Vdm_visu_Tetra,D_visu_Tetra) -CALL getHexaBasis(N,nVisu+1,Vdm_visu_Hexa,D_visu_Hexa) +CALL getPrismBasis(N,nVisu+1,Vdm_Visu_Prism,D_Visu_Prism) +CALL getHexaBasis(N,nVisu+1,Vdm_Visu_Hexa,D_Visu_Hexa) ! map edges of triangle surface counterclockwise points to BoundaBasisMappingInv(i,j) ALLOCATE(EdgeToTria(3,N+1)) @@ -220,7 +222,7 @@ SUBROUTINE fillBasisMapping() END DO !q DO q=0,N DO p=0,N-q - MapSideToVol(triaMapInv(p,q),4,6)=PyraMapInv(p,q,N) + MapSideToVol(triaMapInv(p,q),4,6)=PrismMapInv(p,q,N) MapSideToVol(triaMapInv(p,q),5,6)=PrismMapInv(q,p,0) END DO !p END DO !q diff --git a/src/basis/basis_vars.f90 b/src/basis/basis_vars.f90 index 774ce03..4a84adb 100644 --- a/src/basis/basis_vars.f90 +++ b/src/basis/basis_vars.f90 @@ -12,7 +12,7 @@ ! Copyright (C) 2017 Claus-Dieter Munz ! This file is part of HOPR, a software for the generation of high-order meshes. ! -! HOPR is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License +! HOPR is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License ! as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. ! ! HOPR is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty @@ -29,44 +29,47 @@ MODULE MOD_Basis_Vars IMPLICIT NONE PUBLIC !----------------------------------------------------------------------------------------------------------------------------------- -! GLOBAL VARIABLES +! GLOBAL VARIABLES !----------------------------------------------------------------------------------------------------------------------------------- ! Vandermande and D matrices -REAL,ALLOCATABLE,TARGET :: VdM_visu_Tria(:,:) -REAL,ALLOCATABLE,TARGET :: D_visu_Tria(:,:,:) -REAL,ALLOCATABLE,TARGET :: VdM_visu_Quad(:,:) -REAL,ALLOCATABLE,TARGET :: D_visu_Quad(:,:,:) -REAL,ALLOCATABLE,TARGET :: VdM_visu_Tetra(:,:) -REAL,ALLOCATABLE,TARGET :: D_visu_Tetra(:,:,:) -REAL,ALLOCATABLE,TARGET :: VdM_visu_Pyra(:,:) -REAL,ALLOCATABLE,TARGET :: D_visu_Pyra(:,:,:) -REAL,ALLOCATABLE,TARGET :: VdM_visu_Prism(:,:) -REAL,ALLOCATABLE,TARGET :: D_visu_Prism(:,:,:) -REAL,ALLOCATABLE,TARGET :: VdM_visu_Hexa(:,:) -REAL,ALLOCATABLE,TARGET :: D_visu_Hexa(:,:,:) +REAL,ALLOCATABLE,TARGET :: VdM_visu_Tria(:,:) +REAL,ALLOCATABLE,TARGET :: D_visu_Tria(:,:,:) +REAL,ALLOCATABLE,TARGET :: VdM_visu_Quad(:,:) +REAL,ALLOCATABLE,TARGET :: D_visu_Quad(:,:,:) +REAL,ALLOCATABLE,TARGET :: VdM_visu_Tetra(:,:) +REAL,ALLOCATABLE,TARGET :: D_visu_Tetra(:,:,:) +REAL,ALLOCATABLE,TARGET :: VdM_visu_Pyra(:,:) +REAL,ALLOCATABLE,TARGET :: D_visu_Pyra(:,:,:) +REAL,ALLOCATABLE,TARGET :: VdM_visu_Prism(:,:) +REAL,ALLOCATABLE,TARGET :: D_visu_Prism(:,:,:) +REAL,ALLOCATABLE,TARGET :: VdM_visu_Hexa(:,:) +REAL,ALLOCATABLE,TARGET :: D_visu_Hexa(:,:,:) ! Tensorproduct mappings + inverse mappings for all elements -INTEGER,ALLOCATABLE,TARGET :: TriaMap(:,:) +INTEGER,ALLOCATABLE,TARGET :: TriaMap(:,:) INTEGER,ALLOCATABLE,TARGET :: TriaMapInv(:,:) -INTEGER,ALLOCATABLE,TARGET :: QuadMap(:,:) +INTEGER,ALLOCATABLE,TARGET :: QuadMap(:,:) INTEGER,ALLOCATABLE,TARGET :: QuadMapInv(:,:) -INTEGER,ALLOCATABLE,TARGET :: TetraMap(:,:) +INTEGER,ALLOCATABLE,TARGET :: TetraMap(:,:) INTEGER,ALLOCATABLE,TARGET :: TetraMapInv(:,:,:) -INTEGER,ALLOCATABLE,TARGET :: PyraMap(:,:) +INTEGER,ALLOCATABLE,TARGET :: PyraMap(:,:) INTEGER,ALLOCATABLE,TARGET :: PyraMapInv(:,:,:) -INTEGER,ALLOCATABLE,TARGET :: PrismMap(:,:) +INTEGER,ALLOCATABLE,TARGET :: PrismMap(:,:) +INTEGER,ALLOCATABLE,TARGET :: PrismMap2(:,:) INTEGER,ALLOCATABLE,TARGET :: PrismMapInv(:,:,:) -INTEGER,ALLOCATABLE,TARGET :: HexaMap(:,:) +INTEGER,ALLOCATABLE,TARGET :: HexaMap(:,:) INTEGER,ALLOCATABLE,TARGET :: HexaMapInv(:,:,:) INTEGER :: nVisu ! number of 1D Points -1 for visualization INTEGER :: nAnalyze ! number of 1D Points -1 for analysis (Jacobian...) -INTEGER,ALLOCATABLE :: VisuTriaMap(:,:) +INTEGER,ALLOCATABLE :: VisuTriaMap(:,:) INTEGER,ALLOCATABLE :: VisuTriaMapInv(:,:) -INTEGER,ALLOCATABLE :: VisuQuadMap(:,:) +INTEGER,ALLOCATABLE :: VisuQuadMap(:,:) INTEGER,ALLOCATABLE :: VisuQuadMapInv(:,:) -INTEGER,ALLOCATABLE :: VisuHexaMap(:,:) +INTEGER,ALLOCATABLE :: VisuPrismMap(:,:) +INTEGER,ALLOCATABLE :: VisuPrismMapInv(:,:,:) +INTEGER,ALLOCATABLE :: VisuHexaMap(:,:) INTEGER,ALLOCATABLE :: VisuHexaMapInv(:,:,:) INTEGER,ALLOCATABLE :: edgeToTria(:,:) ! mapping from edges of a triangle to surface diff --git a/src/basis/prismBasis.f90 b/src/basis/prismBasis.f90 index 7986635..e5ac7a0 100644 --- a/src/basis/prismBasis.f90 +++ b/src/basis/prismBasis.f90 @@ -12,7 +12,7 @@ ! Copyright (C) 2017 Claus-Dieter Munz ! This file is part of HOPR, a software for the generation of high-order meshes. ! -! HOPR is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License +! HOPR is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License ! as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. ! ! HOPR is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty @@ -30,7 +30,7 @@ MODULE MOD_PrismBasis IMPLICIT NONE PRIVATE !----------------------------------------------------------------------------------------------------------------------------------- -! GLOBAL VARIABLES +! GLOBAL VARIABLES !----------------------------------------------------------------------------------------------------------------------------------- ! Private Part --------------------------------------------------------------------------------------------------------------------- ! Public Part ---------------------------------------------------------------------------------------------------------------------- @@ -46,9 +46,9 @@ MODULE MOD_PrismBasis !=================================================================================================================================== CONTAINS -SUBROUTINE getPrismBasis(Deg,nNodes1D) !,Vdm,GradVdm) +SUBROUTINE getPrismBasis(Deg,nNodes1D,Vdm,GradVdm) !=================================================================================================================================== -! given the degree of the orthogonal basis and the number of 1D nodes, equidistant nodes in the tetrahedron are generated and +! given the degree of the orthogonal basis and the number of 1D nodes, equidistant nodes in the tetrahedron are generated and ! we compute the Vadndermonde matrix and the Vandermondematrix of the gradient of the basis function !=================================================================================================================================== ! MODULES @@ -56,18 +56,18 @@ SUBROUTINE getPrismBasis(Deg,nNodes1D) !,Vdm,GradVdm) IMPLICIT NONE !----------------------------------------------------------------------------------------------------------------------------------- ! INPUT VARIABLES -INTEGER, INTENT(IN) :: Deg ! ? -INTEGER, INTENT(IN) :: nNodes1D ! ? +INTEGER, INTENT(IN) :: Deg ! ? +INTEGER, INTENT(IN) :: nNodes1D ! ? +REAL,ALLOCATABLE,INTENT(OUT) :: Vdm(:,:),GradVdm(:,:,:) ! ? !----------------------------------------------------------------------------------------------------------------------------------- ! OUTPUT VARIABLES -REAL,ALLOCATABLE :: Vdm(:,:),GradVdm(:,:,:) ! ? !----------------------------------------------------------------------------------------------------------------------------------- -! LOCAL VARIABLES +! LOCAL VARIABLES REAL,ALLOCATABLE :: r(:),s(:),t(:) ! coordinates in tetrahedra reference space r,s,t in [-1,1] REAL :: r1D(0:nNodes1D-1) ! equidistant 1D Lobatto nodes in [-1,1] INTEGER,ALLOCATABLE :: bMap(:,:) ! basis mapping iAns=>i,j,k INTEGER,ALLOCATABLE :: nodeMap(:,:) ! mapping for equidistant nodes iNode=>i,j,k -INTEGER :: nNodes ! ? +INTEGER :: nNodes INTEGER :: nAns ! ? INTEGER :: i,j,iNode ! ? !=================================================================================================================================== @@ -78,7 +78,7 @@ SUBROUTINE getPrismBasis(Deg,nNodes1D) !,Vdm,GradVdm) WRITE(*,*)'============================================' WRITE(*,*)'getPrismBasis for Degree', Deg WRITE(*,*)'============================================' -!equidistant nodes +!equidistant nodes r1D=0. DO i=0,nNodes1D-1 r1D(i)=-1.+2.*REAL(i)/REAL(nNodes1D-1) @@ -135,27 +135,28 @@ SUBROUTINE getPrismBasis(Deg,nNodes1D) !,Vdm,GradVdm) END DO WRITE(*,*)' ' END DO -DEALLOCATE(bMap,nodeMap) +DEALLOCATE(bMap,nodeMap) END SUBROUTINE getPrismBasis -SUBROUTINE getBasisMappingPrism(Deg,nAns,bMap,bMapInv) +SUBROUTINE getBasisMappingPrism(Deg,nAns,bMap,bMap2,bMapInv) !=================================================================================================================================== -! mapping from iAns -> i,j in [0,Deg], can be used for nodeMap too: CALL getBasisMapping(nNodes1D-1,nodeMap) +! mapping from iAns -> i,j in [0,Deg], can be used for nodeMap too: CALL getBasisMapping(nNodes1D-1,nodeMap) !=================================================================================================================================== ! MODULES ! IMPLICIT VARIABLE HANDLING IMPLICIT NONE !----------------------------------------------------------------------------------------------------------------------------------- ! INPUT VARIABLES -INTEGER, INTENT(IN) :: Deg ! ? +INTEGER, INTENT(IN) :: Deg ! ? !----------------------------------------------------------------------------------------------------------------------------------- ! OUTPUT VARIABLES -INTEGER, INTENT(OUT) :: nAns ! ? +INTEGER, INTENT(OUT) :: nAns ! ? INTEGER,ALLOCATABLE,INTENT(OUT) :: bMap(:,:) ! ? +INTEGER,ALLOCATABLE,OPTIONAL,INTENT(OUT) :: bMap2(:,:) ! ? INTEGER,ALLOCATABLE,OPTIONAL,INTENT(OUT) :: bMapInv(:,:,:) ! ? !----------------------------------------------------------------------------------------------------------------------------------- -! LOCAL VARIABLES +! LOCAL VARIABLES INTEGER :: iAns,i,j,k ! ? !=================================================================================================================================== nAns=(Deg+1)*(Deg+1)*(Deg+2)/2 @@ -164,11 +165,26 @@ SUBROUTINE getBasisMappingPrism(Deg,nAns,bMap,bMapInv) DO k=0,Deg DO j=0,Deg DO i=0,Deg-j + !DO i=Deg-j,0,-1 iAns=iAns+1 - bMap(iAns,:)=(/i,j,k/) + bMap(iAns,:)=(/k,i,j/) + ! bMap(iAns,:)=(/i,j,k/) END DO END DO END DO +IF(PRESENT(bMap2))THEN + ALLOCATE(bMap2(nAns,3)) + iAns=0 + DO k=0,Deg + DO j=Deg,0,-1 + DO i=Deg-j,Deg + !DO i=Deg,Deg-j,-1 + iAns=iAns+1 + bMap2(iAns,:)=(/k,j,i/) + END DO + END DO + END DO +END IF IF(PRESENT(bMapInv))THEN ALLOCATE(bMapInv(0:Deg,0:Deg,0:Deg)) bMapInv=0 @@ -201,7 +217,7 @@ SUBROUTINE rst2abcPrism(nNodes,r,s,t,a,b,c) ! OUTPUT VARIABLES REAL,INTENT(OUT) :: a(nNodes),b(nNodes),c(nNodes) ! ? !----------------------------------------------------------------------------------------------------------------------------------- -! LOCAL VARIABLES +! LOCAL VARIABLES INTEGER :: iNode ! ? !=================================================================================================================================== WRITE(*,*)'entering rst2abcPrism' @@ -213,14 +229,14 @@ SUBROUTINE rst2abcPrism(nNodes,r,s,t,a,b,c) END IF END DO !iNode b=s -c=t +c=t END SUBROUTINE rst2abcPrism SUBROUTINE VandermondePrism(nNodes,nAns,bMap,r,s,t,VdM) !=================================================================================================================================== -! For a given vector of nNodes in reference coordinates (r,s,t) of the prism, and nAns=(Deg+1)*(Deg+1)*(Deg+2)/2 -! basis functions, we compute the 3D Vandermonde matrix. +! For a given vector of nNodes in reference coordinates (r,s,t) of the prism, and nAns=(Deg+1)*(Deg+1)*(Deg+2)/2 +! basis functions, we compute the 3D Vandermonde matrix. ! The polynomial basis is orthonormal with respect to the reference prism r,s,t>=-1, r+s<=1 , t<=1 !=================================================================================================================================== ! MODULES @@ -236,14 +252,14 @@ SUBROUTINE VandermondePrism(nNodes,nAns,bMap,r,s,t,VdM) ! OUTPUT VARIABLES REAL,INTENT(OUT) :: VdM(nNodes,nAns) ! ? !----------------------------------------------------------------------------------------------------------------------------------- -! LOCAL VARIABLES +! LOCAL VARIABLES REAL,DIMENSION(nNodes) :: a,b,c,Pa,Pb,Pc ! ? INTEGER :: iAns,i,j,k ! ? REAL :: norm ! ? !=================================================================================================================================== norm=SQRT(2.) WRITE(*,*)'entering VandermondePrism' -CALL rst2abcPrism(nNodes,r,s,t,a,b,c) +CALL rst2abcPrism(nNodes,r,s,t,a,b,c) DO iAns=1,nAns i=bMap(iAns,1) j=bMap(iAns,2) @@ -254,14 +270,14 @@ SUBROUTINE VandermondePrism(nNodes,nAns,bMap,r,s,t,VdM) CALL JacobiP(nNodes,c, 0, 0,k,Pc) VdM(:,iAns)=norm*(Pa*Pb*Pc) !orthonormal IF(i.GT.0) VdM(:,iAns)=VdM(:,iAns)*((1.-b)**i) -END DO +END DO END SUBROUTINE VandermondePrism SUBROUTINE GradVandermondePrism(nNodes,nAns,bMap,r,s,t,gradVdM) !=================================================================================================================================== -! For a given vector of nNodes in reference coordinates (r,s,t) of the prism, and nAns=(Deg+1)*(Deg+1)*(Deg+2)/2 -! basis functions, we compute the 3D Gradient Vandermonde matrix. +! For a given vector of nNodes in reference coordinates (r,s,t) of the prism, and nAns=(Deg+1)*(Deg+1)*(Deg+2)/2 +! basis functions, we compute the 3D Gradient Vandermonde matrix. ! The polynomial basis is orthonormal with respect to the reference prism r,s,t>=-1, r+s<=1 , t<=1 !=================================================================================================================================== ! MODULES @@ -277,19 +293,19 @@ SUBROUTINE GradVandermondePrism(nNodes,nAns,bMap,r,s,t,gradVdM) ! OUTPUT VARIABLES REAL,INTENT(OUT) :: gradVdM(nNodes,nAns,3) ! ? !----------------------------------------------------------------------------------------------------------------------------------- -! LOCAL VARIABLES +! LOCAL VARIABLES REAL,DIMENSION(nNodes) :: a,b,c,Pa,Pb,Pc,dPa,dPb,dPc,dPhi_da,dPhi_db,dPhi_dc ! ? INTEGER :: iAns,i,j,k ! ? REAL :: norm ! ? !=================================================================================================================================== norm=SQRT(2.) WRITE(*,*)'entering GradVandermondePrism' -CALL rst2abcPrism(nNodes,r,s,t,a,b,c) +CALL rst2abcPrism(nNodes,r,s,t,a,b,c) DO iAns=1,nAns i=bMap(iAns,1) j=bMap(iAns,2) k=bMap(iAns,3) - + CALL JacobiP(nNodes,a, 0, 0,i,Pa) CALL JacobiP(nNodes,b, 2*i+1, 0,j,Pb) CALL JacobiP(nNodes,c, 0, 0,k,Pc) @@ -297,14 +313,14 @@ SUBROUTINE GradVandermondePrism(nNodes,nAns,bMap,r,s,t,gradVdM) CALL GradJacobiP(nNodes,b, 2*i+1, 0,j,dPb) CALL GradJacobiP(nNodes,c, 0, 0,k,dPc) !dphi/da*(1-b)^-1 - dPhi_da=norm*(dPa*Pb*Pc) + dPhi_da=norm*(dPa*Pb*Pc) IF(i.GT.1) dPhi_da=dPhi_da*((1.-b)**(i-1)) !dphi/db dPhi_db=dPb IF(i.GT.0) dPhi_db=dPhi_db*(1.-b)**i IF(i.EQ.1) dPhi_db=dPhi_db-Pb IF(i.GT.1) dPhi_db=dPhi_db-i*(1.-b)**(i-1)*Pb - dPhi_db=norm*Pa*Pc*dPhi_db + dPhi_db=norm*Pa*Pc*dPhi_db !dphi/dc dPhi_dc=dPc dPhi_dc=norm*Pa*Pb*dPhi_dc @@ -315,7 +331,7 @@ SUBROUTINE GradVandermondePrism(nNodes,nAns,bMap,r,s,t,gradVdM) gradVdM(:,iAns,2)=(1.+a)*dPhi_da+dPhi_db ! t-derivative gradVdM(:,iAns,3)=dPhi_dc -END DO +END DO END SUBROUTINE GradVandermondePrism END MODULE MOD_PrismBasis diff --git a/src/basis/pyraBasis.f90 b/src/basis/pyraBasis.f90 index db69df4..2a3294f 100644 --- a/src/basis/pyraBasis.f90 +++ b/src/basis/pyraBasis.f90 @@ -12,7 +12,7 @@ ! Copyright (C) 2017 Claus-Dieter Munz ! This file is part of HOPR, a software for the generation of high-order meshes. ! -! HOPR is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License +! HOPR is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License ! as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. ! ! HOPR is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty @@ -30,7 +30,7 @@ MODULE MOD_pyraBasis IMPLICIT NONE PRIVATE !----------------------------------------------------------------------------------------------------------------------------------- -! GLOBAL VARIABLES +! GLOBAL VARIABLES !----------------------------------------------------------------------------------------------------------------------------------- ! Private Part --------------------------------------------------------------------------------------------------------------------- ! Public Part ---------------------------------------------------------------------------------------------------------------------- @@ -48,7 +48,7 @@ MODULE MOD_pyraBasis CONTAINS SUBROUTINE getPyraBasis(Deg,nNodes1D) !,Vdm,GradVdm) !=================================================================================================================================== -! given the degree of the orthogonal basis and the number of 1D nodes, equidistant nodes in the tetrahedron are generated and +! given the degree of the orthogonal basis and the number of 1D nodes, equidistant nodes in the tetrahedron are generated and ! we compute the Vadndermonde matrix and the Vandermondematrix of the gradient of the basis function !=================================================================================================================================== ! MODULES @@ -62,7 +62,7 @@ SUBROUTINE getPyraBasis(Deg,nNodes1D) !,Vdm,GradVdm) ! OUTPUT VARIABLES REAL,ALLOCATABLE :: Vdm(:,:),GradVdm(:,:,:) ! ? !----------------------------------------------------------------------------------------------------------------------------------- -! LOCAL VARIABLES +! LOCAL VARIABLES REAL,ALLOCATABLE :: r(:),s(:),t(:) ! coordinates in tetrahedra reference space r,s,t in [-1,1] REAL :: r1D(0:nNodes1D-1) ! equidistant 1D Lobatto nodes in [-1,1] INTEGER,ALLOCATABLE :: bMap(:,:) ! basis mapping iAns=>i,j,k @@ -78,7 +78,7 @@ SUBROUTINE getPyraBasis(Deg,nNodes1D) !,Vdm,GradVdm) WRITE(*,*)'============================================' WRITE(*,*)'getPyraBasis for Degree', Deg WRITE(*,*)'============================================' -!equidistant nodes +!equidistant nodes r1D=0. DO i=0,nNodes1D-1 r1D(i)=-1.+2.*REAL(i)/REAL(nNodes1D-1) @@ -135,13 +135,13 @@ SUBROUTINE getPyraBasis(Deg,nNodes1D) !,Vdm,GradVdm) END DO WRITE(*,*)' ' END DO -DEALLOCATE(bMap,nodeMap) +DEALLOCATE(bMap,nodeMap) END SUBROUTINE getPyraBasis SUBROUTINE getBasisMappingPyra(Deg,nAns,bMap,bMapInv) !=================================================================================================================================== -! mapping from iAns -> i,j in [0,Deg], can be used for nodeMap too: CALL getBasisMapping(nNodes1D-1,nodeMap) +! mapping from iAns -> i,j in [0,Deg], can be used for nodeMap too: CALL getBasisMapping(nNodes1D-1,nodeMap) !=================================================================================================================================== ! MODULES ! IMPLICIT VARIABLE HANDLING @@ -155,7 +155,7 @@ SUBROUTINE getBasisMappingPyra(Deg,nAns,bMap,bMapInv) INTEGER,ALLOCATABLE,INTENT(OUT) :: bMap(:,:) ! ? INTEGER,ALLOCATABLE,OPTIONAL,INTENT(OUT):: bMapInv(:,:,:) ! ? !----------------------------------------------------------------------------------------------------------------------------------- -! LOCAL VARIABLES +! LOCAL VARIABLES INTEGER :: iAns,i,j,k ! ? !=================================================================================================================================== nAns=(Deg+1)*(Deg+2)*(2*Deg+3)/6 @@ -165,7 +165,7 @@ SUBROUTINE getBasisMappingPyra(Deg,nAns,bMap,bMapInv) DO j=0,Deg-k DO i=0,Deg-k iAns=iAns+1 - bMap(iAns,:)=(/i,j,k/) + bMap(iAns,:)=(/k,j,i/) END DO END DO END DO @@ -201,7 +201,7 @@ SUBROUTINE rst2abcPyra(nNodes,r,s,t,a,b,c) ! OUTPUT VARIABLES REAL,INTENT(OUT) :: a(nNodes),b(nNodes),c(nNodes) ! ? !----------------------------------------------------------------------------------------------------------------------------------- -! LOCAL VARIABLES +! LOCAL VARIABLES INTEGER :: iNode ! ? !=================================================================================================================================== WRITE(*,*)'entering rst_to_abc' @@ -214,14 +214,14 @@ SUBROUTINE rst2abcPyra(nNodes,r,s,t,a,b,c) b(iNode) = 2.*(1.+s(iNode))/(1.-t(iNode))-1. END IF END DO !iNode -c=t +c=t END SUBROUTINE rst2abcPyra SUBROUTINE VandermondePyra(nNodes,nAns,bMap,r,s,t,VdM) !=================================================================================================================================== -! For a given vector of nNodes in reference coordinates (r,s,t) of the pyramid, and nAns=(Deg+1)*(Deg+2)*(2*Deg+3)/6 -! basis functions, we compute the 3D Vandermonde matrix. +! For a given vector of nNodes in reference coordinates (r,s,t) of the pyramid, and nAns=(Deg+1)*(Deg+2)*(2*Deg+3)/6 +! basis functions, we compute the 3D Vandermonde matrix. ! The polynomial basis is orthonormal with respect to the reference pyramid r,s,t>=-1, r+t<=1 ,s+t<=1 !=================================================================================================================================== ! MODULES @@ -237,14 +237,14 @@ SUBROUTINE VandermondePyra(nNodes,nAns,bMap,r,s,t,VdM) ! OUTPUT VARIABLES REAL,INTENT(OUT) :: VdM(nNodes,nAns) ! ? !----------------------------------------------------------------------------------------------------------------------------------- -! LOCAL VARIABLES +! LOCAL VARIABLES REAL,DIMENSION(nNodes) :: a,b,c,Pa,Pb,Pc ! ? INTEGER :: iAns,i,j,k,ij ! ? REAL :: norm ! ? !=================================================================================================================================== norm=2. WRITE(*,*)'entering VandermondePyra' -CALL rst2abcPyra(nNodes,r,s,t,a,b,c) +CALL rst2abcPyra(nNodes,r,s,t,a,b,c) DO iAns=1,nAns i=bMap(iAns,1) j=bMap(iAns,2) @@ -259,14 +259,14 @@ SUBROUTINE VandermondePyra(nNodes,nAns,bMap,r,s,t,VdM) CALL JacobiP(nNodes,c, 2*ij+2, 0,k,Pc) VdM(:,iAns)=norm*(Pa*Pb*Pc) !orthonormal IF(ij.GT.0) VdM(:,iAns)=VdM(:,iAns)*((1.-c)**(ij)) -END DO +END DO END SUBROUTINE VandermondePyra SUBROUTINE GradVandermondePyra(nNodes,nAns,bMap,r,s,t,gradVdM) !=================================================================================================================================== -! For a given vector of nNodes in reference coordinates (r,s,t) of the pyramid, and nAns=(Deg+1)*(Deg+2)*(2*Deg+3)/6 -! basis functions, we compute the 3D Gradient Vandermonde matrix. +! For a given vector of nNodes in reference coordinates (r,s,t) of the pyramid, and nAns=(Deg+1)*(Deg+2)*(2*Deg+3)/6 +! basis functions, we compute the 3D Gradient Vandermonde matrix. ! The polynomial basis is orthonormal with respect to the reference pyramid r,s,t>=-1, r+t<=1 ,s+t<=1 !=================================================================================================================================== ! MODULES @@ -282,14 +282,14 @@ SUBROUTINE GradVandermondePyra(nNodes,nAns,bMap,r,s,t,gradVdM) ! OUTPUT VARIABLES REAL,INTENT(OUT) :: gradVdM(nNodes,nAns,3) ! ? !----------------------------------------------------------------------------------------------------------------------------------- -! LOCAL VARIABLES +! LOCAL VARIABLES REAL,DIMENSION(nNodes) :: a,b,c,Pa,Pb,Pc,dPa,dPb,dPc,dPhi_da,dPhi_db,dPhi_dc ! ? INTEGER :: iAns,i,j,k,ij ! ? REAL :: norm ! ? !=================================================================================================================================== norm=2. WRITE(*,*)'entering GradVandermondePyra' -CALL rst2abcPyra(nNodes,r,s,t,a,b,c) +CALL rst2abcPyra(nNodes,r,s,t,a,b,c) DO iAns=1,nAns i=bMap(iAns,1) j=bMap(iAns,2) @@ -298,7 +298,7 @@ SUBROUTINE GradVandermondePyra(nNodes,nAns,bMap,r,s,t,gradVdM) ! ij=max(i,j) ! karniadakis ij=(i+j) - + CALL JacobiP(nNodes,a, 0, 0,i,Pa) CALL JacobiP(nNodes,b, 0, 0,j,Pb) CALL JacobiP(nNodes,c, 2*ij+2, 0,k,Pc) @@ -306,7 +306,7 @@ SUBROUTINE GradVandermondePyra(nNodes,nAns,bMap,r,s,t,gradVdM) CALL GradJacobiP(nNodes,b, 0, 0,j,dPb) CALL GradJacobiP(nNodes,c, 2*ij+2, 0,k,dPc) !dphi/da*(1-c)^-1 - dPhi_da=norm*(dPa*Pb*Pc) + dPhi_da=norm*(dPa*Pb*Pc) IF(ij.GT.1) dPhi_da=dPhi_da*((1.-c)**(ij-1)) !dphi/db*(1-c)^-1 dPhi_db=norm*Pa*dPb*Pc @@ -323,7 +323,7 @@ SUBROUTINE GradVandermondePyra(nNodes,nAns,bMap,r,s,t,gradVdM) gradVdM(:,iAns,2)=2.*dPhi_db ! t-derivative gradVdM(:,iAns,3)=(1.+a)*dPhi_da+(1.+b)*dPhi_db+dPhi_dc -END DO +END DO END SUBROUTINE GradVandermondePyra END MODULE MOD_pyraBasis diff --git a/src/mesh/cartmesh.f90 b/src/mesh/cartmesh.f90 index 96fa658..e63822b 100644 --- a/src/mesh/cartmesh.f90 +++ b/src/mesh/cartmesh.f90 @@ -12,7 +12,7 @@ ! Copyright (C) 2017 Claus-Dieter Munz ! This file is part of HOPR, a software for the generation of high-order meshes. ! -! HOPR is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License +! HOPR is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License ! as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. ! ! HOPR is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty @@ -31,7 +31,7 @@ MODULE MOD_CartMesh IMPLICIT NONE PRIVATE !----------------------------------------------------------------------------------------------------------------------------------- -! GLOBAL VARIABLES +! GLOBAL VARIABLES !----------------------------------------------------------------------------------------------------------------------------------- ! Private Part --------------------------------------------------------------------------------------------------------------------- TYPE tCartesianMesh ! provides data structure for Cartesian mesh @@ -395,22 +395,29 @@ SUBROUTINE GetNewPrism(CornerNode) aElem%nNodes=6 ALLOCATE(aElem%Node(aElem%nNodes)) aElem%Node(1)%NP=>CornerNode(1)%NP -aElem%Node(2)%NP=>CornerNode(4)%NP -aElem%Node(3)%NP=>CornerNode(5)%NP -aElem%Node(4)%NP=>CornerNode(2)%NP -aElem%Node(5)%NP=>CornerNode(3)%NP -aElem%Node(6)%NP=>CornerNode(6)%NP +aElem%Node(2)%NP=>CornerNode(2)%NP +aElem%Node(3)%NP=>CornerNode(4)%NP +aElem%Node(4)%NP=>CornerNode(5)%NP +aElem%Node(5)%NP=>CornerNode(6)%NP +aElem%Node(6)%NP=>CornerNode(8)%NP + +! NodeCoords(1:3,1) = (/0.,0.,0./) +! NodeCoords(1:3,2) = (/1.,0.,0./) +! NodeCoords(1:3,3) = (/0.,1.,0./) +! NodeCoords(1:3,4) = (/0.,0.,1./) +! NodeCoords(1:3,5) = (/1.,0.,1./) +! NodeCoords(1:3,6) = (/0.,1.,1./) CALL getNewElem(aElem%nextElem) aElem%nextElem%prevElem => aElem aElem => aElem%nextElem aElem%nNodes=6 ALLOCATE(aElem%Node(aElem%nNodes)) -aElem%Node(1)%NP=>CornerNode(4)%NP -aElem%Node(2)%NP=>CornerNode(8)%NP -aElem%Node(3)%NP=>CornerNode(5)%NP -aElem%Node(4)%NP=>CornerNode(3)%NP +aElem%Node(1)%NP=>CornerNode(2)%NP +aElem%Node(2)%NP=>CornerNode(3)%NP +aElem%Node(3)%NP=>CornerNode(4)%NP +aElem%Node(4)%NP=>CornerNode(6)%NP aElem%Node(5)%NP=>CornerNode(7)%NP -aElem%Node(6)%NP=>CornerNode(6)%NP +aElem%Node(6)%NP=>CornerNode(8)%NP ! Add elements to list IF(.NOT.ASSOCIATED(FirstElem))THEN @@ -482,7 +489,7 @@ SUBROUTINE CartesianMesh() USE MOD_Mesh_Vars,ONLY:deleteSide,deleteNode USE MOD_Mesh_Vars,ONLY:InnerElemStretch,BoundaryOrder USE MOD_Mesh_Basis,ONLY:CreateSides -USE MOD_CurvedCartMesh,ONLY:GetNewCurvedHexahedron +USE MOD_CurvedCartMesh,ONLY:GetNewCurvedHexahedron,GetNewCurvedPrism,GetNewCurvedPyram IMPLICIT NONE !----------------------------------------------------------------------------------------------------------------------------------- ! INPUT VARIABLES @@ -497,7 +504,7 @@ SUBROUTINE CartesianMesh() !----------------------------------------------------------------------------------------------------------------------------------- ! OUTPUT VARIABLES !----------------------------------------------------------------------------------------------------------------------------------- -! LOCAL VARIABLES +! LOCAL VARIABLES TYPE(tCartesianMesh),POINTER :: CartMesh ! ? TYPE(tNodePtr) :: CornerNode(8) ! ? TYPE(tNodePtr),POINTER :: Mnodes(:,:,:) ! ? @@ -533,9 +540,9 @@ SUBROUTINE CartesianMesh() END SELECT IF(ABS(CartMesh%l0(i_Dim)) .LT. PP_RealTolerance ) THEN !l0 = 0 = inactive IF(ABS(CartMesh%factor(i_Dim)) .LT. PP_RealTolerance ) THEN !fac= 0 = equidistant - fac(i_Dim)=1. + fac(i_Dim)=1. ELSE ! stretched, (nElem,f) given - fac(i_Dim)=(ABS(CartMesh%factor(i_Dim)))**(SIGN(1.,CartMesh%factor(i_Dim))) ! sign for direction + fac(i_Dim)=(ABS(CartMesh%factor(i_Dim)))**(SIGN(1.,CartMesh%factor(i_Dim))) ! sign for direction END IF ELSE !l0 active dx(i_Dim)=e1(i_Dim)/ABS(cartMesh%l0(i_Dim)) ! l/l0 @@ -543,22 +550,22 @@ SUBROUTINE CartesianMesh() CALL abort(__STAMP__,'Stretching error, length l0 longer than grid region, in direction') END IF IF(ABS(CartMesh%factor(i_Dim)) .LT. PP_RealTolerance ) THEN ! fac=0 , (nElem,l0) given, fac calculated - ! + ! ELSE ! (factor, l0) given, change nElems - fac(i_Dim)=(ABS(CartMesh%factor(i_Dim)))**(SIGN(1.,CartMesh%factor(i_Dim)*CartMesh%l0(i_Dim))) ! sign for direction + fac(i_Dim)=(ABS(CartMesh%factor(i_Dim)))**(SIGN(1.,CartMesh%factor(i_Dim)*CartMesh%l0(i_Dim))) ! sign for direction IF(fac(i_Dim) .NE. 1.) THEN - CartMesh%nElems(i_Dim)=NINT(LOG(1.-dx(i_Dim)*(1.-fac(i_Dim))) / LOG(fac(i_Dim))) !nearest Integer + CartMesh%nElems(i_Dim)=NINT(LOG(1.-dx(i_Dim)*(1.-fac(i_Dim))) / LOG(fac(i_Dim))) !nearest Integer ELSE CartMesh%nElems(i_Dim)=NINT(dx(i_Dim)) END IF IF (CartMesh%nElems(i_Dim) .LT. 1) CartMesh%nElems(i_Dim)=1 WRITE(UNIT_stdOut,*)' -Element number in dir',i_Dim,'changed from', & - nElems(i_Dim),'to',CartMesh%nElems(i_Dim) - nElems(i_Dim)=CartMesh%nElems(i_Dim) + nElems(i_Dim),'to',CartMesh%nElems(i_Dim) + nElems(i_Dim)=CartMesh%nElems(i_Dim) END IF - IF(nElems(i_Dim).EQ. 1) THEN - fac(i_Dim)=1.0 - ELSEIF(nElems(i_Dim).EQ. 2) THEN + IF(nElems(i_Dim).EQ. 1) THEN + fac(i_Dim)=1.0 + ELSEIF(nElems(i_Dim).EQ. 2) THEN fac(i_Dim)=dx(i_Dim)-1. ELSE !nElems > 2 fac(i_Dim)=dx(i_Dim)/nElems(i_Dim) !start value for Newton iteration @@ -576,7 +583,7 @@ SUBROUTINE CartesianMesh() fac(i_Dim)=fac(i_Dim)**SIGN(1.,CartMesh%l0(i_Dim)) ! sign for direction END IF END IF - WRITE(UNIT_stdOut,*)' -stretching factor in dir',i_Dim,'is now', fac(i_Dim) + WRITE(UNIT_stdOut,*)' -stretching factor in dir',i_Dim,'is now', fac(i_Dim) END IF IF( ABS((nElems(i_Dim)-1.)*LOG(fac(i_Dim))/LOG(10.)) .GE. 4. )CALL abort(__STAMP__,'Stretching error, length ratio > 1.0E4 in direction') @@ -595,7 +602,7 @@ SUBROUTINE CartesianMesh() WRITE(UNIT_stdOut,formatstr) ' -factor(:) : ', fac(:) END IF - IF(InnerElemStretch.AND.CartMesh%ElemType.EQ.108)THEN + IF(InnerElemStretch)THEN Ngeo=BoundaryOrder-1 ALLOCATE(CurvedNode(0:Ngeo,0:Ngeo,0:Ngeo)) ELSE @@ -624,7 +631,7 @@ SUBROUTINE CartesianMesh() DO l=0,ne(1) ! x1(1)=REAL(l)/REAL(nSplit(1)) x1(2)=x2(2) - dx(2)=e1(2) + dx(2)=e1(2) DO m=0,ne(2) ! x1(2)=REAL(m)/REAL(nSplit(2)) x1(3)=x2(3) @@ -661,16 +668,30 @@ SUBROUTINE CartesianMesh() CASE(104) ! Tetrahedron CALL GetNewTetrahedron(CornerNode,CartMesh,l+1,m+1,n+1,ind=NodeInd) CASE(105) ! Pyramid - CALL GetNewPyramid(CornerNode) + IF(InnerElemStretch)THEN + DO i=0,Ngeo; DO j=0,Ngeo; DO k=0,Ngeo + CurvedNode(i,j,k)%np=>Mnodes(l+i,m+j,n+k)%np + END DO; END DO; END DO; + CALL GetNewCurvedPyram(CurvedNode,Ngeo,iZone) + ELSE + CALL GetNewPyramid(CornerNode) + END IF CASE(106) ! Dreiecks-prismon - CALL GetNewPrism(CornerNode) + IF(InnerElemStretch)THEN + DO i=0,Ngeo; DO j=0,Ngeo; DO k=0,Ngeo + CurvedNode(i,j,k)%np=>Mnodes(l+i,m+j,n+k)%np + END DO; END DO; END DO; + CALL GetNewCurvedPrism(CurvedNode,Ngeo,iZone) + ELSE + CALL GetNewPrism(CornerNode) + END IF CASE(108) ! Hexaeder IF(InnerElemStretch)THEN DO i=0,Ngeo; DO j=0,Ngeo; DO k=0,Ngeo CurvedNode(i,j,k)%np=>Mnodes(l+i,m+j,n+k)%np - END DO; END DO; END DO; + END DO; END DO; END DO; CALL GetNewCurvedHexahedron(CurvedNode,Ngeo,iZone) - ELSE + ELSE CALL GetNewHexahedron(CornerNode) END IF CASE DEFAULT @@ -683,7 +704,7 @@ SUBROUTINE CartesianMesh() DO l=0,ne(1),Ngeo DO m=0,ne(2),Ngeo DO n=0,ne(3),Ngeo - Mnodes(l,m,n)%np%tmp=0 + Mnodes(l,m,n)%np%tmp=0 IF(n.EQ.0) Mnodes(l,m,n)%np%tmp=Mnodes(l,m,n)%np%tmp+1 !zeta minus IF(m.EQ.0) Mnodes(l,m,n)%np%tmp=Mnodes(l,m,n)%np%tmp+20 !eta minus IF(l.EQ.ne(1) ) Mnodes(l,m,n)%np%tmp=Mnodes(l,m,n)%np%tmp+300 !xi plus diff --git a/src/mesh/curved.f90 b/src/mesh/curved.f90 index 1e04be7..26c43f5 100644 --- a/src/mesh/curved.f90 +++ b/src/mesh/curved.f90 @@ -13,7 +13,7 @@ ! Copyright (C) 2017 Claus-Dieter Munz ! This file is part of HOPR, a software for the generation of high-order meshes. ! -! HOPR is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License +! HOPR is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License ! as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. ! ! HOPR is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty @@ -32,7 +32,7 @@ MODULE MOD_Curved IMPLICIT NONE PRIVATE !----------------------------------------------------------------------------------------------------------------------------------- -! GLOBAL VARIABLES +! GLOBAL VARIABLES !----------------------------------------------------------------------------------------------------------------------------------- ! Private Part --------------------------------------------------------------------------------------------------------------------- ! Public Part ---------------------------------------------------------------------------------------------------------------------- @@ -176,13 +176,13 @@ SUBROUTINE readNormals() READ(150,'(I6)',ADVANCE='NO') k !first FaceID IF (k .LT. 0) THEN nPlanars=nPlanars+1 - faceConnectivity(1,nFaces+1+k)=k !actual Face ID - READ(150,*) faceConnectivity(2:j,nFaces+1+k) !Neighbor Face IDs of actual Face + faceConnectivity(1,nFaces+1+k)=k !actual Face ID + READ(150,*) faceConnectivity(2:j,nFaces+1+k) !Neighbor Face IDs of actual Face ELSEIF (k .EQ. 0) THEN CALL abort(__STAMP__,'ERROR: CAD-face connectivity invalid',999,999.) ELSE - faceConnectivity(1,k)=k - READ(150,*) faceConnectivity(2:j,k) + faceConnectivity(1,k)=k + READ(150,*) faceConnectivity(2:j,k) END IF END DO DO i=1,nFaces !all entries have to be filled @@ -211,7 +211,7 @@ SUBROUTINE readNormals() DO j=1,aSide%nNodes aNode=>aSide%node(j)%np IF(aNode%tmp .EQ. 0) THEN - aNode%tmp=999 + aNode%tmp=999 CALL insertNode(searchMesh,aNode) aNode%refCount=aNode%refCount-1 searchMeshNodes=searchMeshNodes+1 @@ -272,9 +272,9 @@ SUBROUTINE readNormals() ToObject=>getNextToObject(searchMesh,.FALSE.) END DO IF (.NOT. pointFound) THEN - IF (useProjections) READ(150,*) + IF (useProjections) READ(150,*) DO j = 1, nVects - READ(150,*) + READ(150,*) END DO END IF END DO !nTotalNodes @@ -319,7 +319,7 @@ SUBROUTINE readNormals() END IF END DO ! mark same Faces - IF (sameFaces(face(1),face(2))) CYCLE + IF (sameFaces(face(1),face(2))) CYCLE sameFaces(face(1),face(2))=.TRUE. sameFaces(face(2),face(1))=.TRUE. ! look if already another face is the same Face @@ -327,12 +327,12 @@ SUBROUTINE readNormals() IF (sameFaces(face(1),i)) THEN sameFaces(i,face(2))=.TRUE. sameFaces(face(2),i)=.TRUE. - END IF + END IF IF (sameFaces(face(2),i)) THEN sameFaces(i,face(1))=.TRUE. sameFaces(face(1),i)=.TRUE. END IF - END DO + END DO END DO END IF aSide=>aSide%nextElemSide @@ -355,7 +355,7 @@ SUBROUTINE readNormals() ! determine for a group of sameFaces the mapping from the actual faceID to the biggest FaceID in the group mergedFaces(:)=0 DO i=1,nFaces ! - DO j=(nFaces-nPlanars),1,-1 !curved Faces, backwards, so j is the biggest FaceID + DO j=(nFaces-nPlanars),1,-1 !curved Faces, backwards, so j is the biggest FaceID IF (sameFaces(i,j)) THEN mergedFaces(i)=j EXIT @@ -388,12 +388,12 @@ SUBROUTINE readNormals() DO j=1, aElem%nNodes IF (aElem%node(j)%np%tmp .NE. 0) CYCLE i=i+1 - aElem%node(j)%np%tmp=i + aElem%node(j)%np%tmp=i aNormal=>aElem%node(j)%np%firstNormal DO WHILE (ASSOCIATED(aNormal)) DO k=1, SIZE(aNormal%FaceID) IF (aNormal%FaceID(k) .LT. 0) THEN !planar face - aNormal%FaceID(k)=mergedFaces(nFaces+1+aNormal%FaceID(k)) !replace old FaceId with new one + aNormal%FaceID(k)=mergedFaces(nFaces+1+aNormal%FaceID(k)) !replace old FaceId with new one ELSE aNormal%FaceID(k)=mergedFaces(aNormal%FaceID(k)) !replace old FaceId with new one END IF @@ -417,7 +417,7 @@ SUBROUTINE readNormals() DO j=1, aElem%nNodes IF (aElem%node(j)%np%tmp .NE. 0) CYCLE i=i+1 - aElem%node(j)%np%tmp=i + aElem%node(j)%np%tmp=i !remove duplicates inside each normal aNormal=>aElem%node(j)%np%firstNormal DO WHILE (ASSOCIATED(aNormal)) @@ -443,7 +443,7 @@ SUBROUTINE readNormals() DuplicatesInd(l)=aNormal%FaceID(k) END IF END DO - DEALLOCATE(aNormal%FaceID) + DEALLOCATE(aNormal%FaceID) ALLOCATE(aNormal%FaceID(SIZE(duplicatesInd))) aNormal%FaceID=DuplicatesInd DEALLOCATE(duplicatesInd) @@ -471,7 +471,7 @@ SUBROUTINE readNormals() ALLOCATE(DuplicatesInd(SIZE(aNormal%FaceID)+SIZE(bNormal%FaceID)-Counter)) DuplicatesInd(1:SIZE(aNormal%FaceID))=aNormal%FaceID(:) l=SIZE(aNormal%FaceID) - DO k=1, SIZE(bNormal%FaceID) + DO k=1, SIZE(bNormal%FaceID) IF (bNormal%FaceID(k) .EQ. 0) THEN CYCLE ELSE @@ -479,11 +479,11 @@ SUBROUTINE readNormals() DuplicatesInd(l)=bNormal%FaceID(k) END IF END DO - DEALLOCATE(aNormal%FaceID) + DEALLOCATE(aNormal%FaceID) ALLOCATE(aNormal%FaceID(SIZE(duplicatesInd))) aNormal%FaceID=DuplicatesInd DEALLOCATE(duplicatesInd) - + !average normals, remove bNormal and reorder pointers aNormal%normal=aNormal%normal+bNormal%normal IF (ASSOCIATED(bNormal%nextNormal)) THEN @@ -506,7 +506,7 @@ SUBROUTINE readNormals() aElem=>aElem%nextElem END DO -! If each node of a side has same negative FaceID => CAD face is not curved => Side ist not curved! +! If each node of a side has same negative FaceID => CAD face is not curved => Side ist not curved! ! => set aSide%curveIndex = 0 => no spline reconstruction => better performance (in work) l=0 aElem=>firstElem @@ -518,7 +518,7 @@ SUBROUTINE readNormals() aNormal=>aSide%node(1)%np%firstNormal outer:DO WHILE(ASSOCIATED(aNormal)) !search neg FaceID in first node, check all normals. If no neg FaceIDs, side is curved! DO i=1,SIZE(aNormal%FaceID) - IF (aNormal%FaceID(i) .GT. 0) CYCLE !loop till negative FaceID is found in aNormal. + IF (aNormal%FaceID(i) .GT. 0) CYCLE !loop till negative FaceID is found in aNormal. DO j=2,aSide%nNodes !neg FaceID found, search in other nodes for that ID changeIndex=.FALSE. bNormal=>aSide%node(j)%np%firstNormal @@ -532,7 +532,7 @@ SUBROUTINE readNormals() IF (.NOT. changeIndex) EXIT !one node does not have neg FaceID, don't search in other nodes END DO IF (changeIndex) EXIT outer !all nodes have same neg FaceID, side is not curved => don't continue with first Node - END DO + END DO aNormal=>aNormal%nextNormal END DO outer IF (changeIndex) THEN @@ -588,7 +588,7 @@ SUBROUTINE checkNormals() DO WHILE(ASSOCIATED(aElem)) DO j=1, aElem%nNodes IF (aElem%node(j)%np%tmp .NE. -1) CYCLE - aElem%node(j)%np%tmp=1 + aElem%node(j)%np%tmp=1 aNormal=>aElem%node(j)%np%firstNormal DO WHILE (ASSOCIATED(aNormal)) nFaceID=min(5,SIZE(aNormal%FaceID)) @@ -629,7 +629,7 @@ END SUBROUTINE getNewNormal SUBROUTINE reconstructNormals() !=================================================================================================================================== -! create Normals on curved boundaries from the actual mesh, so just approximated normals! +! create Normals on curved boundaries from the actual mesh, so just approximated normals! !=================================================================================================================================== ! MODULES USE MOD_Mesh_Vars,ONLY:tElem,tSide,tNode,tNormal @@ -733,7 +733,7 @@ END SUBROUTINE reconstructNormals SUBROUTINE deleteDuplicateNormals() !=================================================================================================================================== -! ? +! ? !=================================================================================================================================== ! MODULES USE MOD_Mesh_Vars,ONLY:tElem,tSide,tNode,tNormal @@ -819,7 +819,7 @@ SUBROUTINE deleteDuplicateNormals() aNormal=>aNormal%nextNormal END DO DEALLOCATE(tempFaceIDs) - END IF !double Normals + END IF !double Normals END DO END IF !curveInd > 0 aSide=>aSide%nextElemSide @@ -830,7 +830,7 @@ END SUBROUTINE deleteDuplicateNormals SUBROUTINE buildCurvedElementsFromVolume() !=================================================================================================================================== -! set surface curvednode pointers from existing curved volume +! set surface curvednode pointers from existing curved volume !=================================================================================================================================== ! MODULES USE MOD_Mesh_Vars,ONLY:tElem,tSide,FirstElem,N @@ -862,7 +862,7 @@ END SUBROUTINE buildCurvedElementsFromVolume SUBROUTINE buildCurvedElementsFromBoundarySides(nCurvedBoundaryLayers) !=================================================================================================================================== -! set surface curvednode pointers from existing curved volume +! set surface curvednode pointers from existing curved volume !=================================================================================================================================== ! MODULES USE MOD_Mesh_Vars,ONLY:tElem,tSide,FirstElem,N,deleteNode @@ -1010,8 +1010,8 @@ SUBROUTINE curvedVolToSurf(Elem,onlyBoundarySides) ! All other elements remain linear (except for elements adjacent to curved edges). !=================================================================================================================================== ! MODULES -USE MOD_Mesh_Vars,ONLY:tElem,tSide,tNodePtr -USE MOD_Mesh_Vars,ONLY:N +USE MOD_Mesh_Vars ,ONLY:tElem,tSide,tNodePtr +USE MOD_Mesh_Vars ,ONLY:N USE MOD_Basis_Vars,ONLY:MapSideToVol,TriaMapInv,QuadMapInv !USE MOD_Basis_Vars,ONLY:TetraMapInv,PyraMapInv,PrismMapInv,HexaMapInv ! IMPLICIT VARIABLE HANDLING @@ -1064,7 +1064,7 @@ END SUBROUTINE curvedVolToSurf SUBROUTINE referenceSideToFlipped(refSide,Side) !=================================================================================================================================== -! set surface curvednode pointers from existing curved volume +! set surface curvednode pointers from existing curved volume !=================================================================================================================================== ! MODULES USE MOD_Mesh_Vars,ONLY:tSide,tNodePtr,N @@ -1089,7 +1089,7 @@ SUBROUTINE referenceSideToFlipped(refSide,Side) CASE(3) Side%nCurvedNodes=((N+1)*(N+2))/2 ALLOCATE(Side%curvedNode(Side%nCurvedNodes)) - DO q=0,N + DO q=0,N DO p=0,N-q SELECT CASE(flip) CASE(0) @@ -1100,14 +1100,14 @@ SUBROUTINE referenceSideToFlipped(refSide,Side) Side%curvedNode(TriaMapInv(p,q))%np => refSide(N-q-p,q)%np CASE(3) Side%curvedNode(TriaMapInv(p,q))%np => refSide(p,N-q-p)%np - END SELECT + END SELECT Side%curvedNode(TriaMapInv(p,q))%np%refCount=Side%curvedNode(TriaMapInv(p,q))%np%refCount+1 END DO !p END DO !q CASE(4) Side%nCurvedNodes=((N+1)**2) ALLOCATE(Side%curvedNode(Side%nCurvedNodes)) - DO q=0,N + DO q=0,N DO p=0,N SELECT CASE(flip) CASE(0) @@ -1177,7 +1177,7 @@ END FUNCTION ElemIsCurved SUBROUTINE curvedSurfToEdges(aSide) !=================================================================================================================================== -! set edge curvednode pointers from existing curved surface +! set edge curvednode pointers from existing curved surface !=================================================================================================================================== ! MODULES USE MOD_Basis_Vars,ONLY:EdgeToTria,EdgeToQuad @@ -1219,7 +1219,7 @@ SUBROUTINE curvedSurfToEdges(aSide) END DO END IF END IF -END DO !iEdge +END DO !iEdge END SUBROUTINE curvedSurfToEdges @@ -1386,7 +1386,7 @@ SUBROUTINE ProjectToExactSurfaces() EXIT END IF END DO - aSide%tmp=exactFunction + aSide%tmp=exactFunction DO i=1,aSide%nNodes aSide%node(i)%np%tmp=1 END DO @@ -1466,12 +1466,12 @@ SUBROUTINE exactSurfaceFunction(xold,exactFunction,xnew) xnew(1:2)=(xold(1:2))*R/rloc xnew(3)=xold(3) CASE(11) !naca0012 profile,zero TE thickness,c=1, - ! x=t^2,y=sign(yold)*0.12/0.2*(0.2969*t-0.1260*t**2-0.3516*t**4+0.2843*t**6-0.1036*t**8) + ! x=t^2,y=sign(yold)*0.12/0.2*(0.2969*t-0.1260*t**2-0.3516*t**4+0.2843*t**6-0.1036*t**8) c(0:4)=(/0.2969,-0.1260,-0.3516,0.2843,-0.1036/) t=MIN(1.,SQRT(MAX(0.,xold(1)))) !start for newton iteration d=0.12/0.2 !find t with minimum distance to xold (xt-x)^2+(yt-y)^2 -> min - DO i=1,1000 + DO i=1,1000 tt(1)=t DO j=2,8 tt(j)=tt(j-1)*t @@ -1490,7 +1490,7 @@ SUBROUTINE exactSurfaceFunction(xold,exactFunction,xnew) xnew(1)=xt xnew(2)=yt*SIGN(1.,xold(2)) xnew(3)=xold(3) - IF(i.GE.1000) THEN + IF(i.GE.1000) THEN WRITE(UNIT_stdOut,'(A)')' Warning: Newton iteration not converged in exactSurfaceFunction:' WRITE(UNIT_stdOut,'(A,3E21.11)')' old position',xold WRITE(UNIT_stdOut,'(A,3E21.11)')' new position',xnew @@ -1503,7 +1503,7 @@ END SUBROUTINE exactSurfaceFunction SUBROUTINE create3DSplines() !=================================================================================================================================== -! Creates Spline Coefficients. +! Creates Spline Coefficients. !=================================================================================================================================== ! MODULES USE MOD_Mesh_Vars,ONLY:tElem,tSide,tEdge @@ -1542,7 +1542,7 @@ SUBROUTINE create3DSplines() aEdge=>aSide%Edge(i)%edp IF(ASSOCIATED(aEdge%curvedNode))CYCLE ALLOCATE(aEdge%curvedNode(4)) - CALL getTangentialVectors(aSide,i,v,normalCaseCount) + CALL getTangentialVectors(aSide,i,v,normalCaseCount) aEdge%CurvedNode(1)%np=>aEdge%Node(1)%np aEdge%CurvedNode(4)%np=>aEdge%Node(2)%np CALL getNewNode(aEdge%CurvedNode(2)%np,1) @@ -1711,7 +1711,7 @@ SUBROUTINE getTangentialVectors(aSide,nodeInd,v,normalCaseCount) aNormal(2)%np=>aNormal(2)%np%nextNormal END DO ERRWRITE(*,'(132("~"))') - + prevNodeInd=nodeInd-1 IF (prevNodeInd .EQ. 0) prevNodeInd=aSide%nNodes vTemp(:)=aSide%OrientedNode(prevNodeInd)%np%x(:)-aSide%OrientedNode(nodeInd)%np%x(:) !calculate normal of aSide @@ -1736,7 +1736,7 @@ SUBROUTINE getTangentialVectors(aSide,nodeInd,v,normalCaseCount) vTemp=cross(aNormal(1)%np%normal,aNormal(2)%np%normal) vTemp=vTemp/NORM2(vTemp) AngleTmp=ABS(SUM(vTemp(:)*approxNormal(:))) !angle to approxnormal not too big - IF (AngleTmp .LE. 0.4) THEN + IF (AngleTmp .LE. 0.4) THEN AngleTmp=ABS(SUM(vTemp(:)*sideVect(:,n)))/length ! result is cos of angle IF (AngleTmp .GT. AngleMin) THEN AngleMin=AngleTmp @@ -1748,19 +1748,19 @@ SUBROUTINE getTangentialVectors(aSide,nodeInd,v,normalCaseCount) END DO aNormal(1)%np=>aNormal(1)%np%nextNormal END DO - IF ((AngleMin .GT. AngleLimit) .AND. (ASSOCIATED(foundNormals(1,1)%np))) THEN !if small angle use cross product + IF ((AngleMin .GT. AngleLimit) .AND. (ASSOCIATED(foundNormals(1,1)%np))) THEN !if small angle use cross product normalCaseCount(2) = normalCaseCount(2)+1 vTemp(:)=cross(foundNormals(1,1)%np%normal,foundNormals(1,2)%np%normal) vTemp(:)=vTemp(:)/SQRT(SUM(vTemp(:)*vTemp(:))) v(:,n)=vTemp(:)*SUM(vTemp(:)*sideVect(:,n)) tryProjection=.FALSE. ELSE !if big angle use projection - tryProjection=.TRUE. + tryProjection=.TRUE. END IF ELSE tryProjection=.TRUE. END IF - + IF (tryProjection) THEN normalCaseCount(1) = normalCaseCount(1)+1 AngleMin=-1. !cos of angle 1=min, 0=max @@ -1783,7 +1783,7 @@ END SUBROUTINE getTangentialVectors SUBROUTINE curvedEdgesToSurf(keepExistingCurveds) !=================================================================================================================================== -! Blend edges of a side to a curved side +! Blend edges of a side to a curved side !=================================================================================================================================== ! MODULES USE MOD_Mesh_Vars,ONLY:tElem,tSide,FirstElem @@ -1842,7 +1842,7 @@ SUBROUTINE curvedEdgesToSurf(keepExistingCurveds) ! compute the splines between the edges IF(Side%nNodes.EQ.3)THEN CALL curvedEdgesToTriaSurf(Side) - ELSE + ELSE CALL curvedEdgesToQuadSurf(Side) END IF @@ -1912,7 +1912,7 @@ SUBROUTINE curvedEdgesToTriaSurf(Side) REAL,DIMENSION(3) :: xa,xb,xc,xt,sN ! ? !=================================================================================================================================== sN=1./REAL(N) -! allocate Side +! allocate Side nAns=(N+1)*(N+2)/2 Side%nCurvedNodes=nAns ALLOCATE(Side%CurvedNode(nAns)) @@ -1920,7 +1920,7 @@ SUBROUTINE curvedEdgesToTriaSurf(Side) DO i=1,nAns NULLIFY(xpos(i)%np) END DO -! fill edge points +! fill edge points DO iEdge=1,Side%nNodes Edge=>Side%Edge(iEdge)%edp IF(.NOT.ASSOCIATED(Edge%CurvedNode))THEN !create curved nodes for linear edge @@ -1933,7 +1933,7 @@ SUBROUTINE curvedEdgesToTriaSurf(Side) END DO END IF - IF(Side%EdgeOrientation(iEdge))THEN !oriented + IF(Side%EdgeOrientation(iEdge))THEN !oriented DO i=1,N xpos(EdgeToTria(iEdge,i))%np=>Edge%CurvedNode(i)%np END DO @@ -1949,8 +1949,8 @@ SUBROUTINE curvedEdgesToTriaSurf(Side) DO i=1,N-j-1 CALL getNewNode(xpos(TriaMapInv(i,j))%np) xa=(REAL(N-i-j)*xpos(TriaMapInv(i, 0))%np%x +REAL(j)*xpos(TriaMapInv(i, N-i))%np%x)/REAL(N-i) - xb=(REAL(N-i-j)*xpos(TriaMapInv(0, j))%np%x +REAL(i)*xpos(TriaMapInv(N-j,j ))%np%x)/REAL(N-j) - xc=(REAL(i) *xpos(TriaMapInv(i+j,0))%np%x +REAL(j)*xpos(TriaMapInv(0, i+j))%np%x)/REAL(i+j) + xb=(REAL(N-i-j)*xpos(TriaMapInv(0, j))%np%x +REAL(i)*xpos(TriaMapInv(N-j,j ))%np%x)/REAL(N-j) + xc=(REAL(i) *xpos(TriaMapInv(i+j,0))%np%x +REAL(j)*xpos(TriaMapInv(0, i+j))%np%x)/REAL(i+j) xt=(REAL(N-i-j)*xpos(TriaMapInv(0, 0))%np%x +REAL(i)*xpos(TriaMapInv(N, 0 ))%np%x+ & REAL(j)*xpos(TriaMapInv(0,N))%np%x)*sN xpos(TriaMapInv(i,j))%np%x=0.5*(xa+xb+xc-xt) @@ -1961,7 +1961,7 @@ SUBROUTINE curvedEdgesToTriaSurf(Side) Side%CurvedNode(i)%np=>xpos(i)%np NULLIFY(xpos(i)%np) END DO -DEALLOCATE(xpos) +DEALLOCATE(xpos) END SUBROUTINE curvedEdgesToTriaSurf @@ -1989,7 +1989,7 @@ SUBROUTINE curvedEdgesToQuadSurf(aSide) !=================================================================================================================================== sN=1./REAL(N) nNodes=aSide%nNodes -! allocate aSide +! allocate aSide nAns=(N+1)**2 aSide%nCurvedNodes=nAns ALLOCATE(aSide%CurvedNode(nAns)) @@ -2011,8 +2011,8 @@ SUBROUTINE curvedEdgesToQuadSurf(aSide) aEdge%curvedNode(i)%np%x=REAL(i-1)*sN*(aEdge%Node(2)%np%x-aEdge%Node(1)%np%x)+aEdge%Node(1)%np%x END DO END IF - - IF(aSide%EdgeOrientation(iEdge))THEN !oriented + + IF(aSide%EdgeOrientation(iEdge))THEN !oriented DO i=1,N xpos(EdgeToQuad(iEdge,i))%np=>aEdge%CurvedNode(i)%np END DO @@ -2041,12 +2041,12 @@ SUBROUTINE curvedEdgesToQuadSurf(aSide) aSide%CurvedNode(i)%np=>xpos(i)%np NULLIFY(xpos(i)%np) END DO -DEALLOCATE(xpos) +DEALLOCATE(xpos) END SUBROUTINE curvedEdgesToQuadSurf SUBROUTINE curvedSurfacesToElem() !=================================================================================================================================== -! Blend curved surfaces of an element to a curved volume (TODO: currently only hexas working) +! Blend curved surfaces of an element to a curved volume (TODO: currently only hexas working) !=================================================================================================================================== ! MODULES USE MOD_Mesh_Vars,ONLY:tElem,FirstElem @@ -2103,7 +2103,7 @@ SUBROUTINE curvedSurfacesToHexa(Elem) INTEGER :: i,j,k,nAns,p,q,flip ! ? REAL,DIMENSION(3) :: xFaces,xEdges,xNodes,xi,rxi ! ? !=================================================================================================================================== -! allocate +! allocate nAns=(N+1)**3 Elem%nCurvedNodes=nAns ALLOCATE(Elem%CurvedNode(nAns)) @@ -2121,9 +2121,9 @@ SUBROUTINE curvedSurfacesToHexa(Elem) xpos(N,0,N)%np=>Elem%Node(6)%np xpos(N,N,N)%np=>Elem%Node(7)%np xpos(0,N,N)%np=>Elem%Node(8)%np -! fill side points +! fill side points Side=>Elem%firstSide -DO WHILE(ASSOCIATED(Side)) +DO WHILE(ASSOCIATED(Side)) flip=getFlip(Side) DO q=0,N; DO p=0,N @@ -2139,23 +2139,23 @@ SUBROUTINE curvedSurfacesToHexa(Elem) Xside(p,q)%np=>Side%CurvedNode(QuadMapInv(N-q,N-p))%np CASE(4) Xside(p,q)%np=>Side%CurvedNode(QuadMapInv(p,N-q))%np - END SELECT + END SELECT END DO; END DO ! p,q DO q=0,N; DO p=0,N SELECT CASE(Side%LocSide) CASE(5) !XI_MINUS) IF(.NOT.ASSOCIATED(xpos(0,p,q)%np)) & - xpos(0,p,q)%np => Xside(q,p)%np + xpos(0,p,q)%np => Xside(q,p)%np CASE(3) !XI_PLUS) IF(.NOT.ASSOCIATED(xpos(N,p,q)%np)) & xpos(N,p,q)%np => Xside(p,q)%np CASE(2) !ETA_MINUS) IF(.NOT.ASSOCIATED(xpos(p,0,q)%np)) & - xpos(p,0,q)%np => Xside(p,q)%np + xpos(p,0,q)%np => Xside(p,q)%np CASE(4) !ETA_PLUS) IF(.NOT.ASSOCIATED(xpos(p,N,q)%np)) & - xpos(p,N,q)%np => Xside(N-p,q)%np + xpos(p,N,q)%np => Xside(N-p,q)%np CASE(1) !ZETA_MINUS) IF(.NOT.ASSOCIATED(xpos(p,q,0)%np)) & xpos(p,q,0)%np => Xside(q,p)%np @@ -2171,13 +2171,13 @@ SUBROUTINE curvedSurfacesToHexa(Elem) DO k=1,N-1; DO j=1,N-1; DO i=1,N-1 CALL getNewNode(xpos(i,j,k)%np,0) xi = REAL((/i,j,k/))/REAL(N) - rxi = 1.-xi + rxi = 1.-xi xFaces = rxi(1)*xpos( 0, j, k)%np%x + & xi(1)*xpos( N, j, k)%np%x + & rxi(2)*xpos( i, 0, k)%np%x + & xi(2)*xpos( i, N, k)%np%x + & rxi(3)*xpos( i, j, 0)%np%x + & - xi(3)*xpos( i, j, N)%np%x + xi(3)*xpos( i, j, N)%np%x xEdges = rxi(1)*rxi(2)*xpos( 0, 0, k)%np%x + & rxi(1)* xi(2)*xpos( 0, N, k)%np%x + & @@ -2208,7 +2208,7 @@ SUBROUTINE curvedSurfacesToHexa(Elem) Elem%CurvedNode(HexaMapInv(i,j,k))%np=>xpos(i,j,k)%np NULLIFY(xpos(i,j,k)%np) END DO; END DO; END DO -DEALLOCATE(xpos,xside) +DEALLOCATE(xpos,xside) END SUBROUTINE curvedSurfacesToHexa @@ -2220,7 +2220,7 @@ SUBROUTINE SplitToSpline() ! 2) for each element side having a curveind>0 => aSide ! a) find a split element attached to the first node of the aSide. Until now it is only known that they share a node, we ! have to find out if the other corners of the element side are corresponding too... -! b) using the connections between the splitElements, we find all splitElements belonging to one side +! b) using the connections between the splitElements, we find all splitElements belonging to one side ! .p.e. for quads we have a (bOrd-1)*(bOrd-1) splitelement array (the 'left' side of each element is saved=>localSides) ! c) if we can buld up the whole array (reachedcorner=T), we can check, if the other aSide Nodes share corners ! d) then the points areinterpolated to the monomial spline basis of aSide by a inverted Vandermonde matrix (sVdMquad). @@ -2289,8 +2289,8 @@ SUBROUTINE SplitToSpline() END DO END IF !found aElem=>aElem%nextElem -END DO -! enlarge box +END DO +! enlarge box xmin=xmin-SQRT(dxmax) xmax=xmax+SQRT(dxmax) !build local search mesh of split sides @@ -2309,14 +2309,14 @@ SUBROUTINE SplitToSpline() END IF aSplitElem=>aSplitElem%nextElem END DO -WRITE(*,*)'number of splitElems in searchmesh',nSplitElems +WRITE(*,*)'number of splitElems in searchmesh',nSplitElems !========================================== Start search of split elements ===================================================== aElem=>firstElem DO WHILE(ASSOCIATED(aElem)) aSide=>aElem%firstSide DO WHILE(ASSOCIATED(aSide)) - IF(aSide%curveIndex.NE.0) THEN + IF(aSide%curveIndex.NE.0) THEN aNode=>aSide%OrientedNode(1)%np !local element side tolerance tol=SQRT(SUM((aSide%OrientedNode(2)%np%x-aSide%OrientedNode(1)%np%x)**2)) @@ -2339,7 +2339,7 @@ SUBROUTINE SplitToSpline() DO k=1,aSplitElem%nNodes IF (SAMEPOINT(aSplitSide%Node(2)%np%x,aNode%x,tol)) THEN pointFound=.TRUE. - EXIT + EXIT END IF aSplitSide=>aSplitSide%nextElemSide END DO @@ -2354,18 +2354,18 @@ SUBROUTINE SplitToSpline() localSides(1,1)%sp=>aSplitSide ! save left lying sides in local sides, bOrd=3 --> 4 sides ! - ! ^ j-dir - ! | - ! | + ! ^ j-dir + ! | + ! | ! (4)---(+)---(3) - ! | | + ! | | ! (+)---(+)---(+) - ! | | + ! | | ! (1)---(+)---(2) --> i-dir ! tria=0 IF(aSplitElem%nNodes.EQ.3) tria=1 -! WRITE(*,'(A,6F10.5)')'DEBUG tria',tria +! WRITE(*,'(A,6F10.5)')'DEBUG tria',tria !now fill up the whole array outer: DO j=1,bOrd-1 ! j direction, for j=1 do nothing @@ -2376,11 +2376,11 @@ SUBROUTINE SplitToSpline() DO l=1,1+tria ! l=1,1 for quads, l=1,2 for trias aSplitSide=>aSplitSide%Connection IF(.NOT.ASSOCIATED(aSplitSide)) THEN - reachedcorner=.FALSE. + reachedcorner=.FALSE. EXIT !no corner found, because no connection (end of local domain) END IF aSplitSide=>GoToSide(aSplitSide,-1) - END DO + END DO localSides(1,j)%sp=>aSplitSide END IF DO i=2,bOrd-1-(j-1)*tria @@ -2390,7 +2390,7 @@ SUBROUTINE SplitToSpline() !go to neighbor element aSplitSide=>aSplitSide%connection IF(.NOT.ASSOCIATED(aSplitSide)) THEN - reachedcorner=.FALSE. + reachedcorner=.FALSE. EXIT outer !no corner found, because no connection (end of local domain) END IF END DO @@ -2400,8 +2400,8 @@ SUBROUTINE SplitToSpline() END DO outer !j IF(.NOT.reachedcorner)THEN WRITE(*,*)'DEBUG corner not reached, no mpi -> PROBLEM!' - END IF ! corner not reached - + END IF ! corner not reached + IF(reachedcorner)THEN ! all local Sides are associated, check if localSides belong to aSide nodekminus=>localSides(1,bord-1)%sp%node(1)%np @@ -2422,14 +2422,14 @@ SUBROUTINE SplitToSpline() DO j=1,bOrd NULLIFY(localNodes(i,j)%np) END DO - END DO + END DO DO j=1,bOrd-1 DO i=1,bOrd-1-(j-1)*tria localNodes(i,j)%np=>localSides(i,j)%sp%Node(2)%np END DO aSplitSide=>GoToSide(localSides(bOrd-1-(j-1)*tria,j)%sp,1) localNodes(bOrd-(j-1)*tria,j)%np=>aSplitSide%Node(2)%np - END DO + END DO IF(tria.EQ.0)THEN DO i=1,bOrd-1 localNodes(i,bOrd)%np=>localSides(i,bOrd-1)%sp%Node(1)%np @@ -2506,19 +2506,19 @@ SUBROUTINE SplitToSpline() aSide=>aSide%nextElemSide END DO aElem=>aElem%nextElem -END DO +END DO CALL deleteSearchMesh(SearchMesh,.FALSE.) CALL deleteSearchMesh(redundantNodes,.TRUE.) -!delete all split elements and nodes +!delete all split elements and nodes aSplitElem=>firstSplitElem -DO WHILE (ASSOCIATED(aSplitElem)) +DO WHILE (ASSOCIATED(aSplitElem)) DO i=1,aSplitElem%nNodes aSplitElem%Node(i)%np%refcount=0 END DO aSplitElem=>aSplitElem%nextElem END DO aSplitElem=>firstSplitElem -DO WHILE (ASSOCIATED(aSplitElem)) +DO WHILE (ASSOCIATED(aSplitElem)) DO i=1,aSplitElem%nNodes aSplitElem%Node(i)%np%refcount=aSplitElem%Node(i)%np%refcount+1 END DO @@ -2552,7 +2552,7 @@ SUBROUTINE SplitToSpline() i=i+1 DEALLOCATE(aSplitElem) aSplitElem=>firstSplitElem -END DO +END DO DEALLOCATE(localNodes,localSides) WRITE(*,'(A,I8)')'nSplitElems deleted',i WRITE(*,'(A,I8)')'nNodes of SplitElems deleted',l @@ -2562,7 +2562,7 @@ END SUBROUTINE SplitToSpline FUNCTION GoToSide(Side,nJumps) !=================================================================================================================================== -! repointer side to the next nJumps +! repointer side to the next nJumps !=================================================================================================================================== ! MODULES USE MOD_Mesh_Vars,ONLY:tSide @@ -2592,7 +2592,7 @@ FUNCTION GoToSide(Side,nJumps) IF(.NOT.ASSOCIATED(GoToSide))GoToSide=>Side%Elem%firstSide END DO END IF -END FUNCTION GoToSide +END FUNCTION GoToSide SUBROUTINE RebuildMortarGeometry() @@ -2623,7 +2623,7 @@ SUBROUTINE RebuildMortarGeometry() ALLOCATE(M_0_1_T(0:N,0:N)) ALLOCATE(M_0_2_T(0:N,0:N)) -CALL GetMortarVandermonde(N, M_0_1_T, M_0_2_T) +CALL GetMortarVandermonde(N, M_0_1_T, M_0_2_T) M_0_1_T=TRANSPOSE(M_0_1_T) M_0_2_T=TRANSPOSE(M_0_2_T) @@ -2812,9 +2812,9 @@ RECURSIVE SUBROUTINE MapBigEdgeToSmall(edge) ELSEIF(ASSOCIATED(Edge%Node(p)%np,smallEdge%Node(3-p)%np))THEN XGeo1DTmp=XGeo1Dsmall(:,:,p) DO l=0,N - XGeo1Dsmall(:,N-l,p)=XGeo1DTmp(:,l) + XGeo1Dsmall(:,N-l,p)=XGeo1DTmp(:,l) END DO - ELSE + ELSE CALL abort(__STAMP__,'Error: Edges of mortar master and slave do not conform!') END IF CALL UnpackGeo(N,XGeo1DSmall(:,:,p),smallEdge) diff --git a/src/mesh/curvedcartmesh.f90 b/src/mesh/curvedcartmesh.f90 index 8376eee..84a028e 100644 --- a/src/mesh/curvedcartmesh.f90 +++ b/src/mesh/curvedcartmesh.f90 @@ -12,7 +12,7 @@ ! Copyright (C) 2017 Claus-Dieter Munz ! This file is part of HOPR, a software for the generation of high-order meshes. ! -! HOPR is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License +! HOPR is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License ! as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. ! ! HOPR is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty @@ -30,7 +30,7 @@ MODULE MOD_CurvedCartMesh IMPLICIT NONE PRIVATE !----------------------------------------------------------------------------------------------------------------------------------- -! GLOBAL VARIABLES +! GLOBAL VARIABLES !----------------------------------------------------------------------------------------------------------------------------------- INTERFACE CurvedCartesianMesh @@ -41,8 +41,18 @@ MODULE MOD_CurvedCartMesh MODULE PROCEDURE GetNewCurvedHexahedron END INTERFACE +INTERFACE GetNewCurvedPrism + MODULE PROCEDURE GetNewCurvedPrism +END INTERFACE + +INTERFACE GetNewCurvedPyram + MODULE PROCEDURE GetNewCurvedPyram +END INTERFACE + PUBLIC::CurvedCartesianMesh PUBLIC::GetNewCurvedHexahedron +PUBLIC::GetNewCurvedPrism +PUBLIC::GetNewCurvedPyram !=================================================================================================================================== CONTAINS @@ -102,6 +112,277 @@ SUBROUTINE GetNewCurvedHexahedron(CurvedNode,Ngeo,Zone) END SUBROUTINE GetNewCurvedHexahedron +SUBROUTINE GetNewCurvedPrism(CurvedNode,Ngeo,Zone) +!=================================================================================================================================== +! Build new prism for cartesian mesh. +!=================================================================================================================================== +! MODULES +USE MOD_Mesh_Vars ,ONLY:tNodePtr,tElem +USE MOD_Mesh_Vars ,ONLY:FirstElem +USE MOD_Mesh_Vars ,ONLY:getNewElem +USE MOD_Mesh_Basis ,ONLY:CreateSides +USE MOD_Basis_Vars ,ONLY:PrismMap,PrismMap2 +! IMPLICIT VARIABLE HANDLING +IMPLICIT NONE +!----------------------------------------------------------------------------------------------------------------------------------- +! INPUT VARIABLES +INTEGER,INTENT(IN) :: NGeo ! ? +INTEGER,INTENT(IN) :: Zone ! ? +TYPE(tNodePtr),INTENT(IN) :: CurvedNode(0:Ngeo,0:NGeo,0:Ngeo) ! ? +!----------------------------------------------------------------------------------------------------------------------------------- +! OUTPUT VARIABLES +!----------------------------------------------------------------------------------------------------------------------------------- +! LOCAL VARIABLES +TYPE(tElem),POINTER :: aElem ! ? +INTEGER :: i ! ? +!=================================================================================================================================== +CALL getNewElem(aElem) +aElem%nNodes=6 +ALLOCATE(aElem%Node(aElem%nNodes)) +aElem%Node(1)%NP=>CurvedNode( 0, 0, 0)%NP +aElem%Node(2)%NP=>CurvedNode(NGeo, 0, 0)%NP +aElem%Node(3)%NP=>CurvedNode( 0,NGeo, 0)%NP +aElem%Node(4)%NP=>CurvedNode( 0, 0,NGeo)%NP +aElem%Node(5)%NP=>CurvedNode(NGeo, 0,NGeo)%NP +aElem%Node(6)%NP=>CurvedNode( 0,NGeo,NGeo)%NP +aElem%splitind = 1 +CALL CreateSides(aElem,.TRUE.) +IF(Ngeo.GT.1)THEN !curved + aElem%nCurvedNodes=NINT((Ngeo+2)*(Ngeo+1)**2/2.) + ALLOCATE(aElem%CurvedNode(aElem%nCurvedNodes)) + DO i=1,aElem%nCurvedNodes + aElem%CurvedNode(i)%NP=>CurvedNode(PrismMap(i,1),PrismMap(i,2),PrismMap(i,3))%NP + !print *,'Elem1',PrismMap(i,:), aElem%CurvedNode(i)%NP%x(:) + END DO +END IF !curved + +CALL getNewElem(aElem%nextElem) +aElem%nextElem%prevElem => aElem +aElem => aElem%nextElem +aElem%nNodes=6 +ALLOCATE(aElem%Node(aElem%nNodes)) +aElem%Node(1)%NP=>CurvedNode(NGeo, 0, 0)%NP +aElem%Node(2)%NP=>CurvedNode(NGeo,NGeo, 0)%NP +aElem%Node(3)%NP=>CurvedNode( 0,NGeo, 0)%NP +aElem%Node(4)%NP=>CurvedNode(NGeo, 0,NGeo)%NP +aElem%Node(5)%NP=>CurvedNode(NGeo,NGeo,NGeo)%NP +aElem%Node(6)%NP=>CurvedNode( 0,NGeo,NGeo)%NP +aElem%splitind = 2 +CALL CreateSides(aElem,.TRUE.) +IF(Ngeo.GT.1)THEN !curved + aElem%nCurvedNodes=NINT((Ngeo+2)*(Ngeo+1)**2/2.) + ALLOCATE(aElem%CurvedNode(aElem%nCurvedNodes)) + DO i=1,aElem%nCurvedNodes + aElem%CurvedNode(i)%NP=>CurvedNode(PrismMap2(i,1),PrismMap2(i,2),PrismMap2(i,3))%NP + !print *,'Elem2',PrismMap2(i,:), aElem%CurvedNode(i)%NP%x(:) + END DO +END IF !curved + +! Add elements to list +IF(.NOT.ASSOCIATED(FirstElem))THEN + FirstElem=>aElem + DO WHILE(ASSOCIATED(FirstElem%prevElem)) + FirstElem=>FirstElem%prevElem + END DO +ELSE + aElem%nextElem => FirstElem + aElem%nextElem%prevElem => aElem + FirstElem => aElem + DO WHILE(ASSOCIATED(FirstElem%prevElem)) + FirstElem=>FirstElem%prevElem + END DO +END IF +NULLIFY(aElem) + +END SUBROUTINE GetNewCurvedPrism + + +SUBROUTINE GetNewCurvedPyram(CurvedNode,Ngeo,Zone) +!=================================================================================================================================== +! Build new prism for cartesian mesh. +!=================================================================================================================================== +! MODULES +USE MOD_Mesh_Vars ,ONLY:tNodePtr,tElem,tNode +USE MOD_Mesh_Vars ,ONLY:FirstElem +USE MOD_Mesh_Vars ,ONLY:getNewElem,getNewNode +USE MOD_Mesh_Basis ,ONLY:CreateSides +USE MOD_Basis_Vars ,ONLY:PyraMap +! IMPLICIT VARIABLE HANDLING +IMPLICIT NONE +!----------------------------------------------------------------------------------------------------------------------------------- +! INPUT VARIABLES +INTEGER,INTENT(IN) :: NGeo ! ? +INTEGER,INTENT(IN) :: Zone ! ? +TYPE(tNodePtr),INTENT(IN) :: CurvedNode(0:Ngeo,0:NGeo,0:Ngeo) ! ? +!----------------------------------------------------------------------------------------------------------------------------------- +! OUTPUT VARIABLES +!----------------------------------------------------------------------------------------------------------------------------------- +! LOCAL VARIABLES +TYPE(tElem),POINTER :: aElem ! ? +TYPE(tNode),POINTER :: CenterNode ! ? +REAL :: xm(3) +REAL :: CornerNode(8,3) ! ? +INTEGER :: i ! ? +!=================================================================================================================================== +! Get center node +CornerNode(1,:)=CurvedNode( 0, 0, 0)%NP%x +CornerNode(2,:)=CurvedNode(NGeo, 0, 0)%NP%x +CornerNode(3,:)=CurvedNode(NGeo,NGeo, 0)%NP%x +CornerNode(4,:)=CurvedNode( 0,NGeo, 0)%NP%x +CornerNode(5,:)=CurvedNode( 0, 0,NGeo)%NP%x +CornerNode(6,:)=CurvedNode(NGeo, 0,NGeo)%NP%x +CornerNode(7,:)=CurvedNode(NGeo,NGeo,NGeo)%NP%x +CornerNode(8,:)=CurvedNode( 0,NGeo,NGeo)%NP%x + +xm=0. +DO i=1,8 + xm=xm+CornerNode(i,:) +END DO +xm=0.125*xm +CALL getNewNode(CenterNode) +CenterNode%x=xm + +! 1. Pyramid +CALL getNewElem(aElem) +aElem%nNodes=5 +ALLOCATE(aElem%Node(aElem%nNodes)) +aElem%Node(1)%NP=>CurvedNode( 0, 0, 0)%NP +aElem%Node(2)%NP=>CurvedNode(NGeo, 0, 0)%NP +aElem%Node(3)%NP=>CurvedNode(NGeo,NGeo, 0)%NP +aElem%Node(4)%NP=>CurvedNode( 0,NGeo, 0)%NP +aElem%Node(5)%NP=>CenterNode +CALL CreateSides(aElem,.TRUE.) +IF(Ngeo.GT.1)THEN !curved + aElem%nCurvedNodes=NINT((NGeo+1)*(NGeo+2)*(2*NGeo+3)/6.) + ALLOCATE(aElem%CurvedNode(aElem%nCurvedNodes)) + DO i=1,aElem%nCurvedNodes + aElem%CurvedNode(i)%NP=> CurvedNode(PyraMap(i,1),PyraMap(i,2),PyraMap(i,3))%NP + END DO +END IF !curved + +! 2. Pyramid +CALL getNewElem(aElem%nextElem) +aElem%nextElem%prevElem => aElem +aElem => aElem%nextElem +aElem%nNodes=5 +ALLOCATE(aElem%Node(aElem%nNodes)) +aElem%Node(1)%NP=>CurvedNode( 0, 0, 0)%NP +aElem%Node(2)%NP=>CurvedNode( 0, 0,NGeo)%NP +aElem%Node(3)%NP=>CurvedNode(NGeo, 0,NGeo)%NP +aElem%Node(4)%NP=>CurvedNode(NGeo, 0, 0)%NP +aElem%Node(5)%NP=>CenterNode +CALL CreateSides(aElem,.TRUE.) +IF(Ngeo.GT.1)THEN !curved + aElem%nCurvedNodes=NINT((NGeo+1)*(NGeo+2)*(2*NGeo+3)/6.) + ALLOCATE(aElem%CurvedNode(aElem%nCurvedNodes)) + DO i=1,aElem%nCurvedNodes + aElem%CurvedNode(i)%NP=>CurvedNode(PyraMap(i,1),PyraMap(i,2),PyraMap(i,3))%NP + END DO +END IF !curved + +! 3. Pyramid +CALL getNewElem(aElem%nextElem) +aElem%nextElem%prevElem => aElem +aElem => aElem%nextElem +aElem%nNodes=5 +ALLOCATE(aElem%Node(aElem%nNodes)) +aElem%Node(1)%NP=>CurvedNode(NGeo, 0, 0)%NP +aElem%Node(2)%NP=>CurvedNode(NGeo, 0,NGeo)%NP +aElem%Node(3)%NP=>CurvedNode(NGeo,NGeo,NGeo)%NP +aElem%Node(4)%NP=>CurvedNode(NGeo,NGeo, 0)%NP +aElem%Node(5)%NP=>CenterNode +CALL CreateSides(aElem,.TRUE.) +IF(Ngeo.GT.1)THEN !curved + aElem%nCurvedNodes=NINT((NGeo+1)*(NGeo+2)*(2*NGeo+3)/6.) + ALLOCATE(aElem%CurvedNode(aElem%nCurvedNodes)) + DO i=1,aElem%nCurvedNodes + aElem%CurvedNode(i)%NP=>CurvedNode(PyraMap(i,1),PyraMap(i,2),PyraMap(i,3))%NP + !print *,'Elem2',PrismMap2(i,:), aElem%CurvedNode(i)%NP%x(:) + END DO +END IF !curved + +! 4. Pyramid +CALL getNewElem(aElem%nextElem) +aElem%nextElem%prevElem => aElem +aElem => aElem%nextElem +aElem%nNodes=5 +ALLOCATE(aElem%Node(aElem%nNodes)) +aElem%Node(1)%NP=>CurvedNode( 0, 0, 0)%NP +aElem%Node(2)%NP=>CurvedNode( 0,NGeo, 0)%NP +aElem%Node(3)%NP=>CurvedNode( 0,NGeo,NGeo)%NP +aElem%Node(4)%NP=>CurvedNode( 0, 0,NGeo)%NP +aElem%Node(5)%NP=>CenterNode +CALL CreateSides(aElem,.TRUE.) +IF(Ngeo.GT.1)THEN !curved + aElem%nCurvedNodes=NINT((NGeo+1)*(NGeo+2)*(2*NGeo+3)/6.) + ALLOCATE(aElem%CurvedNode(aElem%nCurvedNodes)) + DO i=1,aElem%nCurvedNodes + aElem%CurvedNode(i)%NP=>CurvedNode(PyraMap(i,1),PyraMap(i,2),PyraMap(i,3))%NP + !print *,'Elem2',PrismMap2(i,:), aElem%CurvedNode(i)%NP%x(:) + END DO +END IF !curved + +! 5. Pyramid +CALL getNewElem(aElem%nextElem) +aElem%nextElem%prevElem => aElem +aElem => aElem%nextElem +aElem%nNodes=5 +ALLOCATE(aElem%Node(aElem%nNodes)) +aElem%Node(1)%NP=>CurvedNode( 0, 0,NGeo)%NP +aElem%Node(2)%NP=>CurvedNode( 0,NGeo,NGeo)%NP +aElem%Node(3)%NP=>CurvedNode(NGeo,NGeo,NGeo)%NP +aElem%Node(4)%NP=>CurvedNode(NGeo, 0,NGeo)%NP +aElem%Node(5)%NP=>CenterNode +CALL CreateSides(aElem,.TRUE.) +IF(Ngeo.GT.1)THEN !curved + aElem%nCurvedNodes=NINT((NGeo+1)*(NGeo+2)*(2*NGeo+3)/6.) + ALLOCATE(aElem%CurvedNode(aElem%nCurvedNodes)) + DO i=1,aElem%nCurvedNodes + aElem%CurvedNode(i)%NP=>CurvedNode(PyraMap(i,1),PyraMap(i,2),PyraMap(i,3))%NP + !print *,'Elem2',PrismMap2(i,:), aElem%CurvedNode(i)%NP%x(:) + END DO +END IF !curved + +! 6. Pyramid +CALL getNewElem(aElem%nextElem) +aElem%nextElem%prevElem => aElem +aElem => aElem%nextElem +aElem%nNodes=5 +ALLOCATE(aElem%Node(aElem%nNodes)) +aElem%Node(1)%NP=>CurvedNode(NGeo,NGeo,NGeo)%NP +aElem%Node(2)%NP=>CurvedNode( 0,NGeo,NGeo)%NP +aElem%Node(3)%NP=>CurvedNode( 0,NGeo, 0)%NP +aElem%Node(4)%NP=>CurvedNode(NGeo,NGeo, 0)%NP +aElem%Node(5)%NP=>CenterNode +CALL CreateSides(aElem,.TRUE.) +IF(Ngeo.GT.1)THEN !curved + aElem%nCurvedNodes=NINT((NGeo+1)*(NGeo+2)*(2*NGeo+3)/6.) + ALLOCATE(aElem%CurvedNode(aElem%nCurvedNodes)) + DO i=1,aElem%nCurvedNodes + aElem%CurvedNode(i)%NP=>CurvedNode(PyraMap(i,1),PyraMap(i,2),PyraMap(i,3))%NP + !print *,'Elem2',PrismMap2(i,:), aElem%CurvedNode(i)%NP%x(:) + END DO +END IF !curved + +! Add elements to list +IF(.NOT.ASSOCIATED(FirstElem))THEN + FirstElem=>aElem + DO WHILE(ASSOCIATED(FirstElem%prevElem)) + FirstElem=>FirstElem%prevElem + END DO +ELSE + aElem%nextElem => FirstElem + aElem%nextElem%prevElem => aElem + FirstElem => aElem + DO WHILE(ASSOCIATED(FirstElem%prevElem)) + FirstElem=>FirstElem%prevElem + END DO +END IF +NULLIFY(aElem) + +END SUBROUTINE GetNewCurvedPyram + + SUBROUTINE CurvedCartesianMesh() !=================================================================================================================================== ! Builds cartesian mesh. Called by fillMesh. @@ -119,7 +400,7 @@ SUBROUTINE CurvedCartesianMesh() !----------------------------------------------------------------------------------------------------------------------------------- ! INPUT/OUTPUT VARIABLES !----------------------------------------------------------------------------------------------------------------------------------- -! LOCAL VARIABLES +! LOCAL VARIABLES TYPE(tNodePtr) :: CurvedNode(0:BoundaryOrder-1,0:BoundaryOrder-1,0:BoundaryOrder-1) ! ? TYPE(tNodePtr),POINTER :: Mnodes(:,:,:) ! ? TYPE(tElem),POINTER :: aElem ! ? @@ -156,7 +437,7 @@ SUBROUTINE CurvedCartesianMesh() DO n=0,np(3) CALL GetNewNode(Mnodes(l,m,n)%np) NodeInd=NodeInd+1 - Mnodes(l,m,n)%np%ind=NodeInd + Mnodes(l,m,n)%np%ind=NodeInd END DO !n END DO !m END DO !l @@ -181,8 +462,8 @@ SUBROUTINE CurvedCartesianMesh() DR(iDir)%d(ii)=DR(iDir)%d(ii-1)*fac(iDir) END DO CASE(3) ! bell shaped exp(-s^2), s=[-1,1] - ! DR= 1+ (dxmax/dxmin-2)*(exp(-(s*fac)^2)-exp(-(1*fac)^2))/(exp(-(0*fac)^2)-exp(-(1*fac)^2)) - ! fac>1, more stretching outwards, fac<1 less stretching + ! DR= 1+ (dxmax/dxmin-2)*(exp(-(s*fac)^2)-exp(-(1*fac)^2))/(exp(-(0*fac)^2)-exp(-(1*fac)^2)) + ! fac>1, more stretching outwards, fac<1 less stretching s=-1 ds=2./(nElems(iDir)-1) DO ii=1,nElems(iDir) @@ -191,8 +472,8 @@ SUBROUTINE CurvedCartesianMesh() s=s+ds END DO CASE(4) ! half bell shaped exp(-s^2), s=[-1,0] - ! DR= 1+ (dxmax/dxmin-2)*(exp(-(s*fac)^2)-exp(-(1*fac)^2))/(exp(-(0*fac)^2)-exp(-(1*fac)^2)) - ! fac>1, more stretching outwards, fac<1 less stretching + ! DR= 1+ (dxmax/dxmin-2)*(exp(-(s*fac)^2)-exp(-(1*fac)^2))/(exp(-(0*fac)^2)-exp(-(1*fac)^2)) + ! fac>1, more stretching outwards, fac<1 less stretching ds=(DxMaxToDxmin(iDir)-1.)/(1.-DxMaxToDxmin(iDir)*exp((-1+1./REAL(nElems(iDir)))*fac(iDir))) DO ii=1,nElems(iDir) DR(iDir)%d(ii)=1.+ds*exp((-1.+REAL(ii)/REAL(nElems(iDir)))*fac(iDir)) @@ -223,7 +504,7 @@ SUBROUTINE CurvedCartesianMesh() l=i+(ii-1)*Ngeo m=j+(jj-1)*Ngeo n=k+(kk-1)*Ngeo - Mnodes(l,m,n)%np%x(:)=X_Ngeo(:,i,j,k,ii,jj,kk) + Mnodes(l,m,n)%np%x(:)=X_Ngeo(:,i,j,k,ii,jj,kk) CurvedNode(i,j,k)%np=>Mnodes(l,m,n)%np END DO !i END DO !j @@ -240,7 +521,7 @@ SUBROUTINE CurvedCartesianMesh() DO l=0,np(1) DO m=0,np(2) DO n=0,np(3) - Mnodes(l,m,n)%np%tmp=0 + Mnodes(l,m,n)%np%tmp=0 IF(n.EQ.0) Mnodes(l,m,n)%np%tmp=Mnodes(l,m,n)%np%tmp+1 !zeta minus IF(m.EQ.0) Mnodes(l,m,n)%np%tmp=Mnodes(l,m,n)%np%tmp+20 !eta minus IF(l.EQ.np(1) ) Mnodes(l,m,n)%np%tmp=Mnodes(l,m,n)%np%tmp+300 !xi plus @@ -371,7 +652,7 @@ SUBROUTINE BuildRefDomain(Nin,nElems,R_glob,DR1,DR2,DR3) END DO !jj=1,nElems(2) END DO !kk=1,nElems(3) END IF - + CASE(7) ! sin symmetric mesh Pi = 4.*ATAN(1.) IF(.NOT.((fac(dd).EQ.1.).AND.(fac2(dd).EQ.1))) THEN @@ -385,7 +666,7 @@ SUBROUTINE BuildRefDomain(Nin,nElems,R_glob,DR1,DR2,DR3) r1 = R_glob(dd,i,j,k,ii,jj,kk) !coordinate in stretching direction r2 = R_glob(ee,i,j,k,ii,jj,kk) !perpendicular coordinate tmpfac=0.5*(fac2(dd)*(1+r2)+fac(dd)*(1-r2)) - R_glob(dd,i,j,k,ii,jj,kk)= r1-2*(tmpfac-1.)/(tmpfac+1.)*SIN(Pi*r1)/(2.*Pi) + R_glob(dd,i,j,k,ii,jj,kk)= r1-2*(tmpfac-1.)/(tmpfac+1.)*SIN(Pi*r1)/(2.*Pi) END DO !i=0:Nin END DO !j=0:Nin END DO !k=0:Nin @@ -398,7 +679,7 @@ SUBROUTINE BuildRefDomain(Nin,nElems,R_glob,DR1,DR2,DR3) ! its the normalized integral of the function in Case(3): ! r = \int_{-1}^s 1+(dxmaxtodxmin-1)/(1-e^(-f^2))*( e^(-s^2 f^2)-e^(-f^2) ) ds ! => r= x*[1-e^(-f^2)*(1+(dxmaxtodxmin-1)/(1-e^(-f^2)))] +(dxmaxtodxmin-1)/(1-e^(-f^2))*(\int_{-1}^s e^(-s^2 f^2)ds) - ! the integral of the exponential function is the error function erf, + ! the integral of the exponential function is the error function erf, ! and the integral is the normalized, so that r is in [-1,1] IF(.NOT.(DxMaxToDxMin(dd).EQ.1.)) THEN Pi = 4.*ATAN(1.) @@ -413,7 +694,7 @@ SUBROUTINE BuildRefDomain(Nin,nElems,R_glob,DR1,DR2,DR3) DO j=0,Nin DO i=0,Nin r1 = R_glob(dd,i,j,k,ii,jj,kk) !coordinate in stretching direction - R_glob(dd,i,j,k,ii,jj,kk)=((r1+1.)*c1+c2*(ERF(r1*fac(dd))-ERF(-fac(dd))))/(c1+c2*ERF(fac(dd)))-1. + R_glob(dd,i,j,k,ii,jj,kk)=((r1+1.)*c1+c2*(ERF(r1*fac(dd))-ERF(-fac(dd))))/(c1+c2*ERF(fac(dd)))-1. END DO !i=0:Nin END DO !j=0:Nin END DO !k=0:Nin @@ -436,34 +717,34 @@ SUBROUTINE BuildPhysicalDomain(Nin,nElems,X) IMPLICIT NONE !----------------------------------------------------------------------------------------------------------------------------------- ! INPUT VARIABLES -INTEGER,INTENT(IN) :: Nin,nElems(3) +INTEGER,INTENT(IN) :: Nin,nElems(3) ! We use the CGNS notation of the Hexeader: http://www.grc.nasa.gov/WWW/cgns/sids/conv.html#unst_hexa ! P8------P7 ! /| /| ! P5------P6| ! | P4----|-P3 ! |/ |/ - ! P1------P2 + ! P1------P2 ! Faces - !F1 P1,P4,P3,P2 - !F2 P1,P2,P6,P5 - !F3 P2,P3,P7,P6 - !F4 P3,P4,P8,P7 - !F5 P1,P5,P8,P4 + !F1 P1,P4,P3,P2 + !F2 P1,P2,P6,P5 + !F3 P2,P3,P7,P6 + !F4 P3,P4,P8,P7 + !F5 P1,P5,P8,P4 !F6 P5,P6,P7,P8 ! Edges - !E1 P1,P2 - !E2 P2,P3 - !E3 P3,P4 - !E4 P4,P1 - !E5 P1,P5 - !E6 P2,P6 - !E7 P3,P7 - !E8 P4,P8 - !E9 P5,P6 - !E10 P6,P7 - !E11 P7,P8 - !E12 P8,P5 + !E1 P1,P2 + !E2 P2,P3 + !E3 P3,P4 + !E4 P4,P1 + !E5 P1,P5 + !E6 P2,P6 + !E7 P3,P7 + !E8 P4,P8 + !E9 P5,P6 + !E10 P6,P7 + !E11 P7,P8 + !E12 P8,P5 !----------------------------------------------------------------------------------------------------------------------------------- ! INPUT/OUTPUT VARIABLES REAL,INTENT(INOUT) :: X(3,0:Nin,0:Nin,0:Nin,nElems(1),nElems(2),nElems(3)) ! IN: contains the Ref Domain coords r,s,t @@ -488,7 +769,7 @@ SUBROUTINE BuildPhysicalDomain(Nin,nElems,X) DO j=0,Nin DO i=0,Nin ! overwrite r,s,t with x,y,z - X(:,i,j,k,ii,jj,kk)=X0(:)+(1.+X(:,i,j,K,ii,jj,kk))*DXh(:) + X(:,i,j,k,ii,jj,kk)=X0(:)+(1.+X(:,i,j,K,ii,jj,kk))*DXh(:) END DO !i=0:Nin END DO !j=0:Nin END DO !k=0:Nin @@ -530,7 +811,7 @@ SUBROUTINE BuildPhysicalDomain(Nin,nElems,X) +XP(:,8)*(1.-R_loc(1)) & *(1.+R_loc(2)) & *(1.+R_loc(3)) ) - END DO !i=0:Nin + END DO !i=0:Nin END DO !j=0:Nin END DO !k=0:Nin END DO !ii=1,nElems(1) @@ -556,8 +837,8 @@ SUBROUTINE BuildPhysicalDomain(Nin,nElems,X) ! store the local coordinates r,s,t R_loc = X(:,i,j,k,ii,jj,kk) ! overwrite r,s,t with x,y,z - ! As we have to subtract the contribution of the edges, the contribution of the - ! corner nodes is subtracted two times.... we just add another contribution of + ! As we have to subtract the contribution of the edges, the contribution of the + ! corner nodes is subtracted two times.... we just add another contribution of ! corner points to compensate X(:,i,j,k,ii,jj,kk)=0.125*( XP(:,1)*(1.-R_loc(1)) & *(1.-R_loc(2)) & @@ -598,7 +879,7 @@ SUBROUTINE BuildPhysicalDomain(Nin,nElems,X) +FP(:,3)*(1.+R_loc(1)) & +FP(:,4)*(1.+R_loc(2)) & +FP(:,5)*(1.-R_loc(1)) & - +FP(:,6)*(1.+R_loc(3))) + +FP(:,6)*(1.+R_loc(3))) ! Determine the coordinates of the projection of the ref Point on the 12 Edges EP CALL EvaluateFace(EP(:, 1),(/R_loc(1),-1./),DX,XP,WhichMapping,1) CALL EvaluateFace(EP(:, 2),(/+1.,R_loc(2)/),DX,XP,WhichMapping,1) @@ -638,8 +919,8 @@ SUBROUTINE BuildPhysicalDomain(Nin,nElems,X) +EP(:,11)*(1.+R_loc(2)) & *(1.+R_loc(3)) & +EP(:,12)*(1.-R_loc(1)) & - *(1.+R_loc(3))) - END DO !i=0:Nin + *(1.+R_loc(3))) + END DO !i=0:Nin END DO !j=0:Nin END DO !k=0:Nin END DO !ii=1,nElems(1) @@ -652,18 +933,18 @@ SUBROUTINE BuildPhysicalDomain(Nin,nElems,X) ! loop over Chebyshev points DO k=0,Nin;DO j=0,Nin;DO i=0,Nin ! overwrite r,s,t with x,y,z - X(1,i,j,k,ii,jj,kk)=X0(1)+(fac(1)**(0.5*(1.+X(1,i,j,K,ii,jj,kk)))-1.)*DXh(1) - END DO; END DO; END DO - END DO; END DO; END DO + X(1,i,j,k,ii,jj,kk)=X0(1)+(fac(1)**(0.5*(1.+X(1,i,j,K,ii,jj,kk)))-1.)*DXh(1) + END DO; END DO; END DO + END DO; END DO; END DO ELSE DXh(1)=0.5*DX(1) DO kk=1,nElems(3); DO jj=1,nElems(2); DO ii=1,nElems(1) ! loop over Chebyshev points DO k=0,Nin;DO j=0,Nin;DO i=0,Nin ! overwrite r,s,t with x,y,z - X(1,i,j,k,ii,jj,kk)=X0(1)+(1.+X(1,i,j,K,ii,jj,kk))*DXh(1) - END DO; END DO; END DO - END DO; END DO; END DO + X(1,i,j,k,ii,jj,kk)=X0(1)+(1.+X(1,i,j,K,ii,jj,kk))*DXh(1) + END DO; END DO; END DO + END DO; END DO; END DO END IF IF(fac(2).GT.1.) THEN DXh(2)=DX(2)/(fac(2)-1.) @@ -671,18 +952,18 @@ SUBROUTINE BuildPhysicalDomain(Nin,nElems,X) ! loop over Chebyshev points DO k=0,Nin;DO j=0,Nin;DO i=0,Nin ! overwrite r,s,t with x,y,z - X(2,i,j,k,ii,jj,kk)=X0(2)+(fac(2)**(0.5*(1.+X(2,i,j,K,ii,jj,kk)))-1.)*DXh(2) - END DO; END DO; END DO - END DO; END DO; END DO + X(2,i,j,k,ii,jj,kk)=X0(2)+(fac(2)**(0.5*(1.+X(2,i,j,K,ii,jj,kk)))-1.)*DXh(2) + END DO; END DO; END DO + END DO; END DO; END DO ELSE DXh(2)=0.5*DX(2) DO kk=1,nElems(3); DO jj=1,nElems(2); DO ii=1,nElems(1) ! loop over Chebyshev points DO k=0,Nin;DO j=0,Nin;DO i=0,Nin ! overwrite r,s,t with x,y,z - X(2,i,j,k,ii,jj,kk)=X0(2)+(1.+X(2,i,j,K,ii,jj,kk))*DXh(2) - END DO; END DO; END DO - END DO; END DO; END DO + X(2,i,j,k,ii,jj,kk)=X0(2)+(1.+X(2,i,j,K,ii,jj,kk))*DXh(2) + END DO; END DO; END DO + END DO; END DO; END DO END IF IF(fac(3).GT.1.) THEN DXh(3)=DX(3)/(fac(3)-1.) @@ -690,18 +971,18 @@ SUBROUTINE BuildPhysicalDomain(Nin,nElems,X) ! loop over Chebyshev points DO k=0,Nin;DO j=0,Nin;DO i=0,Nin ! overwrite r,s,t with x,y,z - X(3,i,j,k,ii,jj,kk)=X0(3)+(fac(3)**(0.5*(1.+X(3,i,j,K,ii,jj,kk)))-1.)*DXh(3) - END DO; END DO; END DO - END DO; END DO; END DO + X(3,i,j,k,ii,jj,kk)=X0(3)+(fac(3)**(0.5*(1.+X(3,i,j,K,ii,jj,kk)))-1.)*DXh(3) + END DO; END DO; END DO + END DO; END DO; END DO ELSE DXh(3)=0.5*DX(3) DO kk=1,nElems(3); DO jj=1,nElems(2); DO ii=1,nElems(1) ! loop over Chebyshev points DO k=0,Nin;DO j=0,Nin;DO i=0,Nin ! overwrite r,s,t with x,y,z - X(3,i,j,k,ii,jj,kk)=X0(3)+(1.+X(3,i,j,K,ii,jj,kk))*DXh(3) - END DO; END DO; END DO - END DO; END DO; END DO + X(3,i,j,k,ii,jj,kk)=X0(3)+(1.+X(3,i,j,K,ii,jj,kk))*DXh(3) + END DO; END DO; END DO + END DO; END DO; END DO END IF END SELECT END SUBROUTINE BuildPhysicalDomain @@ -727,13 +1008,13 @@ SUBROUTINE EvaluateFace(X,Para,DX,XP,WhichMapping,WhichFace) REAL :: PI,Z_Value,Fac,Para2(2) ! ? !=================================================================================================================================== SELECT CASE(WhichMapping) -CASE(1) ! TriLinearMapping +CASE(1) ! TriLinearMapping SELECT CASE(WhichFace) CASE(1) ! Face 1, Zeta=-1, Para = Xi,Eta - X = 0.25*( XP(:,1)*(1.-Para(1))*(1.-Para(2))& - +XP(:,2)*(1.+Para(1))*(1.-Para(2))& - +XP(:,3)*(1.+Para(1))*(1.+Para(2))& - +XP(:,4)*(1.-Para(1))*(1.+Para(2))) + X = 0.25*( XP(:,1)*(1.-Para(1))*(1.-Para(2))& + +XP(:,2)*(1.+Para(1))*(1.-Para(2))& + +XP(:,3)*(1.+Para(1))*(1.+Para(2))& + +XP(:,4)*(1.-Para(1))*(1.+Para(2))) CASE(2) ! Face 2, Eta=-1, Para = Xi,Zeta X = 0.25*( XP(:,1)*(1.-Para(1))*(1.-Para(2))& +XP(:,2)*(1.+Para(1))*(1.-Para(2))& @@ -768,39 +1049,39 @@ SUBROUTINE EvaluateFace(X,Para,DX,XP,WhichMapping,WhichFace) +XP(:,2)*(1.+Para(1))*(1.-Para(2))& +XP(:,3)*(1.+Para(1))*(1.+Para(2))& +XP(:,4)*(1.-Para(1))*(1.+Para(2)))& - +(/0.,0.,DX(3)/)*(1.-Para(1)**2)*(1.-Para(2)**2) + +(/0.,0.,DX(3)/)*(1.-Para(1)**2)*(1.-Para(2)**2) CASE(2) ! Face 2, Eta=-1, Para = Xi,Zeta X = 0.25*( XP(:,1)*(1.-Para(1))*(1.-Para(2))& +XP(:,2)*(1.+Para(1))*(1.-Para(2))& +XP(:,5)*(1.-Para(1))*(1.+Para(2))& +XP(:,6)*(1.+Para(1))*(1.+Para(2)))& - +(/0.,DX(2),0./)*(1.-Para(1)**2)*(1.-Para(2)**2) + +(/0.,DX(2),0./)*(1.-Para(1)**2)*(1.-Para(2)**2) CASE(3) ! Face 3, Xi=+1, Para = Eta, Zeta X = 0.25*( XP(:,2)*(1.-Para(1))*(1.-Para(2))& +XP(:,3)*(1.+Para(1))*(1.-Para(2))& +XP(:,6)*(1.-Para(1))*(1.+Para(2))& +XP(:,7)*(1.+Para(1))*(1.+Para(2)))& - +(/DX(1),0.,0./)*(1.-Para(1)**2)*(1.-Para(2)**2) + +(/DX(1),0.,0./)*(1.-Para(1)**2)*(1.-Para(2)**2) CASE(4) ! Face 4, Eta=+1, Para = Xi, Zeta X = 0.25*( XP(:,3)*(1.+Para(1))*(1.-Para(2))& +XP(:,4)*(1.-Para(1))*(1.-Para(2))& +XP(:,7)*(1.+Para(1))*(1.+Para(2))& +XP(:,8)*(1.-Para(1))*(1.+Para(2)))& - +(/0.,DX(2),0./)*(1.-Para(1)**2)*(1.-Para(2)**2) + +(/0.,DX(2),0./)*(1.-Para(1)**2)*(1.-Para(2)**2) CASE(5) ! Face 5, Xi=-1, Para = Eta, Zeta X = 0.25*( XP(:,1)*(1.-Para(1))*(1.-Para(2))& +XP(:,4)*(1.+Para(1))*(1.-Para(2))& +XP(:,5)*(1.-Para(1))*(1.+Para(2))& +XP(:,8)*(1.+Para(1))*(1.+Para(2)))& - +(/DX(1),0.,0./)*(1.-Para(1)**2)*(1.-Para(2)**2) + +(/DX(1),0.,0./)*(1.-Para(1)**2)*(1.-Para(2)**2) CASE(6) ! Face 6, Zeta=+1, Xi, Eta X = 0.25*( XP(:,5)*(1.-Para(1))*(1.-Para(2))& +XP(:,6)*(1.+Para(1))*(1.-Para(2))& +XP(:,7)*(1.+Para(1))*(1.+Para(2))& +XP(:,8)*(1.-Para(1))*(1.+Para(2)))& - +(/0.,0.,DX(3)/)*(1.-Para(1)**2)*(1.-Para(2)**2) + +(/0.,0.,DX(3)/)*(1.-Para(1)**2)*(1.-Para(2)**2) END SELECT CASE(3) ! half cylinder, Kopriva Book pg. 263 (we change y and z direction!! Changed to segment: PHI=180deg for half cylinder) PI=ACOS(-1.) @@ -814,7 +1095,7 @@ SUBROUTINE EvaluateFace(X,Para,DX,XP,WhichMapping,WhichFace) ! Scale azimuth of half cylinder Para2(2)= Fac*(Para(2)+1.)-1. ! Determine the values of the face and move it to y=-dy - ! GammaI + ! GammaI Gamma1 = -r_0 * COS(PI*(Para2(2) + 1.)*0.5)*e1 + r_0*sin(PI*(Para2(2) + 1.)*0.5)*e2 Gamma2 = (r_0 + (r_inf-r_0)*(Para2(1)+1.)*0.5)*e1 Gamma3 = -r_inf * COS(PI*(Para2(2) + 1.)*0.5)*e1 + r_inf*sin(PI*(Para2(2) + 1.)*0.5)*e2 @@ -824,7 +1105,7 @@ SUBROUTINE EvaluateFace(X,Para,DX,XP,WhichMapping,WhichFace) Corner2 = (r_0 + (r_inf-r_0)*(-1.+1.)*0.5)*e1 != Gamma2 mit Para = -1.: Para2=-1 Corner3 = -r_inf * COS(PI*((2.*Fac-1.) + 1.)*0.5)*e1 + r_inf*sin(PI*((2.*Fac-1.) + 1.)*0.5)*e2 != Gamma3 mit Para = +1.: Para2=2*Fac-1 Corner4 = (-r_0 - (r_inf-r_0)*((2.*Fac-1.)+1.)*0.5)*e1 != Gamma4 mit Para = +1.: Para2=2*Fac-1 - ! Now use the classic "2D" coons mapping to determine the face, Kopriva Book pg. 230 + ! Now use the classic "2D" coons mapping to determine the face, Kopriva Book pg. 230 X = 0.5 * ( Gamma4*(1.-Para2(1)) + Gamma2*(1.+Para2(1)) + Gamma1*(1.-Para2(2)) + Gamma3*(1.+Para2(2)))& - 0.25* ( (1.-Para2(1))*(Corner1*(1.-Para2(2))+Corner4*(1.+Para2(2)))+(1.+Para2(1))*(Corner2*(1.-Para2(2))+Corner3*(1.+Para2(2)))) X(3) = -dy @@ -858,7 +1139,7 @@ SUBROUTINE EvaluateFace(X,Para,DX,XP,WhichMapping,WhichFace) ! Scale azimuth of half cylinder Para2(2)= Fac*(Para(2)+1.)-1. ! Determine the values of the face and move it to y=+dy - ! GammaI + ! GammaI Gamma1 = -r_0 * COS(PI*(Para2(2) + 1.)*0.5)*e1 + r_0*sin(PI*(Para2(2) + 1.)*0.5)*e2 Gamma2 = (r_0 + (r_inf-r_0)*(Para2(1)+1.)*0.5)*e1 Gamma3 = -r_inf * COS(PI*(Para2(2) + 1.)*0.5)*e1 + r_inf*sin(PI*(Para2(2) + 1.)*0.5)*e2 @@ -868,7 +1149,7 @@ SUBROUTINE EvaluateFace(X,Para,DX,XP,WhichMapping,WhichFace) Corner2 = (r_0 + (r_inf-r_0)*(-1.+1.)*0.5)*e1 != Gamma2 mit Para = -1.: Para2=-1 Corner3 = -r_inf * COS(PI*((2.*Fac-1.) + 1.)*0.5)*e1 + r_inf*sin(PI*((2.*Fac-1.) + 1.)*0.5)*e2 != Gamma3 mit Para = +1.: Para2=2*Fac-1 Corner4 = (-r_0 - (r_inf-r_0)*((2.*Fac-1.)+1.)*0.5)*e1 != Gamma4 mit Para = +1.: Para2=2*Fac-1 - ! Now use the classic "2D" coons mapping to determine the face, Kopriva Book pg. 230 + ! Now use the classic "2D" coons mapping to determine the face, Kopriva Book pg. 230 X = 0.5 * ( Gamma4*(1.-Para2(1)) + Gamma2*(1.+Para2(1)) + Gamma1*(1.-Para2(2)) + Gamma3*(1.+Para2(2)))& - 0.25* ( (1.-Para2(1))*(Corner1*(1.-Para2(2))+Corner4*(1.+Para2(2)))+(1.+Para2(1))*(Corner2*(1.-Para2(2))+Corner3*(1.+Para2(2)))) X(3)=dy @@ -882,7 +1163,7 @@ SUBROUTINE EvaluateFace(X,Para,DX,XP,WhichMapping,WhichFace) SELECT CASE(WhichFace) CASE(1) ! Face 1, Zeta=-1, Para = Xi,Eta ! Determine the values of the face and move it to y=-dy - ! GammaI + ! GammaI Gamma1 = +r_0 * COS(PI*(Para(2) + 1.)*0.5)*e1 - r_0*sin(PI*(Para(2) + 1.)*0.5)*e2 Gamma2 = (r_0 + (r_inf-r_0)*(Para(1)+1.)*0.5)*e1 Gamma3 = +r_inf * COS(PI*(Para(2) + 1.)*0.5)*e1 - r_inf*sin(PI*(Para(2) + 1.)*0.5)*e2 @@ -892,7 +1173,7 @@ SUBROUTINE EvaluateFace(X,Para,DX,XP,WhichMapping,WhichFace) Corner2 = (/+r_0,0.,0./) != Gamma2 mit Para(2) = -1. = XP2 without y component Corner3 = (/+r_inf,0.,0./) != Gamma3 mit Para(1) = +1. = XP6 without y component Corner4 = (/+r_inf,0.,0./) != Gamma4 mit Para(2) = +1. = XP7 without y component - ! Now use the classic "2D" coons mapping to determine the face, Kopriva Book pg. 230 + ! Now use the classic "2D" coons mapping to determine the face, Kopriva Book pg. 230 X = 0.5 * ( Gamma4*(1.-Para(1)) + Gamma2*(1.+Para(1)) + Gamma1*(1.-Para(2)) + Gamma3*(1.+Para(2)))& - 0.25* ( (1.-Para(1))*(Corner1*(1.-Para(2))+Corner4*(1.+Para(2)))+(1.+Para(1))*(Corner2*(1.-Para(2))+Corner3*(1.+Para(2)))) X(3) = -dy @@ -917,7 +1198,7 @@ SUBROUTINE EvaluateFace(X,Para,DX,XP,WhichMapping,WhichFace) CASE(6) ! Face 6, Zeta=+1, Xi, Eta ! The same as Face 2, except that we have now +dy instead of -dy in y direction! ! Determine the values of the face and move it to y=+dy - ! GammaI + ! GammaI Gamma1 = r_0 * COS(PI*(Para(2) + 1.)*0.5)*e1 - r_0*sin(PI*(Para(2) + 1.)*0.5)*e2 Gamma2 = (r_0 + (r_inf-r_0)*(Para(1)+1.)*0.5)*e1 Gamma3 = r_inf * COS(PI*(Para(2) + 1.)*0.5)*e1 - r_inf*sin(PI*(Para(2) + 1.)*0.5)*e2 @@ -927,11 +1208,11 @@ SUBROUTINE EvaluateFace(X,Para,DX,XP,WhichMapping,WhichFace) Corner2 = (/r_0,0.,0./) != Gamma2 mit Para(2) = -1. Corner3 = (/r_inf,0.,0./) != Gamma3 mit Para(1) = +1. Corner4 = (/r_inf,0.,0./) != Gamma4 mit Para(2) = +1. - ! Now use the classic "2D" coons mapping to determine the face, Kopriva Book pg. 230 + ! Now use the classic "2D" coons mapping to determine the face, Kopriva Book pg. 230 X = 0.5 * ( Gamma4*(1.-Para(1)) + Gamma2*(1.+Para(1)) + Gamma1*(1.-Para(2)) + Gamma3*(1.+Para(2)))& - 0.25* ( (1.-Para(1))*(Corner1*(1.-Para(2))+Corner4*(1.+Para(2)))+(1.+Para(1))*(Corner2*(1.-Para(2))+Corner3*(1.+Para(2)))) X(3)=dy - + END SELECT !whichface CASE(5) !SINE BUMP @@ -949,42 +1230,42 @@ SUBROUTINE EvaluateFace(X,Para,DX,XP,WhichMapping,WhichFace) Z_Value=R_0*(1. - MIN(DX(1)/R_inf,1.)**2)**6 gamma1= DX(1)*Para(1)*e1 + R_0*(1. - MIN(ABS(DX(1)*Para(1))/R_inf,1.)**2)**6*e3 gamma2= DX(1)*e1+ (Z_Value*(1.-Para(2))*0.5 + DX(3)*(1+Para(2)))*e3 - gamma3= DX(1)*Para(1)*e1 + 2.*DX(3)*e3 + gamma3= DX(1)*Para(1)*e1 + 2.*DX(3)*e3 gamma4=-DX(1)*e1+ (Z_Value*(1.-Para(2))*0.5 + DX(3)*(1+Para(2)))*e3 ! Determine the 4 corner points of the face Corner1 = (/-DX(1),0.,Z_Value/) != Gamma1 mit Para(1) = -1. Corner2 = (/+DX(1),0.,Z_Value/) != Gamma2 mit Para(2) = -1. Corner3 = (/+DX(1),0.,2.*DX(3)/) != Gamma3 mit Para(1) = +1. Corner4 = (/-DX(1),0.,2.*DX(3)/) != Gamma4 mit Para(2) = +1. - ! Now use the classic "2D" coons mapping to determine the face, Kopriva Book pg. 230 + ! Now use the classic "2D" coons mapping to determine the face, Kopriva Book pg. 230 X = 0.5 * ( Gamma4*(1.-Para(1)) + Gamma2*(1.+Para(1)) + Gamma1*(1.-Para(2)) + Gamma3*(1.+Para(2)))& - 0.25* ( (1.-Para(1))*(Corner1*(1.-Para(2))+Corner4*(1.+Para(2)))+(1.+Para(1))*(Corner2*(1.-Para(2))+Corner3*(1.+Para(2)))) X(2)=DX(2) CASE(3) ! Face 3, Xi=+1, Para = Eta, Zeta Z_Value=R_0*(1. - MIN(DX(1)/R_inf,1.)**2)**6 X(1)=DX(1) - X(2)=DX(2)*Para(1) + X(2)=DX(2)*Para(1) X(3)=Z_Value*(1.-Para(2))*0.5+DX(3)*(1.+Para(2)) CASE(4) ! Face 4, Eta=+1, Para = Xi, Zeta !gamma1= DX(1)*Para(1)*e1 + R_0*0.5*(1.+cos(MIN(ABS(DX(1)*Para(1))/R_inf,1.)*PI))*e3 Z_Value=R_0*(1. - MIN(DX(1)/R_inf,1.)**2)**6 gamma1= DX(1)*Para(1)*e1 + R_0*(1. - MIN(ABS(DX(1)*Para(1))/R_inf,1.)**2)**6*e3 gamma2= DX(1)*e1+ (Z_Value*(1.-Para(2))*0.5 + DX(3)*(1+Para(2)))*e3 - gamma3= DX(1)*Para(1)*e1 + 2.*DX(3)*e3 + gamma3= DX(1)*Para(1)*e1 + 2.*DX(3)*e3 gamma4=-DX(1)*e1+ (Z_Value*(1.-Para(2))*0.5 + DX(3)*(1+Para(2)))*e3 ! Determine the 4 corner points of the face Corner1 = (/-DX(1),0.,Z_Value/) != Gamma1 mit Para(1) = -1. Corner2 = (/+DX(1),0.,Z_Value/) != Gamma2 mit Para(2) = -1. Corner3 = (/+DX(1),0.,2.*DX(3)/) != Gamma3 mit Para(1) = +1. Corner4 = (/-DX(1),0.,2.*DX(3)/) != Gamma4 mit Para(2) = +1. - ! Now use the classic "2D" coons mapping to determine the face, Kopriva Book pg. 230 + ! Now use the classic "2D" coons mapping to determine the face, Kopriva Book pg. 230 X = 0.5 * ( Gamma4*(1.-Para(1)) + Gamma2*(1.+Para(1)) + Gamma1*(1.-Para(2)) + Gamma3*(1.+Para(2)))& - 0.25* ( (1.-Para(1))*(Corner1*(1.-Para(2))+Corner4*(1.+Para(2)))+(1.+Para(1))*(Corner2*(1.-Para(2))+Corner3*(1.+Para(2)))) X(2)=-DX(2) CASE(5) ! Face 5, Xi=-1, Para = Eta, Zeta Z_Value=R_0*(1. - MIN(DX(1)/R_inf,1.)**2)**6 X(1)=-DX(1) - X(2)=DX(2)*Para(1) + X(2)=DX(2)*Para(1) X(3)=Z_Value*(1.-Para(2))*0.5+DX(3)*(1.+Para(2)) CASE(6) ! Face 6, Zeta=+1, Xi, Eta X(1:2)=DX(1:2)*Para(1:2) @@ -1013,7 +1294,7 @@ SUBROUTINE EvaluateFace(X,Para,DX,XP,WhichMapping,WhichFace) Corner2 = (/1.,0.,1./) != Gamma2 mit Para(2) = -1. Corner3 = (/2./3.,0.,7./6./) != Gamma3 mit Para(1) = +1. Corner4 = (/0.,0.,0.5/) != Gamma4 mit Para(2) = +1. - ! Now use the classic "2D" coons mapping to determine the face, Kopriva Book pg. 230 + ! Now use the classic "2D" coons mapping to determine the face, Kopriva Book pg. 230 X = 0.5 * ( Gamma4*(1.-Para(1)) + Gamma2*(1.+Para(1)) + Gamma1*(1.-Para(2)) + Gamma3*(1.+Para(2)))& - 0.25* ( (1.-Para(1))*(Corner1*(1.-Para(2))+Corner4*(1.+Para(2)))+(1.+Para(1))*(Corner2*(1.-Para(2))+Corner3*(1.+Para(2)))) X(2)=0.25 @@ -1031,7 +1312,7 @@ SUBROUTINE EvaluateFace(X,Para,DX,XP,WhichMapping,WhichFace) Corner2 = (/1.,0.,1./) != Gamma2 mit Para(2) = -1. Corner3 = (/2./3.,0.,7./6./) != Gamma3 mit Para(1) = +1. Corner4 = (/0.,0.,0.5/) != Gamma4 mit Para(2) = +1. - ! Now use the classic "2D" coons mapping to determine the face, Kopriva Book pg. 230 + ! Now use the classic "2D" coons mapping to determine the face, Kopriva Book pg. 230 X = 0.5 * ( Gamma4*(1.-Para(1)) + Gamma2*(1.+Para(1)) + Gamma1*(1.-Para(2)) + Gamma3*(1.+Para(2)))& - 0.25* ( (1.-Para(1))*(Corner1*(1.-Para(2))+Corner4*(1.+Para(2)))+(1.+Para(1))*(Corner2*(1.-Para(2))+Corner3*(1.+Para(2)))) X(2)=-0.25 diff --git a/src/mesh/mesh.f90 b/src/mesh/mesh.f90 index 8361d47..fe72899 100644 --- a/src/mesh/mesh.f90 +++ b/src/mesh/mesh.f90 @@ -650,7 +650,6 @@ SUBROUTINE fillMesh() IF(useCurveds.AND.Logging) CALL CountSplines() ! In case of restart there can be splines END IF CALL buildEdges() -IF(generateFEMconnectivity) CALL buildFEMconnectivity() ! check if sides to be curved exist curvedFound=.FALSE. @@ -734,7 +733,7 @@ SUBROUTINE fillMesh() IF(mortarFound) EXIT !do loop END DO !iElem -IF(mortarFound.AND.generateFEMconnectivity) CALL abort(__STAMP__,"generate FEM connectivity not yet implemented for mortar meshes!") +! IF(mortarFound.AND.generateFEMconnectivity) CALL abort(__STAMP__,"generate FEM connectivity not yet implemented for mortar meshes!") IF(doExactSurfProjection) CALL ProjectToExactSurfaces() ! get element types @@ -787,6 +786,7 @@ SUBROUTINE fillMesh() !after rebuild , mortars should be fine, but checking is better: CALL CheckMortarWaterTight() END IF !mortarFound +IF(generateFEMconnectivity) CALL buildFEMconnectivity() ! apply meshscale before output (default) IF( (doShift.AND.postShift) .OR. (doScale.AND.postScale) ) THEN diff --git a/src/mesh/mesh_postdeform.f90 b/src/mesh/mesh_postdeform.f90 index dfa693d..f166ee6 100644 --- a/src/mesh/mesh_postdeform.f90 +++ b/src/mesh/mesh_postdeform.f90 @@ -52,7 +52,7 @@ SUBROUTINE PostDeform() USE MOD_Globals USE MOD_Mesh_Vars ,ONLY: tElem,Elems,MeshPostDeform,PostDeform_useGL USE MOD_Mesh_Vars ,ONLY: N,nMeshElems -USE MOD_Basis_Vars ,ONLY: HexaMap +USE MOD_Basis_Vars ,ONLY: HexaMap,PrismMap,PyraMap USE MOD_Basis1D ,ONLY: LegGaussLobNodesAndWeights,BarycentricWeights,InitializeVandermonde USE MOD_ChangeBasis ,ONLY: ChangeBasis3D !MODULE OUTPUT VARIABLES @@ -75,6 +75,10 @@ SUBROUTINE PostDeform() REAL,DIMENSION(0:N,0:N) :: Vdm_EQtoGL, Vdm_GLtoEQ REAL :: xElem( 3,0:N,0:N,0:N,nMeshElems) REAL :: xElemDeform(3,0:N,0:N,0:N,nMeshElems) +REAL :: xElem1D( 3,1:NINT((N+1)**2*(N+2)/2.),nMeshElems) +REAL :: xElemDeform1D( 3,1:NINT((N+1)**2*(N+2)/2.),nMeshElems) +!REAL :: xElem1D( 3,1:(N+1)**3,nMeshElems) +!REAL :: xElemDeform1D( 3,1:(N+1)**3,nMeshElems) INTEGER :: HexaMapN1(8,3) !=================================================================================================================================== IF(MeshPostDeform.EQ.0) RETURN @@ -100,6 +104,7 @@ SUBROUTINE PostDeform() CALL InitializeVandermonde(N,N,wBary_GL,xi_GL,xi_EQ,Vdm_GLtoEQ) xElem=0. +xElem1D=0. !Copy Equidist. element nodes DO iElem=1,nMeshElems aElem=>Elems(iElem)%ep @@ -111,27 +116,32 @@ SUBROUTINE PostDeform() END DO !iNode ELSE DO iNode=1,aElem%nCurvedNodes - ijk(:)=HexaMap(iNode,:) - xElem(1:3,ijk(1),ijk(2),ijk(3),aElem%ind)= aElem%CurvedNode(iNode)%np%x(:) + !ijk(:)=HexaMap(iNode,:) + !ijk(:)=PrismMap(iNode,:) + ijk(:)=PyraMap(iNode,:) + xElem1D(1:3,iNode,aElem%ind)= aElem%CurvedNode(iNode)%np%x(:) + !IF(xElem1D(3,iNode,aElem%ind).EQ.1.) print *,xElem1D(1:3,iNode,aElem%ind) aElem%CurvedNode(iNode)%np%tmp=-1 END DO !iNode END IF !N=1 END DO ! iElem !transform Equidist. to Gauss-Lobatto points -IF((PostDeform_useGL).AND.(N.GT.2))THEN - CALL ChangeBasis3D(3,nMeshElems,N,N,Vdm_EQtoGL,xElem,xElem,.FALSE.) -END IF !PostDeform_useGL +!IF((PostDeform_useGL).AND.(N.GT.2))THEN + !CALL ChangeBasis3D(3,nMeshElems,N,N,Vdm_EQtoGL,xElem,xElem,.FALSE.) +!END IF !PostDeform_useGL !transform (all nodes are marked from -2 to -1) -nTotal=(N+1)**3*nMeshElems -CALL PostDeformFunc(nTotal,xElem,xElemDeform) -xElem = xElemDeform +!nTotal=(N+1)**3*nMeshElems +!nTotal=NINT((N+1)**2*(N+2)*0.5)*nMeshElems +nTotal=NINT((N+1)*(N+2)*(2*N+3)/6.)*nMeshElems +CALL PostDeformFunc(nTotal,xElem1D,xElemDeform1D) +xElem1D = xElemDeform1D -IF((PostDeform_useGL).AND.(N.GT.2))THEN - !transform back from GL to EQ - CALL ChangeBasis3D(3,nMeshElems,N,N,Vdm_GLtoEQ,xElem,xElem,.FALSE.) -END IF +!IF((PostDeform_useGL).AND.(N.GT.2))THEN + !!transform back from GL to EQ + !CALL ChangeBasis3D(3,nMeshElems,N,N,Vdm_GLtoEQ,xElem,xElem,.FALSE.) +!END IF ! copy back (all nodes are marked from -1 to 0) DO iElem=1,nMeshElems @@ -147,8 +157,11 @@ SUBROUTINE PostDeform() ELSE DO iNode=1,aElem%nCurvedNodes IF(aElem%CurvedNode(iNode)%np%tmp.EQ.-1)THEN - ijk(:)=HexaMap(iNode,:) - aElem%CurvedNode(iNode)%np%x(:)= xElem(1:3,ijk(1),ijk(2),ijk(3),aElem%ind) + !ijk(:)=HexaMap(iNode,:) + !ijk(:)=PrismMap(iNode,:) + ijk(:)=PyraMap(iNode,:) + !aElem%CurvedNode(iNode)%np%x(:)= xElem(1:3,ijk(1),ijk(2),ijk(3),aElem%ind) + aElem%CurvedNode(iNode)%np%x(:)= xElem1D(1:3,iNode,aElem%ind) aElem%CurvedNode(iNode)%np%tmp=0 END IF END DO !iNode @@ -198,6 +211,7 @@ SUBROUTINE PostDeformFunc(nTotal,X_in,X_out) CASE(1,11,12) DO i=1,nTotal x(:)=X_in(:,i) + IF(x(3).EQ.1) print *, x ! 2D box, x,y in [-1,1]^2, to cylinder with radius PostDeform_R0 (with PostDeform_Rtorus>0 to a torus, with zperiodic [0,1]) ! all points outside [-1,1]^2 will be mapped directly to a circle (p.e. 2,2 => sqrt(0.5)*PostDeform_R0*(2,2) ) ! inside [-1,1]^2 and outside [-0.5,0.5]^2 there will be a blending from a circle to a square diff --git a/src/mesh/mesh_tools.f90 b/src/mesh/mesh_tools.f90 index 8747131..449a872 100644 --- a/src/mesh/mesh_tools.f90 +++ b/src/mesh/mesh_tools.f90 @@ -12,7 +12,7 @@ ! Copyright (C) 2017 Claus-Dieter Munz ! This file is part of HOPR, a software for the generation of high-order meshes. ! -! HOPR is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License +! HOPR is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License ! as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. ! ! HOPR is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty @@ -30,7 +30,7 @@ MODULE MOD_Mesh_Tools IMPLICIT NONE PRIVATE !----------------------------------------------------------------------------------------------------------------------------------- -! GLOBAL VARIABLES +! GLOBAL VARIABLES !----------------------------------------------------------------------------------------------------------------------------------- ! Private Part --------------------------------------------------------------------------------------------------------------------- ! Public Part ---------------------------------------------------------------------------------------------------------------------- @@ -85,13 +85,13 @@ SUBROUTINE CountSplines() ! IMPLICIT VARIABLE HANDLING IMPLICIT NONE !----------------------------------------------------------------------------------------------------------------------------------- -! LOCAL VARIABLES +! LOCAL VARIABLES TYPE(tElem),POINTER :: Elem ! ? TYPE(tSide),POINTER :: Side ! ? -INTEGER :: splineCounter(3) ! 1: Boundary splines, 2: periodic splines, 3: +INTEGER :: splineCounter(3) ! 1: Boundary splines, 2: periodic splines, 3: !----------------------------------------------------------------------------------------------------------------------------------- ! INPUT VARIABLES -INTEGER, SAVE :: CallCount=0 ! Counter for calls of this subroutine, used for +INTEGER, SAVE :: CallCount=0 ! Counter for calls of this subroutine, used for !----------------------------------------------------------------------------------------------------------------------------------- ! OUTPUT VARIABLES LOGICAL :: hasChanged ! ? @@ -161,7 +161,7 @@ SUBROUTINE NetVisu() !----------------------------------------------------------------------------------------------------------------------------------- ! OUTPUT VARIABLES !----------------------------------------------------------------------------------------------------------------------------------- -! LOCAL VARIABLES +! LOCAL VARIABLES TYPE(tElem),POINTER :: Elem ! ? REAL :: BaryScale,BaryCoords(3) ! ? INTEGER :: nVal,iNode,iElem,nElems ! ? @@ -241,10 +241,10 @@ SUBROUTINE FEMNetVisu() !----------------------------------------------------------------------------------------------------------------------------------- ! OUTPUT VARIABLES !----------------------------------------------------------------------------------------------------------------------------------- - ! LOCAL VARIABLES - TYPE(tElem),POINTER :: Elem + ! LOCAL VARIABLES + TYPE(tElem),POINTER :: Elem TYPE(tSide),POINTER :: Side - TYPE(tLocalEdge),POINTER :: lEdge + TYPE(tLocalEdge),POINTER :: lEdge TYPE(tNode),POINTER :: aNode,bNode TYPE(tVertex),POINTER :: avert,bVert INTEGER :: nVal,iNode,nElems,nSidesElem,nSides,nEdges,iSide,iElemSide,iEdge,iLocEdge ! ? @@ -274,19 +274,19 @@ SUBROUTINE FEMNetVisu() nVal=2 VarNames(1)='elemind' VarNames(2)='UniqueFaceID' - + NodeMap=0 !mapping from Side nNodes to i,j [0;1] NodeMap(:,3)=(/1,2,3,3/) !tri NodeMap(:,4)=(/1,2,4,3/) !quad - + ALLOCATE(Coord(3,1:4,nSides)) ALLOCATE(Solution(nVal,1:4,nSides)) - + iSide=0 Elem=>firstElem DO WHILE(ASSOCIATED(elem)) - + nSidesElem=nSides_from_nNodes(Elem%nNodes) Side=>Elem%firstSide DO iElemSide=1,nSidesElem @@ -301,7 +301,7 @@ SUBROUTINE FEMNetVisu() Elem=>elem%nextElem END DO CALL Visualize(2,nVal,1,nSides,VarNames(1:nVal),Coord,Solution,FileString) - + DEALLOCATE(Coord,Solution) filestring=TRIM(ProjectName)//'_'//'Debugmesh_Edges' @@ -314,7 +314,7 @@ SUBROUTINE FEMNetVisu() VarNames(4)='EdgeIsMaster' VarNames(5)='VertexIsMaster' - + iEdge=0 Elem=>firstElem DO WHILE(ASSOCIATED(elem)) @@ -338,10 +338,10 @@ SUBROUTINE FEMNetVisu() Elem=>elem%nextElem END DO CALL Visualize(1,nVal,1,nEdges,VarNames(1:nVal),Coord,Solution,FileString) - + DEALLOCATE(Coord,Solution) CALL Timer(.FALSE.) - + END SUBROUTINE FEMnetVisu SUBROUTINE BCVisu() @@ -360,7 +360,7 @@ SUBROUTINE BCVisu() !----------------------------------------------------------------------------------------------------------------------------------- ! OUTPUT VARIABLES !----------------------------------------------------------------------------------------------------------------------------------- -! LOCAL VARIABLES +! LOCAL VARIABLES TYPE(tElem),POINTER :: Elem ! ? TYPE(tSide),POINTER :: Side ! ? REAL :: BaryScale,BaryCoords(3) ! ? @@ -592,7 +592,8 @@ SUBROUTINE chkSpl_Vol() USE MOD_Mesh_Vars ,ONLY:FirstElem USE MOD_Mesh_Vars ,ONLY:N USE MOD_Basis_Vars,ONLY:Vdm_Visu_Hexa,D_Visu_Hexa -USE MOD_Basis_Vars,ONLY:VisuHexaMapInv +USE MOD_Basis_Vars,ONLY:VdM_visu_Prism,D_visu_Prism +USE MOD_Basis_Vars,ONLY:VisuHexaMapInv,VisuPrismMapInv USE MOD_Basis_Vars,ONLY:nVisu USE MOD_Output ,ONLY:Visualize USE MOD_Output_vars,ONLY:Visu_sJ_limit @@ -640,14 +641,16 @@ SUBROUTINE chkSpl_Vol() Nplot=nVisu !number of equidistant points for visualization mesh Nplot_p1=Nplot+1 -Nplot_p1_3=Nplot_p1*Nplot_p1*Nplot_p1 +!Nplot_p1_3=(Nplot_p1)**3 +Nplot_p1_3=NINT((Nplot_p1)**2*(Nplot_p1+1)/2.) -ALLOCATE(xNode((N+1)**3,3)) -ALLOCATE(x( Nplot_p1_3,3)) -ALLOCATE(xt1(Nplot_p1_3,3)) -ALLOCATE(xt2(Nplot_p1_3,3)) -ALLOCATE(xt3(Nplot_p1_3,3)) -ALLOCATE(Jac(Nplot_p1_3)) +!ALLOCATE(xNode((N+1)**3,3)) +ALLOCATE(xNode(NINT((N+1)**2*(N+2)/2.),3)) +ALLOCATE(x (Nplot_p1_3,3)) +ALLOCATE(xt1 (Nplot_p1_3,3)) +ALLOCATE(xt2 (Nplot_p1_3,3)) +ALLOCATE(xt3 (Nplot_p1_3,3)) +ALLOCATE(Jac (Nplot_p1_3)) nVal=4 ALLOCATE(VarNames(nVal)) @@ -664,6 +667,40 @@ SUBROUTINE chkSpl_Vol() DO WHILE(ASSOCIATED(Elem)) IF(ASSOCIATED(Elem%CurvedNode)) THEN SELECT CASE(Elem%nNodes) + CASE(6) + iElem=iElem+1 + nNodes=Elem%nCurvedNodes + x=0. + xt1=0. + xt2=0. + xNode=0. + DO iNode=1,nNodes + IF(ASSOCIATED(Elem%curvedNode(iNode)%np))THEN + xNode(iNode,:)=Elem%curvedNode(iNode)%np%x + ELSE + CALL abort(__STAMP__,'Curved node array has not the right size.',999,999.) ! check required due to intel compiler error (02) + END IF + END DO + x = MATMUL(Vdm_Visu_Prism,xNode(1:nNodes,:)) + xt1= MATMUL(D_Visu_Prism(:,:,1),xNode(1:nNodes,:)) + xt2= MATMUL(D_Visu_Prism(:,:,2),xNode(1:nNodes,:)) + xt3= MATMUL(D_Visu_Prism(:,:,3),xNode(1:nNodes,:)) + DO iNode=1,Nplot_p1_3 + Jac(iNode)=SUM(xt1(iNode,:)*CROSS(xt2(iNode,:),xt3(iNode,:))) + END DO + sMaxJac=1./MAXVAL(Jac) + DO k=0,Nplot + DO j=0,Nplot + DO i=0,Nplot-j + l=VisuPrismMapInv(i,j,k) + Coord(:,l,iElem)=x(l,:) + Values(2,l,iElem)=Jac(l) + Values(3,l,iElem)=Jac(l)*sMaxJac + END DO + END DO + END DO + Values(1,:,iElem)=Elem%ind + Values(4,:,iElem)=MINVAL(Values(3,:,iElem)) CASE(8) iElem=iElem+1 nNodes=Elem%nCurvedNodes @@ -745,7 +782,7 @@ SUBROUTINE SetTempMarker(Elem,value,whichMarker) ! IMPLICIT VARIABLE HANDLING IMPLICIT NONE !----------------------------------------------------------------------------------------------------------------------------------- -! LOCAL VARIABLES +! LOCAL VARIABLES TYPE(tElem),POINTER,INTENT(IN) :: Elem ! ? INTEGER,INTENT(IN) :: value ! value to set tmp to LOGICAL,INTENT(IN),OPTIONAL :: whichMarker(8) ! 1/2: elem corner/curved, 3/4: side corner/curved, 5/6: edge corner/curved @@ -813,7 +850,7 @@ END SUBROUTINE SetTempMarker SUBROUTINE checkMortarWatertight() !=================================================================================================================================== ! Checks if surface normals of mortars are defined such that the mesh is freestream preserving =^ watertight -! builds a 1D basis to change equidistant -> gauss points (0:N_GP) and then use tensor-product gauss +! builds a 1D basis to change equidistant -> gauss points (0:N_GP) and then use tensor-product gauss ! for differentiation and integration. n_GP=N should be exact, since normal vector is of degree (2*N-1 ,2*N-1) ! since its a dot product of two polynomials of degree (N-1,N) * (N,N-1) !=================================================================================================================================== @@ -894,17 +931,17 @@ SUBROUTINE checkMortarWatertight() NsurfBig=EvalNsurf(XgeoSide) NsurfSmall=0. - DO p=1,Side%nMortars + DO p=1,Side%nMortars CALL PackGeo(N,Side%MortarSide(p)%sp,XgeoSide) NsurfSmall(:,p)=EvalNsurf(XgeoSide) - END DO + END DO NsurfErr= ABS(NsurfBig(1)-SUM(NsurfSmall(1,:))) & +ABS(NsurfBig(2)-SUM(NsurfSmall(2,:))) & +ABS(NsurfBig(3)-SUM(NsurfSmall(3,:))) IF(NsurfErr.GT.1.0E-12) THEN ERRWRITE(*,*) & '================> Mortar is not watertight, ERROR=',NsurfErr,' >1.0E-12, big side corners:' - ERRWRITE(*,*)'Nsurf: ', NsurfBig + ERRWRITE(*,*)'Nsurf: ', NsurfBig ERRWRITE(*,*)' P1: ', Side%OrientedNode(1)%np%x ERRWRITE(*,*)' P2: ', Side%OrientedNode(2)%np%x ERRWRITE(*,*)' P3: ', Side%OrientedNode(3)%np%x @@ -918,7 +955,7 @@ SUBROUTINE checkMortarWatertight() ERRWRITE(*,*)' P2: ', Side%MortarSide(p)%sp%OrientedNode(2)%np%x ERRWRITE(*,*)' P3: ', Side%MortarSide(p)%sp%OrientedNode(3)%np%x ERRWRITE(*,*)' P4: ', Side%MortarSide(p)%sp%OrientedNode(4)%np%x - END DO + END DO WaterTight=WaterTight+1 maxNsurfErr=max(maxNsurfErr,NsurfErr) END IF ! diffsurf>0 @@ -960,7 +997,7 @@ SUBROUTINE checkMortarWatertight() pf=N-q; qf=N-p; CASE(4) pf=p; qf=N-q; - END SELECT + END SELECT XCheck=XGeoSide(:,p,q) + vec @@ -1036,8 +1073,8 @@ FUNCTION EvalNsurf(Xgeo_in) RESULT(Nsurf) END DO nVec=CROSS(dXdxiGP(:),dXdetaGP(:)) Nsurf(:)=Nsurf(:)+wGP(i)*wGP(j)*nVec(:) - END DO !i - END DO !j + END DO !i + END DO !j END FUNCTION EvalNsurf END SUBROUTINE checkMortarWaterTight diff --git a/src/mesh/mesh_vars.f90 b/src/mesh/mesh_vars.f90 index 8e8d2be..9ed7b79 100644 --- a/src/mesh/mesh_vars.f90 +++ b/src/mesh/mesh_vars.f90 @@ -82,6 +82,7 @@ MODULE MOD_Mesh_Vars INTEGER :: nEdges ! total number of edges for one Elem INTEGER :: nCurvedNodes ! Used for writing curveds to hdf5 mesh format INTEGER :: ind ! unique Element index for each element on all processors + INTEGER :: splitind ! unique Element index for each element on all processors INTEGER :: tmp ! CHARACTER --------------------------------------------------------! ! LOGICAL ----------------------------------------------------------!