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 extrapolate_d_v_simplec(ielfa,xrlfa,adv,advfa, 20 & hfa,icyclic,c,ifatie,vel,nef,volume, 21 & auv,ipnei,nfacea,nfaceb,ncfd) 22! 23! inter/extrapolation of volume/adv at the center of the elements 24! to the center of the faces 25! 26! inter/extrapolation of v* at the center of the elements 27! to the center of the faces; 28! 29 implicit none 30! 31 integer ielfa(4,*),iel1,iel2,iel3,i,j,icyclic,ifatie(*), 32 & nef,ipnei(*),indexf,nfacea,nfaceb,ncfd 33! 34 real*8 xrlfa(3,*),xl1,xl2,advfa(*),adv(*),vel(nef,0:7),hfa(3,*), 35 & c(3,3),volume(*),auv(*),a1,a2,a3 36! 37! 38! 39 do i=nfacea,nfaceb 40 iel1=ielfa(1,i) 41 xl1=xrlfa(1,i) 42 iel2=ielfa(2,i) 43 if(iel2.gt.0) then 44! 45! internal face 46! 47 xl2=xrlfa(2,i) 48! 49 a1=adv(iel1) 50 do indexf=ipnei(iel1)+1,ipnei(iel1+1) 51 a1=a1+auv(indexf) 52 enddo 53! 54 a2=adv(iel2) 55 do indexf=ipnei(iel2)+1,ipnei(iel2+1) 56 a2=a2+auv(indexf) 57 enddo 58! 59 advfa(i)=xl1*volume(iel1)/a1 60 & +xl2*volume(iel2)/a2 61! 62 if((icyclic.eq.0).or.(ifatie(i).eq.0)) then 63 do j=1,ncfd 64 hfa(j,i)=(xl1*vel(iel1,j) 65 & +xl2*vel(iel2,j)) 66 enddo 67 elseif(ifatie(i).gt.0) then 68 do j=1,ncfd 69 hfa(j,i)=(xl1*vel(iel1,j) 70 & +xl2*(c(j,1)*vel(iel2,1)+ 71 & c(j,2)*vel(iel2,2)+ 72 & c(j,3)*vel(iel2,3))) 73 enddo 74 else 75 do j=1,ncfd 76 hfa(j,i)=(xl1*vel(iel1,j) 77 & +xl2*(c(1,j)*vel(iel2,1)+ 78 & c(2,j)*vel(iel2,2)+ 79 & c(3,j)*vel(iel2,3))) 80 enddo 81 endif 82 elseif(ielfa(3,i).ne.0) then 83! 84! external face; linear extrapolation 85! 86 a1=adv(iel1) 87 do indexf=ipnei(iel1)+1,ipnei(iel1+1) 88 a1=a1+auv(indexf) 89 enddo 90! 91 iel3=abs(ielfa(3,i)) 92 a3=adv(iel3) 93 do indexf=ipnei(iel3)+1,ipnei(iel3+1) 94 a3=a3+auv(indexf) 95 enddo 96! 97 advfa(i)=xl1*volume(iel1)/a1 98 & +xrlfa(3,i)*volume(iel3)/a3 99! 100 do j=1,ncfd 101 hfa(j,i)=(xl1*vel(iel1,j)+xrlfa(3,i)*vel(iel3,j)) 102 enddo 103 else 104! 105! external face: constant extrapolation (only one adjacent 106! element layer) 107! 108 a1=adv(iel1) 109 do indexf=ipnei(iel1)+1,ipnei(iel1+1) 110 a1=a1+auv(indexf) 111 enddo 112! 113 advfa(i)=volume(iel1)/a1 114! 115 do j=1,ncfd 116 hfa(j,i)=vel(iel1,j) 117 enddo 118 endif 119 enddo 120! 121 return 122 end 123