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