!============================================================================== ! 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)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