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 extrapol_vel4(ielfa,vfa,vfap,gradvfa,rf,ifabou, 20 & ipnei,xxn,vel,xxi,xle,nef,nfacea,nfaceb,ncfd) 21! 22! inter/extrapolation of v at the center of the elements 23! to the center of the faces: taking the skewness of the 24! elements into account (using the gradient at the center 25! of the faces) 26! 27 implicit none 28! 29 integer ielfa(4,*),ifabou(*),ipnei(*),nef,nfacea,nfaceb,i,j, 30 & indexf,ipointer,iel1,iel2,ncfd 31! 32 real*8 vfa(0:7,*),vfap(0:7,*),gradvfa(3,3,*),rf(3,*),xxn(3,*), 33 & vel(nef,0:7),xxi(3,*),xle(*),dd 34! 35! 36! 37! Moukalled et al. p 279 38! 39 do i=nfacea,nfaceb 40 iel1=ielfa(1,i) 41 iel2=ielfa(2,i) 42 if(iel2.gt.0) then 43! 44! face between two elements 45! 46 do j=1,ncfd 47 vfa(j,i)=vfap(j,i)+gradvfa(j,1,i)*rf(1,i)+ 48 & gradvfa(j,2,i)*rf(2,i)+ 49 & gradvfa(j,3,i)*rf(3,i) 50 enddo 51 elseif(ielfa(3,i).gt.0) then 52! 53! boundary face; no zero gradient 54! 55 ipointer=-iel2 56! 57! x-direction 58! 59 if(ifabou(ipointer+1).gt.0) then 60 vfa(1,i)=vfap(1,i) 61 else 62 vfa(1,i)=vfap(1,i)+gradvfa(1,1,i)*rf(1,i)+ 63 & gradvfa(1,2,i)*rf(2,i)+ 64 & gradvfa(1,3,i)*rf(3,i) 65 endif 66! 67! y-direction 68! 69 if(ifabou(ipointer+2).gt.0) then 70 vfa(2,i)=vfap(2,i) 71 else 72 vfa(2,i)=vfap(2,i)+gradvfa(2,1,i)*rf(1,i)+ 73 & gradvfa(2,2,i)*rf(2,i)+ 74 & gradvfa(2,3,i)*rf(3,i) 75 endif 76! 77! z-direction 78! 79 if(ifabou(ipointer+3).gt.0) then 80 vfa(3,i)=vfap(3,i) 81 else 82 vfa(3,i)=vfap(3,i)+gradvfa(3,1,i)*rf(1,i)+ 83 & gradvfa(3,2,i)*rf(2,i)+ 84 & gradvfa(3,3,i)*rf(3,i) 85 endif 86! 87! correction for sliding boundary conditions 88! 89 if(ifabou(ipointer+5).lt.0) then 90 indexf=ipnei(iel1)+ielfa(4,i) 91 dd=vfa(1,i)*xxn(1,indexf)+ 92 & vfa(2,i)*xxn(2,indexf)+ 93 & vfa(3,i)*xxn(3,indexf) 94 do j=1,ncfd 95 vfa(j,i)=vfa(j,i)-dd*xxn(j,indexf) 96 enddo 97 endif 98 else 99! 100! boundary face; zero gradient 101! 102c indexf=ipnei(iel1)+ielfa(4,i) 103 do j=1,ncfd 104 vfa(j,i)=vel(iel1,j) 105c & +(gradvfa(j,1,i)*xxi(1,indexf)+ 106c & gradvfa(j,2,i)*xxi(2,indexf)+ 107c & gradvfa(j,3,i)*xxi(3,indexf))*xle(indexf) 108 enddo 109 endif 110 enddo 111! 112 return 113 end 114