[feature.FEM.EdgeAndCornerConnectivity] New Feature: Edge and Vertex connectivity, needed for a FEM solver

CMakePresets.json

@@ -0,0 +1,38 @@
+{
+    "version": 3,
+    "cmakeMinimumRequired": {
+        "major": 3,
+        "minor": 22,
+        "patch": 0
+    },
+    "configurePresets": [
+        {
+            "name": "hopr_config_debug_nohdf5_nocgns",
+            "displayName": "hopr configure: default debug build no cgns no hdf5",
+            "hidden": false,
+            "cacheVariables": {
+                "CMAKE_BUILD_TYPE": "Debug",
+                "LIBS_BUILD_HDF5": "Off",
+                "LIBS_USE_CGNS": "Off"
+            }
+        },
+        {
+            "name": "hopr_config_release_nohdf5_nocgns",
+            "displayName": "hopr configure: default release build no cgns no hdf5",
+            "hidden": false,
+            "cacheVariables": {
+                "CMAKE_BUILD_TYPE": "Release",
+                "LIBS_BUILD_HDF5": "Off",
+                "LIBS_USE_CGNS": "Off"
+            }
+        }
+    ],
+    "testPresets": [
+        {
+            "name": "HOPR_ctest",
+            "displayName": "Test HOPR running ctest ",
+            "configurePreset": "HOPR_config_release",
+            "output": {"outputOnFailure": true}
+        }
+    ]
+}

src/basis/basis.f90

@@ -78,8 +78,7 @@ SUBROUTINE InitBasis()
 WRITE(UNIT_StdOut,'(132("-"))')
 WRITE(UNIT_stdOut,'(A)') ' INIT BASIS...'
 IF(.NOT.MeshInitDone)THEN
-  CALL abort(__STAMP__,  &
-             'ERROR: InitMesh has to be called before InitBasis!',999,999.)
+  CALL abort(__STAMP__,'ERROR: InitMesh has to be called before InitBasis!',999,999.)
 END IF
 WRITE(tmpstr,'(I4)')N
 nVisu=GETINT('nVisu',tmpstr)
@@ -284,8 +283,7 @@ SUBROUTINE GetNodesAndWeights(N_in,NodeType_in,xIP,wIP,wIPBary)
     ! first order intergration !!!
     wIP=2./REAL(N_in+1)
   CASE DEFAULT
