Axisymmetric mesh warping capability

idwarp/AxisymmetricMesh_C.py

@@ -0,0 +1,22 @@
+import os
+
+from . import MExt
+from .AxisymmetricMesh import AxisymmetricMesh
+
+
+class AxisymmetricMesh_C(AxisymmetricMesh):
+    """
+    Represents a (Complex) Axisymmetric mesh
+    """
+
+    def __init__(self, *args, **kwargs):
+        """Initialize the object."""
+
+        debug = False
+        if "debug" in kwargs:
+            debug = kwargs["debug"]
+
+        curDir = os.path.basename(os.path.dirname(os.path.realpath(__file__)))
+        self.warp = MExt.MExt("libidwarp_cs", curDir, debug=debug)._module
+        AxisymmetricMesh.__init__(self, *args, **kwargs)
+        self.dtype = "D"

idwarp/__init__.py

@@ -1,8 +1,10 @@
 __version__ = "2.6.2"
 
+from .AxisymmetricMesh import AxisymmetricMesh
+from .AxisymmetricMesh_C import AxisymmetricMesh_C
 from .UnstructuredMesh import USMesh
 from .UnstructuredMesh_C import USMesh_C
 from .MultiUnstructuredMesh import MultiUSMesh
 from .MultiUnstructuredMesh_C import MultiUSMesh_C
 
-__all__ = ["USMesh", "MultiUSMesh", "USMesh_C", "MultiUSMesh_C"]
+__all__ = ["AxisymmetricMesh", "USMesh", "MultiUSMesh", "AxisymmetricMesh_C", "USMesh_C", "MultiUSMesh_C"]

src/IO/readStructuredCGNS.F90

@@ -28,8 +28,9 @@ subroutine readStructuredCGNS(cg)
 #endif
     integer(kind=intType), dimension(:), allocatable :: surfaceNodes, localSurfaceNodes
     integer(kind=intType), dimension(:, :), allocatable :: sizes
-    integer(kind=intType) :: nSurf, nodeCount, nConn, ni, nj, nx, ny, nz
+    integer(kind=intType) :: nSurf, nodeCount, nConn, ni, nj, nx, ny, nz, bocoIdx
     integer(kind=intType) :: status(MPI_STATUS_SIZE)
+    integer(kind=intType) :: nTotalBocos
     logical :: lowFace
 
     ! Only do reading on root proc:
@@ -66,6 +67,8 @@ subroutine readStructuredCGNS(cg)
 
         nSurf = 0
         nConn = 0
+        nTotalBocos = 0
+        bocoIdx = 0
         zoneLoop1: do iZone = 1, nZones
             call cg_zone_read_f(cg, base, iZone, zoneName, dims, ierr)
             if (ierr .eq. CG_ERROR) call cg_error_exit_f
@@ -77,6 +80,7 @@ subroutine readStructuredCGNS(cg)
             call cg_nbocos_f(cg, base, iZone, nbocos, ierr)
             if (ierr .eq. ERROR) call cg_error_exit_f
             allocate (zones(iZone)%bocos(nbocos))
+            nTotalBocos = nTotalBocos + nbocos
 
             bocoLoop1: do boco = 1, nbocos
 
@@ -137,6 +141,9 @@ subroutine readStructuredCGNS(cg)
         allocate (surfacePoints(3 * nSurf))
         allocate (surfaceConn(4 * nConn))
 
+        ! Allocate the flattened BC types storage
+        allocate (bocoTypes(nTotalBocos))
+
         ! Reset the counters here
         nSurf = 0
         offSet = 0
@@ -191,7 +198,9 @@ subroutine readStructuredCGNS(cg)
             end do
 
             bocoLoop2: do boco = 1, size(zones(iZone)%bocos)
+                bocoIdx = bocoIdx + 1
                 bocoType = zones(iZone)%bocos(boco)%type
+                bocoTypes(bocoIdx) = bocoType
                 pts = zones(iZone)%bocos(boco)%ptRange
 
                 ! Flag the surface nodes
@@ -284,14 +293,81 @@ subroutine readStructuredCGNS(cg)
         if (ierr .eq. CG_ERROR) call cg_error_exit_f
         deallocate (sizes)
     else
