๐Ÿ’ป Fortran Source Code Library

We currently offer 172 open-source, production-grade Fortran codes for offline testing. Run calculations locally on your own machine, view code structure, read technical explanations, and download compilation packages including sample input files.

Mesh Cell Aspect Ratio

Core Numerical Engine in Fortran 90 โ€ข 30 total downloads

aspect_ratio.f90
! =========================================================================
! Source File: aspect_ratio.f90
! =========================================================================

!==============================================================================
! ThermoFluidCalc โ€” Calculator #27 : Mesh Aspect Ratio
!==============================================================================
! Physics : The aspect ratio (AR) measures how stretched / elongated a
!           mesh cell is.  Ideal AR = 1 (equilateral tri or square quad).
!
!   Triangle AR (inradius-based):
!     AR = max_edge / (2 * sqrt(3) * r)
!     where r = 2*Area / perimeter   (inradius)
!     For an equilateral triangle AR = 1.
!
!   Quad AR (opposite-edge-pair):
!     avg1 = (|e1| + |e3|) / 2    (pair of opposite edges)
!     avg2 = (|e2| + |e4|) / 2
!     AR   = max(avg1, avg2) / min(avg1, avg2)
!
!   Simple edge-ratio (both):
!     AR_edge = longest_edge / shortest_edge
!
! Quality rating:
!   1.0 - 1.5   Excellent
!   1.5 - 3.0   Good
!   3.0 - 5.0   Acceptable
!   5.0 - 10.0  Poor
!   > 10.0      Bad
!
! Reference : Gupta, ยง5.1.1
!
! Modes:
!   1 = Single triangle  (3 vertices 2-D)
!   2 = Single quad       (4 vertices 2-D)
!   3 = Batch             (N mixed cells)
!
! Build:
!   gfortran -O2 -o aspect_ratio aspect_ratio.f90
!==============================================================================
program aspect_ratio
  implicit none

  integer, parameter :: dp = selected_real_kind(15, 307)
  real(dp), parameter :: SQRT3 = 1.7320508075688772_dp
  integer, parameter :: MAX_CELLS = 10000

  integer  :: mode, N, i, nv
  real(dp) :: vx(4), vy(4)
  real(dp) :: edges(4), perimeter, area
  real(dp) :: emax, emin, AR_edge, AR_main
  real(dp) :: inrad, circumR
  real(dp) :: avg1, avg2, d1, d2, diag_ratio
  character(len=20) :: rating

  ! Batch arrays
  real(dp) :: ar_arr(MAX_CELLS)
  real(dp) :: avg_ar, min_ar, max_ar, std_ar, s1, s2
  integer  :: hist(5)

  read(*,*) mode

  select case (mode)

  !=========================================================================
  ! MODE 1 : Single triangle
  !=========================================================================
  case (1)
    read(*,*) vx(1),vy(1), vx(2),vy(2), vx(3),vy(3)

    ! Edge lengths
    edges(1) = sqrt((vx(2)-vx(1))**2 + (vy(2)-vy(1))**2)
    edges(2) = sqrt((vx(3)-vx(2))**2 + (vy(3)-vy(2))**2)
    edges(3) = sqrt((vx(1)-vx(3))**2 + (vy(1)-vy(3))**2)

    perimeter = edges(1) + edges(2) + edges(3)
    area = 0.5_dp * abs((vx(2)-vx(1))*(vy(3)-vy(1)) - (vx(3)-vx(1))*(vy(2)-vy(1)))

    emax = max(edges(1), edges(2), edges(3))
    emin = min(edges(1), edges(2), edges(3))
    AR_edge = emax / max(emin, 1.0e-30_dp)

    ! Inradius
    inrad = 2.0_dp * area / max(perimeter, 1.0e-30_dp)

    ! Circumradius
    if (area > 0.0_dp) then
      circumR = (edges(1)*edges(2)*edges(3)) / (4.0_dp * area)
    else
      circumR = 0.0_dp
    end if

    ! Inradius-based AR (equilateral = 1)
    if (inrad > 0.0_dp) then
      AR_main = emax / (2.0_dp * SQRT3 * inrad)
    else
      AR_main = 1.0e6_dp
    end if

    call get_rating(AR_main, rating)

    write(*,'(A,I1)')      'MODE=', mode
    write(*,'(A)')         'MODE_NAME=Single Triangle'
    write(*,'(A,ES15.8)')  'EDGE_1=', edges(1)
    write(*,'(A,ES15.8)')  'EDGE_2=', edges(2)
    write(*,'(A,ES15.8)')  'EDGE_3=', edges(3)
    write(*,'(A,ES15.8)')  'PERIMETER=', perimeter
    write(*,'(A,ES15.8)')  'AREA=', area
    write(*,'(A,ES15.8)')  'INRADIUS=', inrad
    write(*,'(A,ES15.8)')  'CIRCUMRADIUS=', circumR
    write(*,'(A,F12.6)')   'EMAX=', emax
    write(*,'(A,F12.6)')   'EMIN=', emin
    write(*,'(A,F12.6)')   'AR_MAIN=', AR_main
    write(*,'(A,F12.6)')   'AR_EDGE=', AR_edge
    write(*,'(A,A)')       'RATING=', trim(rating)

    write(*,'(A)') 'VERTS_START'
    do i = 1, 3
      write(*,'(F12.6,A,F12.6)') vx(i), ',', vy(i)
    end do
    write(*,'(A)') 'VERTS_END'

  !=========================================================================
  ! MODE 2 : Single quad
  !=========================================================================
  case (2)
    read(*,*) vx(1),vy(1), vx(2),vy(2), vx(3),vy(3), vx(4),vy(4)

    ! Edge lengths (sequential)
    do i = 1, 4
      edges(i) = sqrt((vx(mod(i,4)+1)-vx(i))**2 + (vy(mod(i,4)+1)-vy(i))**2)
    end do

    emax = max(edges(1), edges(2), edges(3), edges(4))
    emin = min(edges(1), edges(2), edges(3), edges(4))
    AR_edge = emax / max(emin, 1.0e-30_dp)

    ! Opposite-edge-pair AR
    avg1 = (edges(1) + edges(3)) / 2.0_dp
    avg2 = (edges(2) + edges(4)) / 2.0_dp
    AR_main = max(avg1, avg2) / max(min(avg1, avg2), 1.0e-30_dp)

    ! Diagonals
    d1 = sqrt((vx(3)-vx(1))**2 + (vy(3)-vy(1))**2)
    d2 = sqrt((vx(4)-vx(2))**2 + (vy(4)-vy(2))**2)
    diag_ratio = max(d1,d2) / max(min(d1,d2), 1.0e-30_dp)

    ! Area (shoelace)
    area = 0.5_dp * abs( (vx(1)*vy(2)-vx(2)*vy(1)) + (vx(2)*vy(3)-vx(3)*vy(2)) &
                       + (vx(3)*vy(4)-vx(4)*vy(3)) + (vx(4)*vy(1)-vx(1)*vy(4)) )

    perimeter = edges(1)+edges(2)+edges(3)+edges(4)

    call get_rating(AR_main, rating)

    write(*,'(A,I1)')      'MODE=', mode
    write(*,'(A)')         'MODE_NAME=Single Quad'
    write(*,'(A,ES15.8)')  'EDGE_1=', edges(1)
    write(*,'(A,ES15.8)')  'EDGE_2=', edges(2)
    write(*,'(A,ES15.8)')  'EDGE_3=', edges(3)
    write(*,'(A,ES15.8)')  'EDGE_4=', edges(4)
    write(*,'(A,ES15.8)')  'PERIMETER=', perimeter
    write(*,'(A,ES15.8)')  'AREA=', area
    write(*,'(A,F12.6)')   'AVG_PAIR1=', avg1
    write(*,'(A,F12.6)')   'AVG_PAIR2=', avg2
    write(*,'(A,F12.6)')   'DIAG_1=', d1
    write(*,'(A,F12.6)')   'DIAG_2=', d2
    write(*,'(A,F12.6)')   'DIAG_RATIO=', diag_ratio
    write(*,'(A,F12.6)')   'EMAX=', emax
    write(*,'(A,F12.6)')   'EMIN=', emin
    write(*,'(A,F12.6)')   'AR_MAIN=', AR_main
    write(*,'(A,F12.6)')   'AR_EDGE=', AR_edge
    write(*,'(A,A)')       'RATING=', trim(rating)

    write(*,'(A)') 'VERTS_START'
    do i = 1, 4
      write(*,'(F12.6,A,F12.6)') vx(i), ',', vy(i)
    end do
    write(*,'(A)') 'VERTS_END'

  !=========================================================================
  ! MODE 3 : Batch
  !=========================================================================
  case (3)
    backspace(5)
    read(*,*) mode, N
    if (N < 1 .or. N > MAX_CELLS) then
      write(*,'(A)') 'ERROR=N must be 1-10000.'; stop
    end if

    hist = 0
    do i = 1, N
      read(*,*) nv
      backspace(5)

      if (nv == 3) then
        read(*,*) nv, vx(1),vy(1), vx(2),vy(2), vx(3),vy(3)
        edges(1) = sqrt((vx(2)-vx(1))**2+(vy(2)-vy(1))**2)
        edges(2) = sqrt((vx(3)-vx(2))**2+(vy(3)-vy(2))**2)
        edges(3) = sqrt((vx(1)-vx(3))**2+(vy(1)-vy(3))**2)
        perimeter = edges(1)+edges(2)+edges(3)
        area = 0.5_dp*abs((vx(2)-vx(1))*(vy(3)-vy(1))-(vx(3)-vx(1))*(vy(2)-vy(1)))
        emax = max(edges(1),edges(2),edges(3))
        inrad = 2.0_dp*area/max(perimeter,1.0e-30_dp)
        if (inrad > 0.0_dp) then
          AR_main = emax/(2.0_dp*SQRT3*inrad)
        else
          AR_main = 1.0e6_dp
        end if

      else if (nv == 4) then
        read(*,*) nv, vx(1),vy(1), vx(2),vy(2), vx(3),vy(3), vx(4),vy(4)
        do nv = 1, 4
          edges(nv) = sqrt((vx(mod(nv,4)+1)-vx(nv))**2+(vy(mod(nv,4)+1)-vy(nv))**2)
        end do
        avg1 = (edges(1)+edges(3))/2.0_dp
        avg2 = (edges(2)+edges(4))/2.0_dp
        AR_main = max(avg1,avg2)/max(min(avg1,avg2),1.0e-30_dp)

      else
        read(*,*)
        AR_main = 1.0e6_dp
      end if

      ar_arr(i) = AR_main

      if (AR_main < 1.5_dp) then;      hist(1)=hist(1)+1
      else if (AR_main < 3.0_dp) then;  hist(2)=hist(2)+1
      else if (AR_main < 5.0_dp) then;  hist(3)=hist(3)+1
      else if (AR_main < 10.0_dp) then; hist(4)=hist(4)+1
      else;                              hist(5)=hist(5)+1
      end if
    end do

    ! Stats
    s1=0; s2=0; min_ar=ar_arr(1); max_ar=ar_arr(1)
    do i=1,N
      s1=s1+ar_arr(i); s2=s2+ar_arr(i)**2
      if(ar_arr(i)<min_ar) min_ar=ar_arr(i)
      if(ar_arr(i)>max_ar) max_ar=ar_arr(i)
    end do
    avg_ar=s1/real(N,dp)
    std_ar=sqrt(max(0.0_dp, s2/real(N,dp)-avg_ar**2))

    write(*,'(A,I1)')    'MODE=', mode
    write(*,'(A)')       'MODE_NAME=Batch'
    write(*,'(A,I6)')    'NCELLS=', N
    write(*,'(A,F12.6)') 'AVG_AR=', avg_ar
    write(*,'(A,F12.6)') 'MIN_AR=', min_ar
    write(*,'(A,F12.6)') 'MAX_AR=', max_ar
    write(*,'(A,F12.6)') 'STD_AR=', std_ar

    write(*,'(A)') 'HIST_START'
    write(*,'(A,I6)') '1.0-1.5,', hist(1)
    write(*,'(A,I6)') '1.5-3.0,', hist(2)
    write(*,'(A,I6)') '3.0-5.0,', hist(3)
    write(*,'(A,I6)') '5.0-10,',  hist(4)
    write(*,'(A,I6)') '>10,',     hist(5)
    write(*,'(A)') 'HIST_END'

    write(*,'(A)') 'DATA_START'
    do i=1,N
      call get_rating(ar_arr(i), rating)
      write(*,'(I6,A,F12.6,A,A)') i, ',', ar_arr(i), ',', trim(rating)
    end do
    write(*,'(A)') 'DATA_END'

  case default
    write(*,'(A)') 'ERROR=Invalid mode (must be 1-3).'; stop
  end select

