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_tee(df,pt1,Tt1,xflow1,xflow2,pt2, 21 &Tt2,A1,A2,zeta_fac,kappa,R,ider,iflag,zeta) 22! 23 implicit none 24! 25 integer ider,iflag 26! 27 real*8 28 &df(6), 29 &pt1, 30 &pt2, 31 &Tt1, 32 &Tt2, 33 &xflow1, 34 &xflow2, 35 &A1, 36 &A2, 37 &kappa, 38 &R, 39 &calc_residual_tee, 40! Accuracy of the numerical derivatives 41 &eps, 42! Step for the derivative 43 &h, 44 &f0, 45 &zeta_fac, 46 &zeta 47! 48 eps = 1.0e-8 49! 50 f0 = calc_residual_tee(pt1,Tt1,xflow1,xflow2,pt2, 51 &Tt2,A1,A2,zeta_fac,kappa,R,ider,iflag,zeta) 52! 53 h = eps*dabs(pt1) 54 if(h.eq.0)then 55 h = eps 56 endif 57 df(1) = (calc_residual_tee(pt1+h,Tt1,xflow1,xflow2,pt2, 58 &Tt2,A1,A2,zeta_fac,kappa,R,ider,iflag,zeta)-f0)/h 59! 60 h = eps*dabs(Tt1) 61 if(h.eq.0)then 62 h = eps 63 endif 64 df(2) = (calc_residual_tee(pt1,Tt1+h,xflow1,xflow2,pt2, 65 &Tt2,A1,A2,zeta_fac,kappa,R,ider,iflag,zeta)-f0)/h 66! 67 h = eps*dabs(xflow1) 68 if(h.eq.0)then 69 h = eps 70 endif 71 df(3) = (calc_residual_tee(pt1,Tt1,xflow1+h,xflow2,pt2, 72 &Tt2,A1,A2,zeta_fac,kappa,R,ider,iflag,zeta)-f0)/h 73! 74 h = eps*dabs(xflow2) 75 if(h.eq.0)then 76 h = eps 77 endif 78 df(4) = (calc_residual_tee(pt1,Tt1,xflow1,xflow2+h,pt2, 79 &Tt2,A1,A2,zeta_fac,kappa,R,ider,iflag,zeta)-f0)/h 80! 81 h = eps*dabs(pt2) 82 if(h.eq.0)then 83 h = eps 84 endif 85 df(5) = (calc_residual_tee(pt1,Tt1,xflow1,xflow2,pt2+h, 86 &Tt2,A1,A2,zeta_fac,kappa,R,ider,iflag,zeta)-f0)/h 87! 88 h = eps*dabs(Tt2) 89 if(h.eq.0)then 90 h = eps 91 endif 92 df(6) = (calc_residual_tee(pt1,Tt1,xflow1,xflow2,pt2, 93 &Tt2+h,A1,A2,zeta_fac,kappa,R,ider,iflag,zeta)-f0)/h 94! 95 return 96 end 97