-        ! Allocate these to zero so we can just blindly dealloc later
+        !Allocate these to zero so we can just blindly dealloc later
         allocate (surfacePoints(0), surfaceConn(0))
     end if rootProc
 
     ! Communicate the total number of nodes to everyone
     call MPI_bcast(nNodes, 1, MPI_INTEGER, 0, warp_comm_world, ierr)
     call EChk(ierr, __FILE__, __LINE__)
 
+    ! For axisymm cases we need to communicate some of the zone data to every proc
+    axisymmetric: if (axisymm) then
+        ! Total number of zones
+        call MPI_bcast(nZones, 1, MPI_INTEGER, 0, warp_comm_world, ierr)
+        call EChk(ierr, __FILE__, __LINE__)
+
+        ! Allocate the zones
+        if (myid /= 0) then
+            allocate (zones(nZones))
+        end if
+
+        ! Communicate the zone data
+        do i = 1, nZones
+            ! Communicate the zone data
+            call MPI_bcast(zones(i)%il, 1, MPI_INTEGER, 0, warp_comm_world, ierr)
+            call EChk(ierr, __FILE__, __LINE__)
+            call MPI_bcast(zones(i)%jl, 1, MPI_INTEGER, 0, warp_comm_world, ierr)
+            call EChk(ierr, __FILE__, __LINE__)
+            call MPI_bcast(zones(i)%kl, 1, MPI_INTEGER, 0, warp_comm_world, ierr)
+            call EChk(ierr, __FILE__, __LINE__)
+            call MPI_bcast(zones(i)%name, maxCGNSNameLen, MPI_CHARACTER, 0, warp_comm_world, ierr)
+            call EChk(ierr, __FILE__, __LINE__)
+            call MPI_bcast(zones(i)%nVertices, 1, MPI_INTEGER, 0, warp_comm_world, ierr)
+            call EChk(ierr, __FILE__, __LINE__)
+            call MPI_bcast(zones(i)%nElements, 1, MPI_INTEGER, 0, warp_comm_world, ierr)
+            call EChk(ierr, __FILE__, __LINE__)
+
+            ! Get the number of bocos for this zone from the root proc
+            if (myID == 0) nbocos = size(zones(i)%bocos)
+
+            ! Communicate the number of bocos
+            call MPI_bcast(nbocos, 1, MPI_INTEGER, 0, warp_comm_world, ierr)
+            call EChk(ierr, __FILE__, __LINE__)
+
+            if (myid /= 0) then
+                allocate (zones(i)%bocos(nbocos))
+
+                ! Nullify section pointers since we won't have any for structured mesh
+                nullify (zones(i)%sections)
+            end if
+
+            do j = 1, nbocos
+                ! Communicate all boco data
+                call MPI_bcast(zones(i)%bocos(j)%name, maxCGNSNameLen, MPI_CHARACTER, 0, warp_comm_world, ierr)
+                call EChk(ierr, __FILE__, __LINE__)
+                call MPI_bcast(zones(i)%bocos(j)%type, 1, MPI_INTEGER, 0, warp_comm_world, ierr)
+                call EChk(ierr, __FILE__, __LINE__)
+                call MPI_bcast(zones(i)%bocos(j)%ptRange, 6, MPI_INTEGER, 0, warp_comm_world, ierr)
+                call EChk(ierr, __FILE__, __LINE__)
+                call MPI_bcast(zones(i)%bocos(j)%family, maxCGNSNameLen, MPI_CHARACTER, 0, warp_comm_world, ierr)
+                call EChk(ierr, __FILE__, __LINE__)
+                call MPI_bcast(zones(i)%bocos(j)%nBCElem, 1, MPI_INTEGER, 0, warp_comm_world, ierr)
+                call EChk(ierr, __FILE__, __LINE__)
+                call MPI_bcast(zones(i)%bocos(j)%nBCNodes, 1, MPI_INTEGER, 0, warp_comm_world, ierr)
+                call EChk(ierr, __FILE__, __LINE__)
+
+                ! Allocate empty data so we don't get errors trying to deallocate later
+                if (myid /= 0) then
+                    allocate (zones(i)%bocos(j)%BCElements(0))
+                    nullify (zones(i)%bocos(j)%elemPtr)
+                    nullify (zones(i)%bocos(j)%elemConn)
+                    allocate (zones(i)%bocos(j)%elemNodes(0, 0))
+                end if
+            end do
+        end do
+    end if axisymmetric
+
     rootProc2: if (myid == 0) then
 
         istart = 0 + 1

