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! 10! This program is distributed in the hope that it will be useful, 11! but WITHOUT ANY WARRANTY; without even the implied warranty of 12! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 13! GNU General Public License for more details. 14! 15! You should have received a copy of the GNU General Public License 16! along with this program; if not, write to the Free Software 17! Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. 18! 19 subroutine calcgammav(ielfa,vel,gradvel,gamma,xlet, 20 & xxj,ipnei,betam,nef,flux,nfacea,nfaceb) 21! 22! determine gamma for the velocity: 23! upwind difference: gamma=0 24! central difference: gamma=1 25! 26 implicit none 27! 28 integer ielfa(4,*),i,j,indexf,ipnei(*),iel1,iel2,nef, 29 & nfacea,nfaceb 30! 31 real*8 vel(nef,0:7),gradvel(3,3,*),xxj(3,*),vud,vcd, 32 & gamma(*),phic,xlet(*),betam,flux(*),dvel1,dvel2 33! 34! 35! 36 do i=nfacea,nfaceb 37 iel2=ielfa(2,i) 38! 39! faces with only one neighbor need not be treated 40! 41 if(iel2.le.0) cycle 42 iel1=ielfa(1,i) 43 j=ielfa(4,i) 44 indexf=ipnei(iel1)+j 45! 46 dvel1=dsqrt(vel(iel1,1)**2+vel(iel1,2)**2+vel(iel1,3)**2) 47 dvel2=dsqrt(vel(iel2,1)**2+vel(iel2,2)**2+vel(iel2,3)**2) 48! 49 vcd=dvel2-dvel1 50! 51 if(dabs(vcd).lt.1.d-3*dvel1) vcd=0.d0 52! 53 if(flux(indexf).ge.0.d0) then 54! 55 vud=2.d0*xlet(indexf)*( 56 & (vel(iel1,1)*gradvel(1,1,iel1)+ 57 & vel(iel1,2)*gradvel(2,1,iel1)+ 58 & vel(iel1,3)*gradvel(3,1,iel1))*xxj(1,indexf)+ 59 & (vel(iel1,1)*gradvel(1,2,iel1)+ 60 & vel(iel1,2)*gradvel(2,2,iel1)+ 61 & vel(iel1,3)*gradvel(3,2,iel1))*xxj(2,indexf)+ 62 & (vel(iel1,1)*gradvel(1,3,iel1)+ 63 & vel(iel1,2)*gradvel(2,3,iel1)+ 64 & vel(iel1,3)*gradvel(3,3,iel1))*xxj(3,indexf)) 65 vcd=vcd*dvel1 66 else 67 vud=2.d0*xlet(indexf)*( 68 & (vel(iel2,1)*gradvel(1,1,iel2)+ 69 & vel(iel2,2)*gradvel(2,1,iel2)+ 70 & vel(iel2,3)*gradvel(3,1,iel2))*xxj(1,indexf)+ 71 & (vel(iel2,1)*gradvel(1,2,iel2)+ 72 & vel(iel2,2)*gradvel(2,2,iel2)+ 73 & vel(iel2,3)*gradvel(3,2,iel2))*xxj(2,indexf)+ 74 & (vel(iel2,1)*gradvel(1,3,iel2)+ 75 & vel(iel2,2)*gradvel(2,3,iel2)+ 76 & vel(iel2,3)*gradvel(3,3,iel2))*xxj(3,indexf)) 77 vcd=vcd*dvel2 78 endif 79! 80 if(dabs(vud).lt.1.d-20) then 81 gamma(i)=0.d0 82 cycle 83 endif 84! 85 phic=1.d0-vcd/vud 86! 87 if(phic.ge.1.d0) then 88 gamma(i)=0.d0 89 elseif(phic.le.0.d0) then 90 gamma(i)=0.d0 91 elseif(betam.le.phic) then 92 gamma(i)=1.d0 93 else 94 gamma(i)=phic/betam 95 endif 96 enddo 97! 98 return 99 end 100