🔬
Solver Purpose & Physical Scope
Calculate the aspect ratio for triangular, quadrilateral, tetrahedral, and hexahedral grid cells.
📂 Discipline: Cfd
⚡ Precision: IEEE-754 64-bit Real(`real(8)`)
📥 Total Downloads: 367 times
📄 Source File:
aspect_ratio.f90
📁 calcul/CFD /
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
💻 How to Compile & Run Locally
1. Compilation (GNU Fortran / Intel oneAPI):
gfortran -O3 aspect_ratio.f90 -o aspect_ratio
2. Execution with input.txt redirection:
aspect_ratio < input.txt
📄 Sample input.txt File Structure
Sample Data:
1 0.0 0.0 1.0 0.0 0.5 0.866
Parameter Description:
Element type (1=Triangle, 2=Quadrilateral)\nCorner coordinates (x1 y1 x2 y2 x3 y3)