src/f2py/libidwarp.pyf

@@ -165,12 +165,20 @@ python module libidwarp ! in
          use constants
          integer(kind=inttype) :: warpmeshdof
          integer(kind=inttype) :: solvermeshdof
-         integer(Kind=inttype) :: commonmeshdof
+         integer(kind=inttype) :: commonmeshdof
+         logical, optional :: axisymm=.false.
+         real(kind=realtype) :: axisymmangle
+         real(kind=realtype) dimension(3) :: axisymmaxis
+
+         subroutine setmirrorfamily(famname) ! in :idwarp:../modules/gridData.f90
+           character*(*) intent(in) :: famname
+         end subroutine setmirrorfamily
        end module griddata
 
        module cgnsgrid
          real(kind=realtype), dimension(:), allocatable :: surfacepoints
          integer(kind=inttype), dimension(:), allocatable :: surfaceconn
+         integer(kind=inttype), dimension(:), allocatable :: bocoTypes
          logical, dimension(:), allocatable :: surfaceiswall
          logical, dimension(:), allocatable :: surfaceissymm
          integer(kind=inttype), dimension(:,:), allocatable :: surfacesizes
@@ -204,7 +212,6 @@ python module libidwarp ! in
        module kd_tree
          integer(kind=inttype) :: bucket_size
        end module kd_tree
-
     end interface 
 #ifdef USE_COMPLEX
 end python module libidwarp_cs

src/modules/cgnsGrid.F90

@@ -45,6 +45,9 @@ module cgnsGrid
     ! List of the zones
     type(zoneDataType), dimension(:), allocatable :: zones
 
+    ! All BC types (Need this for axisymmetric meshes)
+    integer(kind=intType), dimension(:), allocatable :: bocoTypes
+
     ! Deduced information of the surfaces
     real(kind=realType), dimension(:), allocatable :: surfacePoints
     logical, dimension(:), allocatable :: surfaceIsWall

src/modules/gridData.F90

@@ -39,7 +39,7 @@ module gridData
 
 #ifndef USE_TAPENADE
     real(kind=realType), pointer, dimension(:) :: XsPtrb, XsPtrd
-    real(kind=realType), pointer, dimension(:) :: XvPtrb, XVPtrd
+    real(kind=realType), pointer, dimension(:) :: XvPtrb, XvPtrd
 #endif
 
     ! Sizes of the three different mesh sizes:
@@ -60,6 +60,18 @@ module gridData
     integer(kind=intType) :: nLoop
     real(kind=realType), dimension(:, :), allocatable :: symmPts, symmNormals
 
+    ! Axisymmetric Information
+    logical :: axiSymm = .false.
+    real(kind=realType) :: axiSymmAngle
+    real(kind=realType) :: axiSymmAxis(3)
+    character(len=maxStringLen) :: mirrorFamily
+
+contains
+
+    subroutine setMirrorFamily(famName)
+        character(len=*), intent(in) :: famName
+        mirrorFamily = famName
+    end subroutine setMirrorFamily
 end module gridData
 
 module plot3dSurface

src/modules/kd_tree.F90

@@ -1637,7 +1637,7 @@ subroutine computeNodalProperties(tp, initialPoint)
             else
                 ! Now get the rotation Matrix
                 if (useRotations .and. .not. tp%isCorner(i)) then
-                    call getRotationMatrix3d(tp%normals0(:, i), tp%normals(:, i), tp%Mi(:, :, i))
+                    call getRotationMatrixVecVec(tp%normals0(:, i), tp%normals(:, i), tp%Mi(:, :, i))
                 else
                     ! Set identity
                     tp%Mi(:, :, i) = zero