-    CALL Abort(__STAMP__,&
-      'NodeType "'//TRIM(NodeType_in)//'" in GetNodesAndWeights not found!')
+    CALL Abort(__STAMP__,'NodeType "'//TRIM(NodeType_in)//'" in GetNodesAndWeights not found!')
   END SELECT
 ELSE
   SELECT CASE(TRIM(NodeType_in))
@@ -306,8 +304,7 @@ SUBROUTINE GetNodesAndWeights(N_in,NodeType_in,xIP,wIP,wIPBary)
       xIP(i) = 1./REAL(N_in+1)+2.*REAL(i)/REAL(N_in+1) - 1.
     END DO
   CASE DEFAULT
-    CALL Abort(__STAMP__,&
-      'NodeType "'//TRIM(NodeType_in)//'" in GetNodesAndWeights not found!')
+    CALL Abort(__STAMP__,'NodeType "'//TRIM(NodeType_in)//'" in GetNodesAndWeights not found!')
   END SELECT
 END IF !present wIP
 IF(PRESENT(wIPBary)) CALL BarycentricWeights(N_in,xIP,wIPBary)

src/commandlinearguments.f90

@@ -83,14 +83,14 @@ SUBROUTINE ParseCommandlineArguments()
     ! Print version and exit
     WRITE(UNIT_stdOut,'(A)')"hopr version "&
         //TRIM(int2strf(MajorVersion))//"."//TRIM(int2strf(MinorVersion))//"."//TRIM(int2strf(PatchVersion))
-    STOP
+    STOP 0
   CASE DEFAULT ! 2+
     ! Print extended version and exit
     WRITE(UNIT_stdOut,'(A)')"hopr version "&
         //TRIM(int2strf(MajorVersion))//"."//TRIM(int2strf(MinorVersion))//"."//TRIM(int2strf(PatchVersion))&
         //" ("//TRIM(GIT_CURRENT_COMMIT)//", "//TRIM(BUILD_DATE)       //")"                                &
         //" ["//TRIM(BUILD_VERSION_GCC) //", "//TRIM(BUILD_VERSION_MPI)//"]"
-    STOP
+    STOP 0
 END SELECT
 
 ! Get all remaining parameters
@@ -163,7 +163,7 @@ SUBROUTINE PrintHelp()
 WRITE(UNIT_stdOut,'(A)') '  -V, --version    display the version number and exit'
 WRITE(UNIT_stdOut,'(A)') '                   when given twice, print more information about the build'
 ! CALL Abort(__STAMP__,'Parameter file not specified!')
-STOP
+STOP 0
 END SUBROUTINE PrintHelp
 
 END MODULE MOD_Commandline_Arguments

src/globals.f90

@@ -9,6 +9,8 @@
 ! /____//   /____//  /______________//  /____//           /____//   |_____/)    ,X`      XXX`
 ! )____)    )____)   )______________)   )____)            )____)    )_____)   ,xX`     .XX`
 !                                                                           xxX`      XXx
+! Copyright (C) 2023  Florian Hindenlang <hindenlang@gmail.com>
+! Copyright (C) 2023  Tobias Ott <tobias.ott@proton.me>
 ! Copyright (C) 2017 Claus-Dieter Munz <munz@iag.uni-stuttgart.de>
 ! This file is part of HOPR, a software for the generation of high-order meshes.
 !
@@ -55,8 +57,8 @@ MODULE MOD_Globals
 LOGICAL                     :: Logging                    ! Set .TRUE. to activate logging function for each processor
 
 INTEGER,PARAMETER           :: MajorVersion = 1           !> HoprVersion saved in each hdf5 file with hdf5 header
-INTEGER,PARAMETER           :: MinorVersion = 2           !> HoprVersion saved in each hdf5 file with hdf5 header
-INTEGER,PARAMETER           :: PatchVersion = 1           !> HoprVersion saved in each hdf5 file with hdf5 header
+INTEGER,PARAMETER           :: MinorVersion = 3           !> HoprVersion saved in each hdf5 file with hdf5 header
+INTEGER,PARAMETER           :: PatchVersion = 0           !> HoprVersion saved in each hdf5 file with hdf5 header
 INTEGER,PARAMETER           :: HoprVersionInt = PatchVersion+MinorVersion*100+MajorVersion*10000 !> Hopr version number saved in each hdf5 file with hdf5 header
 CHARACTER(LEN=10)           :: HoprVersionStr             !> Hopr version string saved in each hdf5 file with hdf5 header
 
@@ -405,8 +407,7 @@ SUBROUTINE SolveLinSys(dim1,nRHS,A,B)
 !-----------------------------------------------------------------------------------------------------------------------------------
 !-----------------------------------------------------------------------------------------------------------------------------------
 ! LOCAL VARIABLES
-INTEGER            :: IPIV(dim1),INFO,lwork ! ?
-REAL               :: WORK(dim1*dim1)  ! ?
+INTEGER            :: IPIV(dim1),INFO  ! ?
 !===================================================================================================================================
 ! DGETRF computes an LU factorization of a general M-by-N matrix A
 ! using partial pivoting with row interchanges. A will be overwritten with LU factorization.

src/hopr.f90

@@ -33,7 +33,7 @@ PROGRAM HOPR
 USE MOD_Mesh,                  ONLY: InitMesh,FillMesh
 USE MOD_Mesh_Vars,             ONLY: negativeJacobians,jacobianTolerance
 USE MOD_Output,                ONLY: InitOutput
-USE MOD_ReadInTools,           ONLY: IgnoredStrings
+USE MOD_ReadInTools,           ONLY: IgnoredStrings,FinalizeStrings
 USE MOD_Search,                ONLY: InitSearch
 #ifdef _OPENMP
 USE omp_lib
@@ -78,6 +78,9 @@ PROGRAM HOPR
 CALL IgnoredStrings()
 ! Now build mesh!
 CALL FillMesh()
+! Finalization
+CALL FinalizeStrings()
+! Output
 WRITE(UNIT_stdOut,'(132("="))')
 IF(negativeJacobians.GT.0) THEN
   WRITE(UNIT_stdOut,'(A,A,A)')' HOPR finished: Mesh "',TRIM(ProjectName)//'_mesh.h5','" written to HDF5 file.'

src/io_hdf5.f90

@@ -9,10 +9,12 @@
 ! /____//   /____//  /______________//  /____//           /____//   |_____/)    ,X`      XXX`
 ! )____)    )____)   )______________)   )____)            )____)    )_____)   ,xX`     .XX`
 !                                                                           xxX`      XXx
