!============================================================================== ! 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(Jdmxx) 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