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