!============================================================================== ! ThermoFluidCalc — Calculator #28 : Orthogonal Quality (Ortho-Skew) !============================================================================== ! Physics : Orthogonal quality measures how well the vector from the cell ! centroid to each face centroid aligns with the face outward ! normal. Perfect alignment (orthogonal grid) gives OQ = 1. ! ! For each face f: ! OQ_f = |dot(A_f, d_f)| / (|A_f| * |d_f|) = |cos(theta)| ! ! where ! A_f = outward face normal vector ! d_f = vector from cell centroid to face midpoint ! ! OQ_cell = min(OQ_f) over all faces of the cell ! ! Quality rating (0 = worst, 1 = best): ! 0.95 - 1.00 Excellent ! 0.70 - 0.95 Good ! 0.50 - 0.70 Acceptable ! 0.25 - 0.50 Poor ! 0.00 - 0.25 Bad ! ! Reference : Gupta, §5.1.1 ! ! Modes: ! 1 = Single triangle (3 vertices) ! 2 = Single quad (4 vertices) ! 3 = Batch (N mixed cells) ! ! Build: ! gfortran -O2 -o ortho_quality ortho_quality.f90 !============================================================================== program ortho_quality implicit none integer, parameter :: dp = selected_real_kind(15, 307) integer, parameter :: MAX_CELLS = 10000 integer :: mode, N, i, j, nv, jp real(dp) :: vx(4), vy(4) real(dp) :: cx, cy ! cell centroid real(dp) :: fmx, fmy ! face midpoint real(dp) :: nx_f, ny_f ! face outward normal real(dp) :: dx_f, dy_f ! centroid-to-face vector real(dp) :: mag_n, mag_d, dot_val, oq_f real(dp) :: oq_cell, oq_min_face real(dp) :: face_oq(4), face_nx(4), face_ny(4), face_mx(4), face_my(4) real(dp) :: face_dx(4), face_dy(4) character(len=20) :: rating ! Batch real(dp) :: oq_arr(MAX_CELLS) real(dp) :: avg_oq, min_oq, max_oq, std_oq, 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) nv = 3 call compute_oq(vx, vy, nv, cx, cy, face_oq, face_nx, face_ny, & face_mx, face_my, face_dx, face_dy, oq_cell) call get_rating(oq_cell, rating) write(*,'(A,I1)') 'MODE=', mode write(*,'(A)') 'MODE_NAME=Single Triangle' write(*,'(A,F12.6)') 'CENTROID_X=', cx write(*,'(A,F12.6)') 'CENTROID_Y=', cy write(*,'(A,F10.6)') 'OQ_CELL=', oq_cell write(*,'(A,A)') 'RATING=', trim(rating) write(*,'(A,I1)') 'NFACES=', nv write(*,'(A)') 'FACES_START' do j = 1, nv write(*,'(I2,5(A,F12.6))') j, & ',', face_nx(j), ',', face_ny(j), & ',', face_mx(j), ',', face_my(j), ',', face_oq(j) end do write(*,'(A)') 'FACES_END' write(*,'(A)') 'VERTS_START' do j = 1, nv write(*,'(F12.6,A,F12.6)') vx(j), ',', vy(j) 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) nv = 4 call compute_oq(vx, vy, nv, cx, cy, face_oq, face_nx, face_ny, & face_mx, face_my, face_dx, face_dy, oq_cell) call get_rating(oq_cell, rating) write(*,'(A,I1)') 'MODE=', mode write(*,'(A)') 'MODE_NAME=Single Quad' write(*,'(A,F12.6)') 'CENTROID_X=', cx write(*,'(A,F12.6)') 'CENTROID_Y=', cy write(*,'(A,F10.6)') 'OQ_CELL=', oq_cell write(*,'(A,A)') 'RATING=', trim(rating) write(*,'(A,I1)') 'NFACES=', nv write(*,'(A)') 'FACES_START' do j = 1, nv write(*,'(I2,5(A,F12.6))') j, & ',', face_nx(j), ',', face_ny(j), & ',', face_mx(j), ',', face_my(j), ',', face_oq(j) end do write(*,'(A)') 'FACES_END' write(*,'(A)') 'VERTS_START' do j = 1, nv write(*,'(F12.6,A,F12.6)') vx(j), ',', vy(j) 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) else if (nv == 4) then read(*,*) nv, vx(1),vy(1), vx(2),vy(2), vx(3),vy(3), vx(4),vy(4) else read(*,*) oq_arr(i) = 0.0_dp hist(5) = hist(5) + 1 cycle end if call compute_oq(vx, vy, nv, cx, cy, face_oq, face_nx, face_ny, & face_mx, face_my, face_dx, face_dy, oq_cell) oq_arr(i) = oq_cell if (oq_cell >= 0.95_dp) then; hist(1)=hist(1)+1 else if (oq_cell >= 0.70_dp) then; hist(2)=hist(2)+1 else if (oq_cell >= 0.50_dp) then; hist(3)=hist(3)+1 else if (oq_cell >= 0.25_dp) then; hist(4)=hist(4)+1 else; hist(5)=hist(5)+1 end if end do s1=0; s2=0; min_oq=oq_arr(1); max_oq=oq_arr(1) do i=1,N s1=s1+oq_arr(i); s2=s2+oq_arr(i)**2 if(oq_arr(i)max_oq) max_oq=oq_arr(i) end do avg_oq=s1/real(N,dp) std_oq=sqrt(max(0.0_dp, s2/real(N,dp)-avg_oq**2)) write(*,'(A,I1)') 'MODE=', mode write(*,'(A)') 'MODE_NAME=Batch' write(*,'(A,I6)') 'NCELLS=', N write(*,'(A,F10.6)') 'AVG_OQ=', avg_oq write(*,'(A,F10.6)') 'MIN_OQ=', min_oq write(*,'(A,F10.6)') 'MAX_OQ=', max_oq write(*,'(A,F10.6)') 'STD_OQ=', std_oq write(*,'(A)') 'HIST_START' write(*,'(A,I6)') '0.95-1.00,', hist(1) write(*,'(A,I6)') '0.70-0.95,', hist(2) write(*,'(A,I6)') '0.50-0.70,', hist(3) write(*,'(A,I6)') '0.25-0.50,', hist(4) write(*,'(A,I6)') '0.00-0.25,', hist(5) write(*,'(A)') 'HIST_END' write(*,'(A)') 'DATA_START' do i=1,N call get_rating(oq_arr(i), rating) write(*,'(I6,A,F10.6,A,A)') i, ',', oq_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 !------------------------------------------------------------------------ ! Compute orthogonal quality for a 2-D polygon cell (nv = 3 or 4) !------------------------------------------------------------------------ subroutine compute_oq(vx, vy, nv, cx, cy, foq, fnx, fny, fmx, fmy, fdx, fdy, oq_min) real(dp), intent(in) :: vx(4), vy(4) integer, intent(in) :: nv real(dp), intent(out) :: cx, cy, foq(4), fnx(4), fny(4) real(dp), intent(out) :: fmx(4), fmy(4), fdx(4), fdy(4), oq_min integer :: j, jp real(dp) :: ex, ey, mn, md, dv ! Centroid cx = 0.0_dp; cy = 0.0_dp do j = 1, nv cx = cx + vx(j); cy = cy + vy(j) end do cx = cx / real(nv, dp); cy = cy / real(nv, dp) oq_min = 1.0_dp do j = 1, nv jp = mod(j, nv) + 1 ! Face midpoint fmx(j) = 0.5_dp * (vx(j) + vx(jp)) fmy(j) = 0.5_dp * (vy(j) + vy(jp)) ! Edge vector ex = vx(jp) - vx(j) ey = vy(jp) - vy(j) ! Outward normal (rotate edge 90° clockwise, then check orientation) fnx(j) = ey fny(j) = -ex ! Ensure outward: normal should point away from centroid ! dot(normal, midpoint - centroid) > 0 for outward if (fnx(j)*(fmx(j)-cx) + fny(j)*(fmy(j)-cy) < 0.0_dp) then fnx(j) = -fnx(j) fny(j) = -fny(j) end if ! d = centroid to face midpoint fdx(j) = fmx(j) - cx fdy(j) = fmy(j) - cy ! OQ_f = |cos(theta)| mn = sqrt(fnx(j)**2 + fny(j)**2) md = sqrt(fdx(j)**2 + fdy(j)**2) if (mn > 1.0e-30_dp .and. md > 1.0e-30_dp) then dv = abs(fnx(j)*fdx(j) + fny(j)*fdy(j)) / (mn * md) foq(j) = min(dv, 1.0_dp) else foq(j) = 0.0_dp end if if (foq(j) < oq_min) oq_min = foq(j) end do ! Zero unused entries do j = nv+1, 4 foq(j) = 0.0_dp; fnx(j) = 0.0_dp; fny(j) = 0.0_dp fmx(j) = 0.0_dp; fmy(j) = 0.0_dp; fdx(j) = 0.0_dp; fdy(j) = 0.0_dp end do end subroutine !------------------------------------------------------------------------ subroutine get_rating(oq, r) real(dp), intent(in) :: oq character(len=20), intent(out) :: r if (oq >= 0.95_dp) then; r = 'Excellent' else if (oq >= 0.70_dp) then; r = 'Good' else if (oq >= 0.50_dp) then; r = 'Acceptable' else if (oq >= 0.25_dp) then; r = 'Poor' else; r = 'Bad' end if end subroutine end program ortho_quality