1! 2! CalculiX - A 3-dimensional finite element program 3! Copyright (C) 1998-2021 Guido Dhondt 4! 5! This program is free software; you can redistribute it and/or 6! modify it under the terms of the GNU General Public License as 7! published by the Free Software Foundation(version 2); 8! 9! This program is distributed in the hope that it will be useful, 10! but WITHOUT ANY WARRANTY; without even the implied warranty of 11! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 12! GNU General Public License for more details. 13! 14! You should have received a copy of the GNU General Public License 15! along with this program; if not, write to the Free Software 16! Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. 17! 18! author: Yannick Muller 19! 20 subroutine calc_ider_wye(df,pt1,Tt1,xflow1,xflow2,pt2, 21 &Tt2,ichan_num,A1,A2,A_s,dh1,dh2,alpha,zeta_fac,kappa,R,ider,iflag 22 &,zeta) 23! 24 implicit none 25! 26 integer ichan_num,ider,iflag 27! 28 real*8 29 &df(6), 30 &pt1, 31 &pt2, 32 &Tt1, 33 &Tt2, 34 &xflow1, 35 &xflow2, 36 &A1, 37 &A2, 38 &A_s, 39 &kappa, 40 &R, 41 &dh1, 42 &dh2, 43 &alpha, 44 &calc_residual_wye, 45 &eps, 46 &h, 47 &f0, 48 &zeta_fac, 49 &zeta 50! 51 eps = 1.0e-4 52! 53 f0 = calc_residual_wye(pt1,Tt1,xflow1,xflow2,pt2, 54 &Tt2,ichan_num,A1,A2,A_s,dh1,dh2,alpha,zeta_fac,kappa, 55 &R,ider,iflag,zeta) 56! 57 h = eps*dabs(pt1) 58 if(h.eq.0)then 59 h = eps 60 endif 61 df(1) = (calc_residual_wye(pt1+h,Tt1,xflow1,xflow2,pt2, 62 &Tt2,ichan_num,A1,A2,A_s,dh1,dh2,alpha,zeta_fac,kappa,R,ider, 63 &iflag,zeta)-f0)/h 64! 65 h = eps*dabs(Tt1) 66 if(h.eq.0)then 67 h = eps 68 endif 69 df(2) = (calc_residual_wye(pt1,Tt1+h,xflow1,xflow2,pt2, 70 &Tt2,ichan_num,A1,A2,A_s,dh1,dh2,alpha,zeta_fac,kappa,R,ider, 71 &iflag,zeta)-f0)/h 72! 73 h = eps*dabs(xflow1) 74 if(h.eq.0)then 75 h = eps 76 endif 77 df(3) = (calc_residual_wye(pt1,Tt1,xflow1+h,xflow2,pt2, 78 &Tt2,ichan_num,A1,A2,A_s,dh1,dh2,alpha,zeta_fac,kappa,R,ider, 79 &iflag,zeta)-f0)/h 80! 81 h = eps*dabs(xflow2) 82 if(h.eq.0)then 83 h = eps 84 endif 85 df(4) = (calc_residual_wye(pt1,Tt1,xflow1,xflow2+h,pt2, 86 &Tt2,ichan_num,A1,A2,A_s,dh1,dh2,alpha,zeta_fac,kappa,R,ider, 87 &iflag,zeta)-f0)/h 88! 89 h = eps*dabs(pt2) 90 if(h.eq.0)then 91 h = eps 92 endif 93 df(5) = (calc_residual_wye(pt1,Tt1,xflow1,xflow2,pt2+h, 94 &Tt2,ichan_num,A1,A2,A_s,dh1,dh2,alpha,zeta_fac,kappa,R,ider, 95 &iflag,zeta)-f0)/h 96! 97 h = eps*dabs(Tt2) 98 if(h.eq.0)then 99 h = eps 100 endif 101 df(6) = (calc_residual_wye(pt1,Tt1,xflow1,xflow2,pt2, 102 &Tt2+h,ichan_num,A1,A2,A_s,dh1,dh2,alpha,zeta_fac,kappa,R,ider, 103 &iflag,zeta)-f0)/h 104! 105 return 106 end 107