🔬
Solver Purpose & Physical Scope
Verify mesh orthogonal quality and face deviation angles for finite volume solvers.
📂 Discipline: Cfd
⚡ Precision: IEEE-754 64-bit Real(`real(8)`)
📥 Total Downloads: 442 times
📄 Source File:
ortho_quality.f90
📁 calcul/CFD /
ortho_quality.f90
!==============================================================================
! 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)<min_oq) min_oq=oq_arr(i)
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
💻 How to Compile & Run Locally
1. Compilation (GNU Fortran / Intel oneAPI):
gfortran -O3 ortho_quality.f90 -o ortho_quality
2. Execution with input.txt redirection:
ortho_quality < 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)