src/utils/releaseMemory.F90

@@ -106,6 +106,11 @@ subroutine releaseMemory
         ! Delete zone array itself
         deallocate (zones)
 
+        ! Delete the boco types array
+        if (allocated(bocoTypes)) then
+            deallocate (bocoTypes)
+        end if
+
         ! Delete the surface-stuff in the cgns grid if necessary
         if (allocated(surfacePoints)) then
             deallocate (surfacePoints)

src/utils/vectorUtils.f90

@@ -28,7 +28,7 @@ subroutine getMag(V, mag)
 
 end subroutine getMag
 
-subroutine getRotationMatrix3d(v1, v2, Mi)
+subroutine getRotationMatrixVecVec(v1, v2, Mi)
     use precision
     use constants
     implicit none
@@ -100,4 +100,54 @@ subroutine getRotationMatrix3d(v1, v2, Mi)
 
     Mi = Mi + sin(angle) * A + (one - cos(angle)) * C
 
-end subroutine getRotationMatrix3d
+end subroutine getRotationMatrixVecVec
+
+subroutine getRotationMatrixAngleAxis(angle, axis, Mi)
+    use precision
+    use constants
+    implicit none
+
+    ! Input variables
+    real(kind=realType), intent(in) :: axis(3), angle
+
+    ! Output variables
+    real(kind=realType), intent(out) :: Mi(3, 3)
+
+    ! Working variables
+    real(kind=realType), dimension(3, 3) :: A, C
+    real(kind=realType) :: overNorm, axisNorm(3), axisMag
+
+    call getMag(axis, axisMag)
+    overNorm = 1 / axisMag
+    axisNorm = axis * overNorm
+
+    ! Compute the rotation matrix
+    A(1, 1) = zero
+    A(1, 2) = -axisNorm(3)
+    A(1, 3) = axisNorm(2)
+    A(2, 1) = axisNorm(3)
+    A(2, 2) = zero
+    A(2, 3) = -axisNorm(1)
+    A(3, 1) = -axisNorm(2)
+    A(3, 2) = axisNorm(1)
+    A(3, 3) = zero
+
+    !C = A*A
+    C(1, 1) = A(1, 1) * A(1, 1) + A(1, 2) * A(2, 1) + A(1, 3) * A(3, 1)
+    C(1, 2) = A(1, 1) * A(1, 2) + A(1, 2) * A(2, 2) + A(1, 3) * A(3, 2)
+    C(1, 3) = A(1, 1) * A(1, 3) + A(1, 2) * A(2, 3) + A(1, 3) * A(3, 3)
+    C(2, 1) = A(2, 1) * A(1, 1) + A(2, 2) * A(2, 1) + A(2, 3) * A(3, 1)
+    C(2, 2) = A(2, 1) * A(1, 2) + A(2, 2) * A(2, 2) + A(2, 3) * A(3, 2)
+    C(2, 3) = A(2, 1) * A(1, 3) + A(2, 2) * A(2, 3) + A(2, 3) * A(3, 3)
+    C(3, 1) = A(3, 1) * A(1, 1) + A(3, 2) * A(2, 1) + A(3, 3) * A(3, 1)
+    C(3, 2) = A(3, 1) * A(1, 2) + A(3, 2) * A(2, 2) + A(3, 3) * A(3, 2)
+    C(3, 3) = A(3, 1) * A(1, 3) + A(3, 2) * A(2, 3) + A(3, 3) * A(3, 3)
+
+    ! Rodrigues formula for the rotation matrix
+    Mi = zero
+    Mi(1, 1) = one
+    Mi(2, 2) = one
+    Mi(3, 3) = one
+
+    Mi = Mi + sin(angle) * A + (one - cos(angle)) * C
+end subroutine getRotationMatrixAngleAxis

src/warp/Makefile

@@ -11,6 +11,7 @@ include ${RULES_FILE}
 vpath %.o $(OBJDIR)
 
 FILES_TO_COMPLEXIFY = \
+	axisymm.F90\
 	getdXs.F90\
 	getElementProps.F90\
 	getSurfaceCoordinates.F90\