Add high-order simplex elements; TODO: tetras

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

src/basis/basis_vars.f90

@@ -12,7 +12,7 @@
 ! Copyright (C) 2017 Claus-Dieter Munz <[email protected]>
 ! 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

src/basis/prismBasis.f90

@@ -12,7 +12,7 @@
 ! Copyright (C) 2017 Claus-Dieter Munz <[email protected]>
 ! 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,28 +46,28 @@ 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
 ! IMPLICIT VARIABLE HANDLING
 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,34 +293,34 @@ 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)
   CALL GradJacobiP(nNodes,a,     0, 0,i,dPa)
   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