getAreas.f90

Go to the documentation of this file.

subroutine getareas(areas, pts, npts, sps_in, axis)
    use constants
    use blockpointers
    use flowvarrefstate
    use inputtimespectral
    use communication
    use inputphysics
    use utils, only: setpointers
    implicit none
    !
    !      Local variables.
    !
    integer(kind=intType), intent(in) :: npts, sps_in
    real(kind=realtype), intent(in) :: pts(3, npts), axis(3)
    real(kind=realtype), intent(out) :: areas(3, npts)
 
    integer(kind=intType) :: mm, nn, i, j, ii, sps
    integer(kind=intType) :: iBeg, iEnd, jBeg, jEnd
    integer(kind=intType) :: lower_left, lower_right, upper_left, upper_right
    real(kind=realtype) :: da, fact, fact2
 
    areas = zero
    sps = sps_in
 
    ! Compute the local forces (or tractions). Take the scaling
    ! factor into account to obtain the forces in SI-units,
    ! i.e. Newton.
    ii = 0
    domains: do nn = 1, ndom
        call setpointers(nn, 1_inttype, sps)
        if (flowdoms(nn, 1_inttype, sps)%rightHanded) then
            fact2 = one
        else
            fact2 = -one
        end if
 
        ! Loop over the number of boundary subfaces of this block.
        bocos: do mm = 1, nbocos
 
            if (bctype(mm) == eulerwall .or. bctype(mm) == nswalladiabatic .or. &
                bctype(mm) == nswallisothermal) then
 
                select case (bcfaceid(mm))
 
                    ! NOTE: The 'fact' here are NOT the same as you will
                    ! find in ForcesAndMoment.f90. The reason is that, we
                    ! are not using points to si, sj, sk. Those have teh
                    ! normals pointing in the direction of increasing
                    ! {i,j,k}. Here we are evaluating the normal from
                    ! directly from the coordinates on the faces. As it
                    ! happens, the normals for the jMin and jMax faces are
                    ! flipped.
                case (imin)
                    fact = one
                case (imax)
                    fact = -one
                case (jmin)
                    fact = -one
                case (jmax)
                    fact = one
                case (kmin)
                    fact = one
                case (kmax)
                    fact = -one
                end select
 
                ! Store the cell range of the subfaces a bit easier.
                ! As only owned faces must be considered the nodal range
                ! in BCData must be used to obtain this data.
 
                jbeg = bcdata(mm)%jnBeg + 1; jend = bcdata(mm)%jnEnd
                ibeg = bcdata(mm)%inBeg + 1; iend = bcdata(mm)%inEnd
 
                ! Compute the dual area at each node. Just store in first dof
                do j = jbeg, jend ! This is a face loop
                    do i = ibeg, iend ! This is a face loop
 
                        ! Compute Normal
 
                        lower_left = ii + (j - jbeg) * (iend - ibeg + 2) + i - ibeg + 1
                        lower_right = lower_left + 1
                        upper_left = lower_right + iend - ibeg + 1
                        upper_right = upper_left + 1
 
                        call quad_area( &
                            pts(:, lower_left), pts(:, lower_right), &
                            pts(:, upper_left), pts(:, upper_right), &
                            axis, da)
                        da = fourth * da * fact * fact2
 
                        if (da > zero) then
                            ! Scatter to nodes
                            areas(1, lower_left) = areas(1, lower_left) + da
                            areas(1, lower_right) = areas(1, lower_right) + da
                            areas(1, upper_left) = areas(1, upper_left) + da
                            areas(1, upper_right) = areas(1, upper_right) + da
                        end if
                    end do
                end do
 
                ! Note how iBeg,iBeg is defined above... it is one MORE
                ! then the starting node (used for looping over faces, not
                ! nodes)
                ii = ii + (jend - jbeg + 2) * (iend - ibeg + 2)
 
            end if
        end do bocos
    end do domains
end subroutine getareas
 