+! Copyright (C) 2023  Florian Hindenlang <hindenlang@gmail.com>
+! Copyright (C) 2023  Tobias Ott <tobias.ott@proton.me>
 ! Copyright (C) 2017 Claus-Dieter Munz <munz@iag.uni-stuttgart.de>
 ! 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
@@ -33,7 +35,7 @@ MODULE MOD_IO_HDF5
 IMPLICIT NONE
 PUBLIC
 !-----------------------------------------------------------------------------------------------------------------------------------
-! GLOBAL VARIABLES 
+! GLOBAL VARIABLES
 !-----------------------------------------------------------------------------------------------------------------------------------
 
 INTEGER, ALLOCATABLE           :: BCType(:,:)
@@ -50,7 +52,7 @@ MODULE MOD_IO_HDF5
 
 INTEGER,PARAMETER              :: ElemInfoSize=6        !number of entry in each line of ElemInfo
 INTEGER,PARAMETER              :: ELEM_Type=1           !entry position in ElemInfo
-INTEGER,PARAMETER              :: ELEM_Zone=2           
+INTEGER,PARAMETER              :: ELEM_Zone=2
 INTEGER,PARAMETER              :: ELEM_FirstSideInd=3
 INTEGER,PARAMETER              :: ELEM_LastSideInd=4
 INTEGER,PARAMETER              :: ELEM_FirstNodeInd=5
@@ -64,28 +66,54 @@ MODULE MOD_IO_HDF5
 INTEGER,PARAMETER              :: SIDE_BCID=5
 INTEGER,PARAMETER              :: nSidesElem(3:8)= & ! number of sides of an element nSides=nSidesElem(Elem%nNodes)
                                         (/3,4,5,5,0,6/) !tria,quad/tetra,pyra,prism,hex
