๐ป 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 Cell Skewness
Core Numerical Engine in Fortran 90 โข 24 total downloads
skewness.f90
! =========================================================================
! Source File: skewness.f90
! =========================================================================
!==============================================================================
! ThermoFluidCalc โ Calculator #26 : Mesh Skewness
!==============================================================================
! Physics : Cell skewness measures how far an element deviates from the
! ideal (equilateral / equiangular) shape.
!
! Angular skewness (any 2-D cell):
! skew = max( (theta_max - theta_eq)/(180 - theta_eq) ,
! (theta_eq - theta_min)/theta_eq )
! theta_eq = (n-2)*180/n (60 for tri, 90 for quad)
!
! Equivolume skewness (triangle):
! skew_vol = (A_opt - A_cell) / A_opt
! A_opt = (3*sqrt(3)/4) * R^2 (equilateral tri with same circumradius R)
!
! Quality rating:
! 0.00 - 0.25 Excellent
! 0.25 - 0.50 Good
! 0.50 - 0.75 Acceptable
! 0.75 - 0.90 Poor
! 0.90 - 0.97 Bad
! 0.97 - 1.00 Degenerate
!
! Reference : Gupta, ยง5.1.1
!
! Modes:
! 1 = Single triangle (2-D, 3 vertices)
! 2 = Single quad (2-D, 4 vertices)
! 3 = Batch (N mixed cells)
!
! Build:
! gfortran -O2 -o skewness skewness.f90
!==============================================================================
program skewness
implicit none
integer, parameter :: dp = selected_real_kind(15, 307)
real(dp), parameter :: PI = 3.141592653589793238_dp
real(dp), parameter :: DEG = 180.0_dp / PI
integer, parameter :: MAX_CELLS = 10000
integer :: mode, N, i, nv
real(dp) :: vx(4), vy(4)
real(dp) :: angles(4), theta_min, theta_max, theta_eq
real(dp) :: skew_ang, skew_vol, area, circumR
real(dp) :: a_side, b_side, c_side, A_opt
character(len=20) :: rating
! Batch arrays
real(dp) :: skew_arr(MAX_CELLS)
real(dp) :: avg_s, min_s, max_s, std_s, sum_s, sum2
integer :: hist(5) ! bins: [0,.25) [.25,.5) [.5,.75) [.75,.9) [.9,1]
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_angles_tri(vx, vy, angles, theta_min, theta_max)
theta_eq = 60.0_dp
skew_ang = max((theta_max - theta_eq)/(180.0_dp - theta_eq), &
(theta_eq - theta_min)/theta_eq)
skew_ang = max(0.0_dp, min(skew_ang, 1.0_dp))
! Area
area = 0.5_dp * abs((vx(2)-vx(1))*(vy(3)-vy(1)) - (vx(3)-vx(1))*(vy(2)-vy(1)))
! Side lengths
a_side = sqrt((vx(2)-vx(1))**2 + (vy(2)-vy(1))**2)
b_side = sqrt((vx(3)-vx(2))**2 + (vy(3)-vy(2))**2)
c_side = sqrt((vx(1)-vx(3))**2 + (vy(1)-vy(3))**2)
! Circumradius R = abc / (4*Area)
if (area > 0.0_dp) then
circumR = (a_side * b_side * c_side) / (4.0_dp * area)
else
circumR = 0.0_dp
end if
! Equivolume skewness
A_opt = (3.0_dp * sqrt(3.0_dp) / 4.0_dp) * circumR**2
if (A_opt > 0.0_dp) then
skew_vol = (A_opt - area) / A_opt
else
skew_vol = 1.0_dp
end if
skew_vol = max(0.0_dp, min(skew_vol, 1.0_dp))
call get_rating(skew_ang, rating)
write(*,'(A,I1)') 'MODE=', mode
write(*,'(A)') 'MODE_NAME=Single Triangle'
write(*,'(A,F10.4)') 'ANGLE_1=', angles(1)
write(*,'(A,F10.4)') 'ANGLE_2=', angles(2)
write(*,'(A,F10.4)') 'ANGLE_3=', angles(3)
write(*,'(A,F10.4)') 'THETA_MIN=', theta_min
write(*,'(A,F10.4)') 'THETA_MAX=', theta_max
write(*,'(A,F10.4)') 'THETA_EQ=', theta_eq
write(*,'(A,F10.6)') 'SKEW_ANGULAR=', skew_ang
write(*,'(A,F10.6)') 'SKEW_EQUIVOL=', skew_vol
write(*,'(A,ES15.8)') 'AREA=', area
write(*,'(A,ES15.8)') 'CIRCUMRADIUS=', circumR
write(*,'(A,ES15.8)') 'A_OPT=', A_opt
write(*,'(A,A)') 'RATING=', trim(rating)
! Vertex data for drawing
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)
nv = 4
call compute_angles_quad(vx, vy, angles, theta_min, theta_max)
theta_eq = 90.0_dp
skew_ang = max((theta_max - theta_eq)/(180.0_dp - theta_eq), &
(theta_eq - theta_min)/theta_eq)
skew_ang = max(0.0_dp, min(skew_ang, 1.0_dp))
! Area via 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)) )
call get_rating(skew_ang, rating)
write(*,'(A,I1)') 'MODE=', mode
write(*,'(A)') 'MODE_NAME=Single Quad'
write(*,'(A,F10.4)') 'ANGLE_1=', angles(1)
write(*,'(A,F10.4)') 'ANGLE_2=', angles(2)
write(*,'(A,F10.4)') 'ANGLE_3=', angles(3)
write(*,'(A,F10.4)') 'ANGLE_4=', angles(4)
write(*,'(A,F10.4)') 'THETA_MIN=', theta_min
write(*,'(A,F10.4)') 'THETA_MAX=', theta_max
write(*,'(A,F10.4)') 'THETA_EQ=', theta_eq
write(*,'(A,F10.6)') 'SKEW_ANGULAR=', skew_ang
write(*,'(A,ES15.8)') 'AREA=', area
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)
call compute_angles_tri(vx, vy, angles, theta_min, theta_max)
theta_eq = 60.0_dp
else if (nv == 4) then
read(*,*) nv, vx(1),vy(1), vx(2),vy(2), vx(3),vy(3), vx(4),vy(4)
call compute_angles_quad(vx, vy, angles, theta_min, theta_max)
theta_eq = 90.0_dp
else
read(*,*) ! skip line
skew_arr(i) = 1.0_dp
cycle
end if
skew_ang = max((theta_max - theta_eq)/(180.0_dp - theta_eq), &
(theta_eq - theta_min)/theta_eq)
skew_ang = max(0.0_dp, min(skew_ang, 1.0_dp))
skew_arr(i) = skew_ang
! Histogram bin
if (skew_ang < 0.25_dp) then
hist(1) = hist(1) + 1
else if (skew_ang < 0.50_dp) then
hist(2) = hist(2) + 1
else if (skew_ang < 0.75_dp) then
hist(3) = hist(3) + 1
else if (skew_ang < 0.90_dp) then
hist(4) = hist(4) + 1
else
hist(5) = hist(5) + 1
end if
end do
! Statistics
sum_s = 0.0_dp; sum2 = 0.0_dp
min_s = skew_arr(1); max_s = skew_arr(1)
do i = 1, N
sum_s = sum_s + skew_arr(i)
sum2 = sum2 + skew_arr(i)**2
if (skew_arr(i) < min_s) min_s = skew_arr(i)
if (skew_arr(i) > max_s) max_s = skew_arr(i)
end do
avg_s = sum_s / real(N, dp)
std_s = sqrt(max(0.0_dp, sum2/real(N,dp) - avg_s**2))
write(*,'(A,I1)') 'MODE=', mode
write(*,'(A)') 'MODE_NAME=Batch'
write(*,'(A,I6)') 'NCELLS=', N
write(*,'(A,F10.6)') 'AVG_SKEW=', avg_s
write(*,'(A,F10.6)') 'MIN_SKEW=', min_s
write(*,'(A,F10.6)') 'MAX_SKEW=', max_s
write(*,'(A,F10.6)') 'STD_SKEW=', std_s
write(*,'(A)') 'HIST_START'
write(*,'(A,I6)') '0.00-0.25,', hist(1)
write(*,'(A,I6)') '0.25-0.50,', hist(2)
write(*,'(A,I6)') '0.50-0.75,', hist(3)
write(*,'(A,I6)') '0.75-0.90,', hist(4)
write(*,'(A,I6)') '0.90-1.00,', hist(5)
write(*,'(A)') 'HIST_END'
! Per-cell data
write(*,'(A)') 'DATA_START'
do i = 1, N
call get_rating(skew_arr(i), rating)
write(*,'(I6,A,F10.6,A,A)') i, ',', skew_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 compute_angles_tri(x, y, ang, amin, amax)
real(dp), intent(in) :: x(4), y(4)
real(dp), intent(out) :: ang(4), amin, amax
real(dp) :: dx1,dy1,dx2,dy2,dot,cross,mag1,mag2
integer :: j, j1, j2
do j = 1, 3
j1 = mod(j, 3) + 1 ! next vertex
j2 = mod(j + 1, 3) + 1 ! prev vertex
dx1 = x(j1) - x(j); dy1 = y(j1) - y(j)
dx2 = x(j2) - x(j); dy2 = y(j2) - y(j)
mag1 = sqrt(dx1**2 + dy1**2)
mag2 = sqrt(dx2**2 + dy2**2)
if (mag1 < 1.0e-30_dp .or. mag2 < 1.0e-30_dp) then
ang(j) = 0.0_dp; cycle
end if
dot = dx1*dx2 + dy1*dy2
ang(j) = acos(max(-1.0_dp, min(1.0_dp, dot/(mag1*mag2)))) * DEG
end do
ang(4) = 0.0_dp
amin = min(ang(1), ang(2), ang(3))
amax = max(ang(1), ang(2), ang(3))
end subroutine
!------------------------------------------------------------------------
subroutine compute_angles_quad(x, y, ang, amin, amax)
real(dp), intent(in) :: x(4), y(4)
real(dp), intent(out) :: ang(4), amin, amax
real(dp) :: dx1,dy1,dx2,dy2,mag1,mag2,dot
integer :: j, jp, jm
do j = 1, 4
jm = mod(j + 2, 4) + 1 ! previous vertex
jp = mod(j, 4) + 1 ! next vertex
dx1 = x(jm) - x(j); dy1 = y(jm) - y(j)
dx2 = x(jp) - x(j); dy2 = y(jp) - y(j)
mag1 = sqrt(dx1**2 + dy1**2)
mag2 = sqrt(dx2**2 + dy2**2)
if (mag1 < 1.0e-30_dp .or. mag2 < 1.0e-30_dp) then
ang(j) = 0.0_dp; cycle
end if
dot = dx1*dx2 + dy1*dy2
ang(j) = acos(max(-1.0_dp, min(1.0_dp, dot/(mag1*mag2)))) * DEG
end do
amin = min(ang(1), ang(2), ang(3), ang(4))
amax = max(ang(1), ang(2), ang(3), ang(4))
end subroutine
!------------------------------------------------------------------------
subroutine get_rating(s, r)
real(dp), intent(in) :: s
character(len=20), intent(out) :: r
if (s < 0.25_dp) then
r = 'Excellent'
else if (s < 0.50_dp) then
r = 'Good'
else if (s < 0.75_dp) then
r = 'Acceptable'
else if (s < 0.90_dp) then
r = 'Poor'
else if (s < 0.97_dp) then
r = 'Bad'
else
r = 'Degenerate'
end if
end subroutine
end program skewness
Solver Description
Evaluate mesh angular and equivolume skewness to avoid numerical solver instability.
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 skewness.f90 -o skewness
Execution Command:
Execute the program by feeding the sample input file into the program using stdin redirection:
skewness < 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