subroutine getareasensitivity(darea, pts, npts, sps_in, axis)
    use constants
    use blockpointers
    use flowvarrefstate
    use inputtimespectral
    use communication
    use inputphysics
    use utils, only: setpointers
    implicit none
    !
    !      Local variables.
    !
    integer(kind=intType), intent(in) :: npts, sps_in
    real(kind=realtype), intent(in) :: pts(3, npts), axis(3)
    real(kind=realtype), intent(out) :: darea(3, npts)
 
    integer(kind=intType) :: mm, nn, i, j, ii, sps
    integer(kind=intType) :: iBeg, iEnd, jBeg, jEnd
    integer(kind=intType) :: lower_left, lower_right, upper_left, upper_right
    real(kind=realtype) :: area, areab, pt1b(3), pt2b(3), pt3b(3), pt4b(3)
    real(kind=realtype) :: fact, fact2, da
 
    darea = zero
    sps = sps_in
 
    ! Compute the local forces (or tractions). Take the scaling
    ! factor into account to obtain the forces in SI-units,
    ! i.e. Newton.
    ii = 0
    domains: do nn = 1, ndom
        call setpointers(nn, 1_inttype, sps)
        if (flowdoms(nn, 1_inttype, sps)%rightHanded) then
            fact2 = one
        else
            fact2 = -one
        end if
 
        ! Loop over the number of boundary subfaces of this block.
        bocos: do mm = 1, nbocos
 
            if (bctype(mm) == eulerwall .or. bctype(mm) == nswalladiabatic .or. &
                bctype(mm) == nswallisothermal) then
 
                select case (bcfaceid(mm))
 
                    ! NOTE: The 'fact' here are NOT the same as you will
                    ! find in ForcesAndMoment.f90. The reason is that, we
                    ! are not using points to si, sj, sk. Those have teh
                    ! normals pointing in the direction of increasing
                    ! {i,j,k}. Here we are evaluating the normal from
                    ! directly from the coordinates on the faces. As it
                    ! happens, the normals for the jMin and jMax faces are
                    ! flipped.
                case (imin)
                    fact = one
                case (imax)
                    fact = -one
                case (jmin)
                    fact = -one
                case (jmax)
                    fact = one
                case (kmin)
                    fact = one
                case (kmax)
                    fact = -one
                end select
 
                ! Store the cell range of the subfaces a bit easier.
                ! As only owned faces must be considered the nodal range
                ! in BCData must be used to obtain this data.
 
                jbeg = bcdata(mm)%jnBeg + 1; jend = bcdata(mm)%jnEnd
                ibeg = bcdata(mm)%inBeg + 1; iend = bcdata(mm)%inEnd
 
                do j = jbeg, jend ! This is a face loop
                    do i = ibeg, iend ! This is a face loop
 
                        ! Extract 4 corner points
                        lower_left = ii + (j - jbeg) * (iend - ibeg + 2) + i - ibeg + 1
                        lower_right = lower_left + 1
                        upper_left = lower_right + iend - ibeg + 1
                        upper_right = upper_left + 1
 
                        ! Compute actual area since we need to know to
                        ! include or not. The reverse mode calc does NOT
                        ! compute area
                        call quad_area( &
                            pts(:, lower_left), pts(:, lower_right), &
                            pts(:, upper_left), pts(:, upper_right), &
                            axis, da)
 
                        da = fourth * da * fact * fact2
                        if (da > zero) then
 
                            areab = one
                            call quad_area_b( &
                                pts(:, lower_left), pt1b, &
                                pts(:, lower_right), pt2b, &
                                pts(:, upper_left), pt3b, &
                                pts(:, upper_right), pt4b, &
                                axis, area, areab)
 
                            darea(:, lower_left) = darea(:, lower_left) + pt1b
                            darea(:, lower_right) = darea(:, lower_right) + pt2b
                            darea(:, upper_left) = darea(:, upper_left) + pt3b
                            darea(:, upper_right) = darea(:, upper_right) + pt4b
                        end if
                    end do
                end do
 
                ! Note how iBeg,iBeg is defined above... it is one MORE
                ! then the starting node (used for looping over faces, not
                ! nodes)
                ii = ii + (jend - jbeg + 2) * (iend - ibeg + 2)
 
            end if
        end do bocos
    end do domains
end subroutine getareasensitivity
 