-INTEGER,PARAMETER              :: LinMap(1:8,4:8)=RESHAPE(        & !CGNS -> IJK ordering for element corner nodes 
+INTEGER,PARAMETER              :: LinMap(1:8,4:8)=RESHAPE(        & !CGNS -> IJK ordering for element corner nodes
                                                                     !jNode=LinMap(iNode,Elem%nNodes)
                                        (/1,2,3,4,0,0,0,0,         & !Tet, one-to-one
-                                         1,2,4,3,5,0,0,0,         & !Pyra 
+                                         1,2,4,3,5,0,0,0,         & !Pyra
                                          1,2,3,4,5,6,0,0,         & !Prism one-to-one
                                          0,0,0,0,0,0,0,0,         & !nothing
                                          1,2,4,3,5,6,8,7/),(/8,5/)) !Hex
 
+INTEGER,PARAMETER              :: FEMElemInfoSize=4        !number of entry in each line of ElemInfo
+INTEGER,PARAMETER              :: FEMELEM_FirstEdgeInd=1
+INTEGER,PARAMETER              :: FEMELEM_LastEdgeInd=2
+INTEGER,PARAMETER              :: FEMELEM_FirstVertexInd=3
+INTEGER,PARAMETER              :: FEMELEM_LastVertexInd=4
+
+INTEGER,PARAMETER              :: EDGEInfoSize=3
+INTEGER,PARAMETER              :: EDGE_FEMEdgeID=1           
+INTEGER,PARAMETER              :: EDGE_offsetIndEdgeConnect=2 !entry position in EdgeConnectInfo
+INTEGER,PARAMETER              :: EDGE_lastIndEdgeConnect=3
+
+INTEGER,PARAMETER              :: EDGEConnectInfoSize=2
+INTEGER,PARAMETER              :: EDGEConnect_nbElemID=1            
+INTEGER,PARAMETER              :: EDGEConnect_nbLocEdgeID=2
+
+INTEGER,PARAMETER              :: VERTEXInfoSize=3
+INTEGER,PARAMETER              :: VERTEX_FEMVertexID=1            
+INTEGER,PARAMETER              :: VERTEX_offsetIndVertexConnect=2 !entry position in VertexConnectInfo
+INTEGER,PARAMETER              :: VERTEX_lastIndVertexConnect=3
+
+INTEGER,PARAMETER              :: VERTEXConnectInfoSize=2
+INTEGER,PARAMETER              :: VERTEXConnect_nbElemID=1            
+INTEGER,PARAMETER              :: VERTEXConnect_nbLocVertexID=2
+
 INTEGER,ALLOCATABLE            :: ElemInfo(:,:),SideInfo(:,:)
+INTEGER,ALLOCATABLE            :: FEMElemInfo(:,:), EdgeInfo(:,:),EdgeConnectInfo(:,:), VertexInfo(:,:),VertexConnectInfo(:,:)
 REAL,ALLOCATABLE               :: ElemWeight(:)
 REAL,ALLOCATABLE               :: ElemBarycenters(:,:)
 INTEGER,ALLOCATABLE            :: GlobalNodeIDs(:)
 REAL,ALLOCATABLE               :: NodeCoords(:,:)
 INTEGER,ALLOCATABLE            :: Elem_IJK(:,:)
 INTEGER                        :: nElems_IJK(3)
 INTEGER                        :: nGlobalElems
-INTEGER                        :: nElems,nSides,nNodes
+INTEGER                        :: nElems,nSides,nNodes,nEdges,nVertices
 INTEGER                        :: ElemCounter(2,11)
-INTEGER                        :: nSideIDs,nNodeIDs
+INTEGER                        :: nSideIDs,nNodeIDs,nEdgeIDs
+INTEGER                        :: nFEMEdgeIDs,nFEMEdgeConnections,nFEMVertexIDs,nFEMVertexConnections
 INTEGER                        :: nBCs
 LOGICAL                        :: curvedfound
-LOGICAL                        :: initMesh=.FALSE. 
+LOGICAL                        :: initMesh=.FALSE.
 
 INTERFACE INVMAP
   MODULE PROCEDURE INVMAP
@@ -106,7 +134,7 @@ MODULE MOD_IO_HDF5
 
 FUNCTION INVMAP(ID,nIDs,ArrID)
 !===================================================================================================================================
-! find the inverse Mapping p.e. NodeID-> entry in NodeMap (a sorted array of unique NodeIDs), using bisection 
+! find the inverse Mapping p.e. NodeID-> entry in NodeMap (a sorted array of unique NodeIDs), using bisection
 ! if Index is not in the range, -1 will be returned, if it is in the range, but is not found, 0 will be returned!!
 !===================================================================================================================================
 ! MODULES
@@ -132,7 +160,7 @@ FUNCTION INVMAP(ID,nIDs,ArrID)
   !WRITE(*,*)'WARNING, Node Index Not in local range -> set to -1'
   INVMAP=-1  ! not in the range!
   RETURN
-END IF 
+END IF
 IF(ID.EQ.ArrID(low))THEN
   INVMAP=low
 ELSEIF(ID.EQ.ArrID(up))THEN
@@ -151,7 +179,7 @@ FUNCTION INVMAP(ID,nIDs,ArrID)
     END IF
   END DO
 END IF
-END FUNCTION INVMAP 
+END FUNCTION INVMAP
 
 ! HFD5 STUFF
 SUBROUTINE OpenHDF5File(FileString,create)
@@ -189,7 +217,7 @@ END SUBROUTINE OpenHDF5File
 
 SUBROUTINE CloseHDF5File()
 !===================================================================================================================================
-! Close HDF5 file and groups 
+! Close HDF5 file and groups
 !===================================================================================================================================
 ! MODULES
 ! IMPLICIT VARIABLE HANDLING

src/mesh/cartmesh.f90

@@ -191,7 +191,7 @@ SUBROUTINE GetNewTetrahedron(CornerNode,cartMesh,l,m,n,ind)
   tetMap(:,3)=(/5,6,7,2/)
   tetMap(:,4)=(/7,8,5,4/)
   tetMap(:,5)=(/2,4,5,7/)
-  IF(ANY(cartmesh%nElems.GT.1)) STOP 'The selected mesh template is not periodic and can only be used for a single element.'
+  IF(ANY(cartmesh%nElems.GT.1)) CALL abort(__STAMP__,'The selected mesh template is not periodic and can only be used for a single element.')
 CASE(5)
   nTets=12
   !Center node of box is required
@@ -219,7 +219,7 @@ SUBROUTINE GetNewTetrahedron(CornerNode,cartMesh,l,m,n,ind)
 
 
 CASE DEFAULT
-  STOP 'The selected mesh template does not exist for tetrahedra.'
+  CALL abort(__STAMP__,'The selected mesh template does not exist for tetrahedra.')
 END SELECT
 
 DO i=1,nTets
@@ -540,8 +540,7 @@ SUBROUTINE CartesianMesh()
     ELSE !l0 active
       dx(i_Dim)=e1(i_Dim)/ABS(cartMesh%l0(i_Dim)) ! l/l0
       IF(dx(i_Dim) .LT. (1.-PP_RealTolerance)) THEN ! l0 > l
-        !CALL abort(__STAMP__,  & 
-         stop  'stretching error, length l0 longer than grid region, in direction '  !,i_Dim,999.)
+         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
         ! 
@@ -573,16 +572,14 @@ SUBROUTINE CartesianMesh()
             fac(i_Dim)= fac(i_Dim) - F/dF
             iter=iter+1
           END DO
-          IF(iter.GT.1000) STOP 'Newton iteration for computing the stretching function has failed.'
+          IF(iter.GT.1000) CALL abort(__STAMP__,'Newton iteration for computing the stretching function has failed.')
           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)  
     END IF
 
-    IF( ABS((nElems(i_Dim)-1.)*LOG(fac(i_Dim))/LOG(10.)) .GE. 4. )                      &
-      !CALL abort(__STAMP__, &
-       stop  'stretching error, length ratio > 1.0E4 in direction '  !,i_Dim,999.)
+    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')
     IF(ABS(fac(i_Dim)-1.) .GT. PP_RealTolerance ) THEN
       dx(i_Dim)=(1.-fac(i_Dim))/(1.-fac(i_Dim)**nElems(i_Dim)) ! first length to start
     ELSE !equidistant case
@@ -677,8 +674,7 @@ SUBROUTINE CartesianMesh()
             CALL GetNewHexahedron(CornerNode)
           END IF
         CASE DEFAULT
-          CALL abort(__STAMP__,&
-            'The specified element type is not known. Valid types: 104,105,106,108',CartMesh%ElemType)
+          CALL abort(__STAMP__,'The specified element type is not known. Valid types: 104,105,106,108',CartMesh%ElemType)
         END SELECT
       END DO !n
     END DO !m