๐ป 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 Orthogonal Quality
Core Numerical Engine in Fortran 90 โข 28 total downloads
ortho_quality.f90
! =========================================================================
! Source File: 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
Solver Description
Verify mesh orthogonal quality and face deviation angles for finite volume solvers.
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 ortho_quality.f90 -o ortho_quality
Execution Command:
Execute the program by feeding the sample input file into the program using stdin redirection:
ortho_quality < 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
1
! Parameter 2
0.0 0.0 1.0 0.0 0.5 0.866