contains

  subroutine get_rating(ar, r)
    real(dp), intent(in) :: ar
    character(len=20), intent(out) :: r
    if (ar < 1.5_dp) then;       r = 'Excellent'
    else if (ar < 3.0_dp) then;  r = 'Good'
    else if (ar < 5.0_dp) then;  r = 'Acceptable'
    else if (ar < 10.0_dp) then; r = 'Poor'
    else;                         r = 'Bad'
    end if
  end subroutine

end program aspect_ratio


Solver Description

Calculate the aspect ratio for triangular, quadrilateral, tetrahedral, and hexahedral grid cells.

Key Numerical Methods & Architecture

  • Input Redirection: Reads parameters sequentially from standard input (`stdin`) using Fortran sequential read (`read(*,*)`), ensuring modular integration.
  • Modular Design: Formulated using pure mathematical routines, separation of equations from output formatting, and precise numerical solvers (e.g. bisection, Newton-Raphson).
  • Standard Compliant: Written in clean, standards-compliant Fortran 90 to ensure cross-compiler compatibility.

๐Ÿ› ๏ธ Local Compilation

To test this code on your machine, compile the source code file(s) using a standard Fortran compiler (e.g., `gfortran`).

Compilation Command:

gfortran -O3 aspect_ratio.f90 -o aspect_ratio

Execution Command:

Execute the program by feeding the sample input file into the program using stdin redirection:

aspect_ratio < input.txt

๐Ÿ“ฅ Downloads & Local Files

Preview of the required input file (input.txt):

! Element type (1=Triangle, 2=Quadrilateral)\nCorner coordinates (x1 y1 x2 y2 x3 y3)
1
! Parameter 2
0.0 0.0 1.0 0.0 0.5 0.866