🔬
Solver Purpose & Physical Scope
Determine the determinant Jacobian ratio of finite element and finite volume grid mesh elements.
📂 Discipline: Cfd
⚡ Precision: IEEE-754 64-bit Real(`real(8)`)
📥 Total Downloads: 415 times
📄 Source File:
jacobian.f90
📁 calcul/CFD /
jacobian.f90
!==============================================================================
! ThermoFluidCalc — Calculator #29 : Jacobian of Transformation
!==============================================================================
! J = x_xi * y_eta - x_eta * y_xi
! Bilinear mapping (xi,eta in [0,1]):
! x = (1-xi)(1-eta)x1 + xi(1-eta)x2 + xi*eta*x3 + (1-xi)*eta*x4
! Reference : Gupta, §5.6
! Build: gfortran -O2 -o jacobian jacobian.f90
!==============================================================================
program jacobian
implicit none
integer, parameter :: dp = selected_real_kind(15,307)
integer, parameter :: MX=200
integer :: mode,N,M,i,j,idx
real(dp) :: vx(4),vy(4),xi,eta
real(dp) :: x_xi,x_eta,y_xi,y_eta,Jd
real(dp) :: Jc,Jcorner(4),Jmin,Jmax,Jrat
real(dp) :: xi_x,xi_y,eta_x,eta_y
real(dp) :: gx(0:MX,0:MX),gy(0:MX,0:MX)
real(dp) :: cx,cy,Jarr(MX*MX),s1,s2,avg,std,mn,mxx
integer :: neg
character(len=20) :: rating
real(dp) :: cxi(4),ceta(4)
data cxi /0.0_dp,1.0_dp,1.0_dp,0.0_dp/
data ceta /0.0_dp,0.0_dp,1.0_dp,1.0_dp/
read(*,*) mode
select case(mode)
!== MODE 1 ==
case(1)
read(*,*) vx(1),vy(1),vx(2),vy(2),vx(3),vy(3),vx(4),vy(4)
! Center
call eval_jac(vx,vy,0.5_dp,0.5_dp,x_xi,x_eta,y_xi,y_eta,Jc)
! Corners
do i=1,4
call eval_jac(vx,vy,cxi(i),ceta(i),x_xi,x_eta,y_xi,y_eta,Jcorner(i))
end do
Jmin=min(Jc,Jcorner(1),Jcorner(2),Jcorner(3),Jcorner(4))
Jmax=max(Jc,Jcorner(1),Jcorner(2),Jcorner(3),Jcorner(4))
if(abs(Jmax)>1e-30_dp) then; Jrat=Jmin/Jmax; else; Jrat=0; end if
call get_rat(Jrat,Jmin,rating)
! Full matrix at center
call eval_jac(vx,vy,0.5_dp,0.5_dp,x_xi,x_eta,y_xi,y_eta,Jd)
if(abs(Jd)>1e-30_dp) then
xi_x= y_eta/Jd; xi_y=-x_eta/Jd; eta_x=-y_xi/Jd; eta_y= x_xi/Jd
else
xi_x=0;xi_y=0;eta_x=0;eta_y=0
end if
write(*,'(A,I1)') 'MODE=',mode
write(*,'(A)') 'MODE_NAME=Single Quad'
write(*,'(A,ES15.8)') 'J_CENTER=',Jc
write(*,'(A,ES15.8)') 'J_MIN=',Jmin
write(*,'(A,ES15.8)') 'J_MAX=',Jmax
write(*,'(A,F12.6)') 'J_RATIO=',Jrat
write(*,'(A,A)') 'RATING=',trim(rating)
write(*,'(A,ES15.8)') 'X_XI=',x_xi
write(*,'(A,ES15.8)') 'X_ETA=',x_eta
write(*,'(A,ES15.8)') 'Y_XI=',y_xi
write(*,'(A,ES15.8)') 'Y_ETA=',y_eta
write(*,'(A,ES15.8)') 'XI_X=',xi_x
write(*,'(A,ES15.8)') 'XI_Y=',xi_y
write(*,'(A,ES15.8)') 'ETA_X=',eta_x
write(*,'(A,ES15.8)') 'ETA_Y=',eta_y
write(*,'(A)') 'CORNERS_START'
write(*,'(A,ES15.8)') 'SW,',Jcorner(1)
write(*,'(A,ES15.8)') 'SE,',Jcorner(2)
write(*,'(A,ES15.8)') 'NE,',Jcorner(3)
write(*,'(A,ES15.8)') 'NW,',Jcorner(4)
write(*,'(A,ES15.8)') 'CTR,',Jc
write(*,'(A)') 'CORNERS_END'
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 2 ==
case(2)
backspace(5); read(*,*) mode,N,M
if(N<1.or.N>MX.or.M<1.or.M>MX) then
write(*,'(A)') 'ERROR=N,M must be 1-200.'; stop
end if
do j=0,M; do i=0,N; read(*,*) gx(i,j),gy(i,j); end do; end do
idx=0; neg=0; s1=0; s2=0; mn=1e30_dp; mxx=-1e30_dp
write(*,'(A,I1)') 'MODE=',mode
write(*,'(A)') 'MODE_NAME=Grid Field'
write(*,'(A,I5)') 'NI=',N
write(*,'(A,I5)') 'NJ=',M
write(*,'(A)') 'DATA_START'
do j=0,M-1; do i=0,N-1
vx(1)=gx(i,j); vy(1)=gy(i,j)
vx(2)=gx(i+1,j); vy(2)=gy(i+1,j)
vx(3)=gx(i+1,j+1);vy(3)=gy(i+1,j+1)
vx(4)=gx(i,j+1); vy(4)=gy(i,j+1)
call eval_jac(vx,vy,0.5_dp,0.5_dp,x_xi,x_eta,y_xi,y_eta,Jd)
cx=0.25_dp*(vx(1)+vx(2)+vx(3)+vx(4))
cy=0.25_dp*(vy(1)+vy(2)+vy(3)+vy(4))
idx=idx+1; Jarr(idx)=Jd
if(Jd<0) neg=neg+1
s1=s1+Jd; s2=s2+Jd**2
if(Jd<mn) mn=Jd; if(Jd>mxx) mxx=Jd
write(*,'(I4,A,I4,A,F12.6,A,F12.6,A,ES15.8)') i,',',j,',',cx,',',cy,',',Jd
end do; end do
write(*,'(A)') 'DATA_END'
avg=s1/real(idx,dp); std=sqrt(max(0.0_dp,s2/real(idx,dp)-avg**2))
write(*,'(A,I6)') 'NCELLS=',idx
write(*,'(A,ES15.8)')'AVG_J=',avg
write(*,'(A,ES15.8)')'MIN_J=',mn
write(*,'(A,ES15.8)')'MAX_J=',mxx
write(*,'(A,ES15.8)')'STD_J=',std
write(*,'(A,I6)') 'NEG_COUNT=',neg
!== MODE 3 ==
case(3)
read(*,*) vx(1),vy(1),vx(2),vy(2),vx(3),vy(3),vx(4),vy(4),xi,eta
call eval_jac(vx,vy,xi,eta,x_xi,x_eta,y_xi,y_eta,Jd)
cx=(1-xi)*(1-eta)*vx(1)+xi*(1-eta)*vx(2)+xi*eta*vx(3)+(1-xi)*eta*vx(4)
cy=(1-xi)*(1-eta)*vy(1)+xi*(1-eta)*vy(2)+xi*eta*vy(3)+(1-xi)*eta*vy(4)
if(abs(Jd)>1e-30_dp) then
xi_x=y_eta/Jd; xi_y=-x_eta/Jd; eta_x=-y_xi/Jd; eta_y=x_xi/Jd
else
xi_x=0;xi_y=0;eta_x=0;eta_y=0
end if
write(*,'(A,I1)') 'MODE=',mode
write(*,'(A)') 'MODE_NAME=Probe Point'
write(*,'(A,F10.6)') 'XI=',xi
write(*,'(A,F10.6)') 'ETA=',eta
write(*,'(A,F12.6)') 'PHYS_X=',cx
write(*,'(A,F12.6)') 'PHYS_Y=',cy
write(*,'(A,ES15.8)') 'J=',Jd
write(*,'(A,ES15.8)') 'X_XI=',x_xi
write(*,'(A,ES15.8)') 'X_ETA=',x_eta
write(*,'(A,ES15.8)') 'Y_XI=',y_xi
write(*,'(A,ES15.8)') 'Y_ETA=',y_eta
write(*,'(A,ES15.8)') 'XI_X=',xi_x
write(*,'(A,ES15.8)') 'XI_Y=',xi_y
write(*,'(A,ES15.8)') 'ETA_X=',eta_x
write(*,'(A,ES15.8)') 'ETA_Y=',eta_y
if(Jd>0) then; write(*,'(A)') 'VALID=YES'; else; write(*,'(A)') 'VALID=NO'; end if
write(*,'(A)') 'VERTS_START'
do i=1,4; write(*,'(F12.6,A,F12.6)') vx(i),',',vy(i); end do
write(*,'(A)') 'VERTS_END'
case default
write(*,'(A)') 'ERROR=Invalid mode (1-3).'; stop
end select
contains
subroutine eval_jac(vx,vy,xi,eta,xxi,xet,yxi,yet,Jd)
real(dp),intent(in)::vx(4),vy(4),xi,eta
real(dp),intent(out)::xxi,xet,yxi,yet,Jd
xxi=-(1-eta)*vx(1)+(1-eta)*vx(2)+eta*vx(3)-eta*vx(4)
xet=-(1-xi)*vx(1)-xi*vx(2)+xi*vx(3)+(1-xi)*vx(4)
yxi=-(1-eta)*vy(1)+(1-eta)*vy(2)+eta*vy(3)-eta*vy(4)
yet=-(1-xi)*vy(1)-xi*vy(2)+xi*vy(3)+(1-xi)*vy(4)
Jd=xxi*yet-xet*yxi
end subroutine
subroutine get_rat(ratio,jmin,r)
real(dp),intent(in)::ratio,jmin
character(len=20),intent(out)::r
if(jmin<0) then; r='Folded/Invalid'
else if(ratio>0.5_dp) then; r='Excellent'
else if(ratio>0.3_dp) then; r='Good'
else if(ratio>0.1_dp) then; r='Acceptable'
else; r='Poor'
end if
end subroutine
end program jacobian
💻 How to Compile & Run Locally
1. Compilation (GNU Fortran / Intel oneAPI):
gfortran -O3 jacobian.f90 -o jacobian
2. Execution with input.txt redirection:
jacobian < input.txt
📄 Sample input.txt File Structure
Sample Data:
1 0.0 0.0 1.0 0.0 1.0 1.0 0.0 1.0
Parameter Description:
Mode (1=Single Quadrilateral, 2=Grid, 3=Single Point evaluation)\nCorner coordinates (x1 y1 x2 y2 x3 y3 x4 y4)