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 mafillpcomp(nef,lakonf,ipnei,neifa,neiel,vfa,area, 20 & advfa,xlet,cosa,volume,au,ad,jq,irow,ap,ielfa,ifabou,xle, 21 & b,xxn,neq,nzs,hfa,gradpel,bp,xxi,neij, 22 & xlen,cosb,ielmatf,mi,a1,a2,a3,velo,veloo,dtimef,shcon, 23 & ntmat_,vel,nactdohinv,xrlfa,flux,nefa,nefb,iau6,xxicn, 24 & gamma,inlet) 25! 26! filling the lhs and rhs to calculate p 27! 28! bc: 29! outflow, p known: diffusion (subsonic) 30! p unknown: convection (supersonic) 31! inflow, p known: none 32! p unknown: convection (subsonic) 33! 34 implicit none 35! 36 character*8 lakonf(*) 37! 38 integer i,nef,indexf,ipnei(*),neifa(*),inlet(*), 39 & neiel(*),iel,ifa,irow(*),ielfa(4,*), 40 & ifabou(*),neq,jq(*),indexb, 41 & neij(*),nzs,imat,nefa,nefb,iau6(6,*), 42 & mi(*),ielmatf(mi(3),*),ntmat_,nactdohinv(*) 43! 44 real*8 coef,vfa(0:7,*),volume(*),area(*),advfa(*),xlet(*), 45 & cosa(*),ad(neq),au(nzs),xle(*),xxn(3,*),ap(*),b(neq),cosb(*), 46 & hfa(3,*),gradpel(3,*),bp(*),xxi(3,*),xlen(*),r, 47 & xflux,velo(nef,0:7),veloo(nef,0:7),dtimef, 48 & shcon(0:3,ntmat_,*),vel(nef,0:7),a1,a2,a3, 49 & xrlfa(3,*),gamma(*),flux(*),xxicn(3,*) 50! 51! 52! 53 do i=nefa,nefb 54 imat=ielmatf(1,i) 55 r=shcon(3,1,imat) 56 do indexf=ipnei(i)+1,ipnei(i+1) 57! 58! loop over all element faces 59! 60 ifa=neifa(indexf) 61 iel=neiel(indexf) 62 xflux=flux(indexf) 63! 64! convection (density correction) 65! 66 if(xflux.ge.0.d0) then 67! 68! outflowing flux 69! 70 if(iel.gt.0) then 71! 72! internal face 73! 74 ad(i)=ad(i)+xflux/(vfa(5,ifa)*r*vfa(0,ifa)) 75! 76 elseif(ielfa(3,ifa).le.0) then 77 indexb=-ielfa(2,ifa) 78 if(indexb.gt.0) then 79 if((ifabou(indexb+4).eq.0).and. 80 & (ifabou(indexb+5).eq.0)) then 81! 82! outflow, pressure not known 83! 84 ad(i)=ad(i)+xflux/(vfa(5,ifa)*r*vfa(0,ifa)) 85 endif 86 else 87! 88! outflow, pressure not known 89! 90 ad(i)=ad(i)+xflux/(vfa(5,ifa)*r*vfa(0,ifa)) 91 endif 92 endif 93 else 94! 95! inflowing flux 96! 97 if(iel.gt.0) then 98! 99! internal face 100! 101 au(indexf)=au(indexf)+xflux/(vfa(5,ifa)*r*vfa(0,ifa)) 102 else 103! 104! external face 105! 106 if(ielfa(2,ifa).lt.0) then 107 indexb=-ielfa(2,ifa) 108 if((inlet(ifa).ne.0).and. 109 & (ifabou(indexb+4).eq.0)) then 110! 111! inlet 112! and no pressure given (typical subsonic inlet) 113! 114 ad(i)=ad(i)+xflux/(vfa(5,ifa)*r*vfa(0,ifa)) 115 endif 116 endif 117 endif 118 endif 119! 120! diffusion (velocity correction) 121! 122 if(iel.gt.0) then 123! 124! internal face 125! 126 coef=vfa(5,ifa)*area(ifa)*advfa(ifa)/ 127 & xlet(indexf) 128 ad(i)=ad(i)+coef 129 au(indexf)=au(indexf)-coef 130 else 131! 132! external face 133! 134 if(ielfa(2,ifa).lt.0) then 135 indexb=-ielfa(2,ifa) 136 if((xflux.ge.0.d0).and. 137 & (ifabou(indexb+4).ne.0)) then 138! 139! not all velocity components known 140! pressure known (typical subsonic outlet) 141! 142 coef=vfa(5,ifa)*area(ifa)*advfa(ifa)/ 143 & xle(indexf) 144 ad(i)=ad(i)+coef 145 endif 146 endif 147 endif 148! 149! flux 150! 151 b(i)=b(i)-flux(indexf) 152 enddo 153! 154! transient term 155! 156 ad(i)=ad(i)+volume(i)/(r*vel(i,0)*dtimef) 157 b(i)=b(i)+(velo(i,5)-vel(i,5))*volume(i)/dtimef 158 enddo 159! 160 return 161 end 162