subroutine quad_area(p1, p2, p3, p4, axis, area)
    ! Kernel-level function to get area of quad defined by 4 points
    ! projected onto plane defined by axis. Only +ve areas are computed.
    use constants
    implicit none
 
    ! I/O
    real(kind=realtype), intent(in) :: p1(3), p2(3), p3(3), p4(3), axis(3)
    real(kind=realtype), intent(out) :: area
 
    ! Working
    real(kind=realtype) :: v1(3), v2(3), sss(3)
    ! Vectors for Cross Product
 
    v1(:) = p4 - p1
    v2(:) = p3 - p2
 
    ! Cross Product
    sss(1) = half * (v1(2) * v2(3) - v1(3) * v2(2))
    sss(2) = half * (v1(3) * v2(1) - v1(1) * v2(3))
    sss(3) = half * (v1(1) * v2(2) - v1(2) * v2(1))
 
    area = sss(1) * axis(1) + sss(2) * axis(2) + sss(3) * axis(3)
 
end subroutine quad_area
 
!        Generated by TAPENADE     (INRIA, Tropics team)
!  Tapenade 3.4 (r3375) - 10 Feb 2010 15:08
!
!  Differentiation of quad_area in reverse (adjoint) mode:
!   gradient     of useful results: area
!   with respect to varying inputs: area p1 p2 p3 p4
!   RW status of diff variables: area:in-zero p1:out p2:out p3:out
!                p4:out
SUBROUTINE quad_area_b(p1, p1b, p2, p2b, p3, p3b, p4, p4b, axis, area, &
&  areab)
    use constants
    IMPLICIT NONE
! Kernel-level function to get area of quad defined by 4 points
! projected onto plane defined by axis. Only +ve areas are computed.
! I/O
    REAL(kind=realtype), INTENT(IN) :: p1(3), p2(3), p3(3), p4(3), axis(3)
    REAL(kind=realtype) :: p1b(3), p2b(3), p3b(3), p4b(3)
    REAL(kind=realtype) :: area
    REAL(kind=realtype) :: areab
! Working
    REAL(kind=realtype) :: v1(3), v2(3), sss(3)
    REAL(kind=realtype) :: v1b(3), v2b(3), sssb(3)
    REAL(kind=realtype) :: tempb1
    REAL(kind=realtype) :: tempb0
    INTRINSIC abs
    REAL(kind=realtype) :: tempb
! Vectors for Cross Product
    v1(:) = p4 - p1
    v2(:) = p3 - p2
! Cross Product
    sss(1) = half * (v1(2) * v2(3) - v1(3) * v2(2))
    sss(2) = half * (v1(3) * v2(1) - v1(1) * v2(3))
    sss(3) = half * (v1(1) * v2(2) - v1(2) * v2(1))
    IF (sss(1) * axis(1) + sss(2) * axis(2) + sss(3) * axis(3) .GE. 0.) THEN
        sssb = 0.0
        sssb(1) = axis(1) * areab
        sssb(2) = axis(2) * areab
        sssb(3) = axis(3) * areab
    ELSE
        sssb = 0.0
        sssb(1) = -(axis(1) * areab)
        sssb(2) = -(axis(2) * areab)
        sssb(3) = -(axis(3) * areab)
    END IF
    v1b = 0.0
    v2b = 0.0
    tempb = half * sssb(3)
    v1b(1) = v2(2) * tempb
    v2b(2) = v1(1) * tempb
    v1b(2) = -(v2(1) * tempb)
    sssb(3) = 0.0
    tempb0 = half * sssb(2)
    v2b(1) = v1(3) * tempb0 - v1(2) * tempb
    v1b(3) = v1b(3) + v2(1) * tempb0
    v1b(1) = v1b(1) - v2(3) * tempb0
    sssb(2) = 0.0
    tempb1 = half * sssb(1)
    v2b(3) = v2b(3) + v1(2) * tempb1 - v1(1) * tempb0
    v1b(2) = v1b(2) + v2(3) * tempb1
    v1b(3) = v1b(3) - v2(2) * tempb1
    v2b(2) = v2b(2) - v1(3) * tempb1
    p2b = 0.0
    p3b = 0.0
    p3b = v2b(:)
    p2b = -v2b(:)
    p1b = 0.0
    p4b = 0.0
    p4b = v1b(:)
    p1b = -v1b(:)
    areab = 0.0
END SUBROUTINE quad_area_b