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 correctflux(nef,ipnei,neifa,neiel,flux,vfa,advfa,area, 20 & vel,alet,ielfa,ale,ifabou,nefa,nefb,xxnj,gradpcfa) 21! 22! correction of v due to the balance of mass 23! the correction is in normal direction to the face 24! 25 implicit none 26! 27 integer i,nef,indexf,ipnei(*),neifa(*),neiel(*),ielfa(4,*), 28 & iel,ifa,ifabou(*),indexb,nefa,nefb 29! 30 real*8 flux(*),vfa(0:7,*),advfa(*),area(*),vel(nef,0:7),alet(*), 31 & ale(*),xxnj(3,*),gradpcfa(3,*) 32! 33! 34! 35 do i=nefa,nefb 36c totflux=0.d0 37 do indexf=ipnei(i)+1,ipnei(i+1) 38 ifa=neifa(indexf) 39 iel=neiel(indexf) 40 if(iel.gt.0) then 41! 42! internal face 43! 44 flux(indexf)=flux(indexf)+vfa(5,ifa)*advfa(ifa)* 45 & ((vel(i,4)-vel(iel,4))*alet(indexf) 46 & -(gradpcfa(1,ifa)*xxnj(1,indexf)+ 47 & gradpcfa(2,ifa)*xxnj(2,indexf)+ 48 & gradpcfa(3,ifa)*xxnj(3,indexf))) 49 else 50 indexb=-ielfa(2,ifa) 51 if(indexb.gt.0) then 52 if(((ifabou(indexb+1).eq.0).or. 53 & (ifabou(indexb+2).eq.0).or. 54 & (ifabou(indexb+3).eq.0)).and. 55 & (ifabou(indexb+4).ne.0)) then 56! 57! external face with pressure boundary conditions 58! 59 flux(indexf)=flux(indexf) 60 & +vfa(5,ifa)*advfa(ifa)*(vel(i,4)*ale(indexf) 61 & -(gradpcfa(1,ifa)*xxnj(1,indexf)+ 62 & gradpcfa(2,ifa)*xxnj(2,indexf)+ 63 & gradpcfa(3,ifa)*xxnj(3,indexf))) 64 endif 65 endif 66 endif 67c totflux=totflux+flux(indexf) 68 enddo 69 enddo 70! 71 return 72 end 73