RBF curving (#6)

* Added Volume curving using a radial basis function interpolation of the displacement of the curved boundary sides to the volume points.

* Add remaining RBF types, scale global support RBF types with SupportRadius (set to 1 to recover original formulation)

* Fixed wrong indizes in RBF coefficient calculations. Some RBFs need special threatment to avoid NANs since calculations of log(0) can occur.

* Also use the corner nodes of all elements as base points for the RBF interpolation. This way, e.g. the cell height in a BL mesh will not be affected.

* Move RBF curving after zCorrection

* RBF curving is only done inside of a user-defined box

* Update and declutterng of comments.

* Also only evaluate the interpolation inside the box

* Use LAPACK routines to directly solve the linear System, not to invert the Matrix and apply using MATMUL.

* Z correction for RBF curving to undo any 3D effects introduced for extruded meshes

* Fix bug in RBF curving, different conditions for the bounding box where applied when counting the RBF control points and when storing them, leading to a singular matrix

* Allow for multiple RBF boxes

* Only curve all or no nodes of an element

* Example for RBF curving

* Correctly count RBF control if the first is outside the bounding box but others inside

* Only add complete elements, remove print statement

* Update .gitlab-ci.yml

* Clean standard values for RBF bounding box limits xlim and ylim and cleaner fomatting of parameter file for RBF tutorial

* Move RBF tutorial to avoid clash with numbering of tutorials added prior.

* Address comments of code review. Remove distance function between points by the NORM2 intrinsic, clean up the if-conditions in the RBF implementation and replace the bary-center with the bounding box of an elemenet to check whether it resides within the current RBF domain

* Reduced image file sizes (changed total file size from 25M to 3.4M).

---------

Co-authored-by: iagkrais <krais@iag.uni-stuttgart.de>
Co-authored-by: Patrick Kopper <patrick.kopper@ila.uni-stuttgart.de>
Co-authored-by: Marcel Blind <marcel.blind@iag.uni-stuttgart.de>
Co-authored-by: Stephen Copplestone <stephen@copplestone.de>

src/globals.f90

@@ -94,6 +94,10 @@ MODULE MOD_Globals
    MODULE PROCEDURE getInverse
 END INTERFACE
 
+INTERFACE SolveLinSys
+   MODULE PROCEDURE SolveLinSys
+END INTERFACE
+
 INTERFACE
   SUBROUTINE setstacksizeunlimited() BIND(C)
   END SUBROUTINE setstacksizeunlimited
@@ -384,4 +388,42 @@ FUNCTION GetInverse(dim1,A) RESULT(Ainv)
   END IF
 END FUNCTION GetInverse
 
+
+SUBROUTINE SolveLinSys(dim1,nRHS,A,B)
+!===================================================================================================================================
+! Solve liner system A*X=B using LAPACK routines. A will be overwritten in the process! Multiple RHS possible.
+!===================================================================================================================================
+! MODULES
+! IMPLICIT VARIABLE HANDLING
+IMPLICIT NONE
+!-----------------------------------------------------------------------------------------------------------------------------------
+! INPUT VARIABLES
+INTEGER, INTENT(IN) :: dim1   ! size of matrix A
+INTEGER,INTENT(IN)  :: nRHS   ! number of right hand sides
+REAL,INTENT(INOUT)  :: A(dim1,dim1) ! System matrix
+REAL,INTENT(INOUT)  :: B(dim1,nRHS) ! Right hand side of the system
+!-----------------------------------------------------------------------------------------------------------------------------------
+!-----------------------------------------------------------------------------------------------------------------------------------
+! LOCAL VARIABLES
+INTEGER            :: IPIV(dim1),INFO,lwork ! ?
+REAL               :: WORK(dim1*dim1)  ! ?
+!===================================================================================================================================
+! 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.
+CALL DGETRF(dim1, dim1, A, dim1, IPIV, INFO)
+
+IF (INFO.NE.0) THEN
+  CALL abort(__STAMP__,'SolveLinSys(dim1,nRHS,A,B): MATRIX A IS NUMERICALLY SINGULAR! INFO = ',IntInfoOpt=INFO)
+END IF
+
+! DGETRES computes the solution to the liner system A*X=B using the factorization
+! computed by DGETRF.
+CALL DGETRS('N', dim1, nRHS, A, dim1, IPIV, B, dim1, INFO)
+
+IF (INFO.NE.0) THEN
+  CALL abort(__STAMP__,'SolveLinSys(dim1,nRHS,A,B): SOLVING LINEAR SYSTEM FAILED! INFO = ',IntInfoOpt=INFO)
+END IF
+END SUBROUTINE SolveLinSys
+
+
 END MODULE MOD_Globals

src/mesh/mesh.f90

@@ -247,6 +247,28 @@ SUBROUTINE InitMesh()
       CALL abort(__STAMP__,&
         'Specified curving method does not exist. =1: NormalVectors, =3: SplitElemFile, =4: SpecElemFile')
     END SELECT
+
+    ! Volume curving by radial basis functions
+    useRBF = GETLOGICAL('useRBF','.FALSE.')
+    IF (useRBF) THEN
+      nRBFBoxes = CNTSTR('RBFType','0')
+
+      ALLOCATE(SupportRadius(nRBFBoxes))
+      ALLOCATE(RBFType(      nRBFBoxes))
+      ALLOCATE(xlim(       2,nRBFBoxes))
+      ALLOCATE(ylim(       2,nRBFBoxes))
+
+      ! If xlim or ylim not given, apply RBF curving to whole mesh
+      WRITE(DefStr,'(E21.11,A1,E21.11)')-HUGE(1.),',',HUGE(1.)
+
+      DO i=1,nRBFBoxes
+        SupportRadius(i) = GETREAL('SupportRadius')
+        RBFType(i) = GETINT('RBFType')
+        xlim(:, i) = GETREALARRAY('xlim',2,TRIM(DefStr))
+        ylim(:, i) = GETREALARRAY('ylim',2,TRIM(DefStr))
+      END DO
+    END IF
+
   END IF
   doExactSurfProjection=GETLOGICAL('doExactSurfProjection','.FALSE.')
   IF(doExactSurfProjection)THEN
@@ -448,6 +470,7 @@ SUBROUTINE fillMesh()
 USE MOD_Mesh_PostDeform,  ONLY: PostDeform
 USE MOD_Output_HDF5,      ONLY: WriteMeshToHDF5
 USE MOD_Mesh_Jacobians,   ONLY: CheckJacobians
+USE MOD_RBF,              ONLY: RBFVolumeCurving
 USE MOD_Readin_ANSA
 USE MOD_Readin_CGNS
 USE MOD_Readin_Gambit
@@ -470,7 +493,7 @@ SUBROUTINE fillMesh()
 TYPE(tElem),POINTER :: Elem  ! ?
 TYPE(tSide),POINTER :: Side    ! ?
 LOGICAL             :: curvedFound  ! ?
-INTEGER             :: i,iElem  ! ?
+INTEGER             :: i,iElem,iRBFBox  ! ?
 !===================================================================================================================================
 CALL Timer(.TRUE.)
 WRITE(UNIT_stdOut,'(132("="))')
@@ -668,7 +691,22 @@ SUBROUTINE fillMesh()
 ! get element types
 CALL FindElemTypes()
 ! correct displacement in z (e.g. for periodic/2.5D)
-IF(doZcorrection) CALL zCorrection()
+IF(doZcorrection) CALL zCorrection(InitZOrient_In=.TRUE.)
+
+! Call RBF curving after GlobalUniqueNodes, since no double entries in the RBF nodes are allowed
+IF (useCurveds.AND.useRBF) THEN
+  DO iRBFBox=1,nRBFBoxes
+    CALL RBFVolumeCurving(iRBFBox)
+  END DO
+  ! Deallocate all RBF quantities
+  DEALLOCATE(RBFType)
+  DEALLOCATE(SupportRadius)
+  DEALLOCATE(xlim)
+  DEALLOCATE(ylim)
+  ! Call z-Correction again to remove 3D effects of interpolation, zPeriodic is already done
+  zPeriodic = .FALSE.
+  IF(doZcorrection) CALL zCorrection(InitZOrient_In=.FALSE.)
+END IF
 
 CALL CheckNodeConnectivity()
 

src/mesh/mesh_vars.f90

@@ -255,6 +255,12 @@ MODULE MOD_Mesh_Vars
 LOGICAL                        :: rebuildCurveds         ! switch .TRUE.= if curveds are already present in the mesh, delete them
                                                          ! and rebuild them using our methods
 LOGICAL                        :: meshIsAlreadyCurved    ! flag: mesh is already curved (GMSH, HDF5, block CGNS)
+LOGICAL                        :: useRBF                 ! Volume curving using interpolation of surface curving by
+                                                         ! radial basis functions
+INTEGER                        :: nRBFBoxes              ! Number of RBF Bounding Boxes
+INTEGER,ALLOCATABLE            :: RBFType(:)             ! Type of radial basis function of each RBFBoundingBox
+REAL,ALLOCATABLE               :: SupportRadius(:)       ! Support radius of radial basis functions of each RBFBoundingBox
+REAL,ALLOCATABLE               :: xlim(:,:),ylim(:,:)    ! Only inside this box the RBF curving is used (2,nRBFBoxes)
 LOGICAL                        :: InnerElemStretch       ! for cartmeshes, apply stretching also to inner element nodes
 !-----------------------------------------------------------------------------------------------------------------------------------
 ! CURVE GRID GENERATOR
@@ -297,6 +303,7 @@ MODULE MOD_Mesh_Vars
 LOGICAL                        :: doZcorrection 
 REAL                           :: zstart 
 LOGICAL                        :: zPeriodic 
+INTEGER,ALLOCATABLE            :: whichdirArr(:),orientArr(:)
 !-----------------------------------------------------------------------------------------------------------------------------------
 ! Splitting of Elements 
 !-----------------------------------------------------------------------------------------------------------------------------------