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 mafillp(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,nefa,nefb,iau6,xxicn,flux)
23!
24!     filling the lhs and rhs to calculate p
25!
26      implicit none
27!
28      character*8 lakonf(*)
29!
30      integer i,nef,indexf,ipnei(*),j,neifa(*),nefa,nefb,
31     &  neiel(*),iel,ifa,irow(*),ielfa(4,*),compressible,
32     &  ifabou(*),neq,jq(*),iel2,indexb,indexf2,
33     &  j2,neij(*),nzs,k,iau6(6,*)
34!
35      real*8 coef,vfa(0:7,*),volume(*),area(*),advfa(*),xlet(*),
36     &  cosa(*),ad(*),au(*),xle(*),xxn(3,*),ap(*),b(*),cosb(*),
37     &  hfa(3,*),gradpel(3,*),bp(*),xxi(3,*),xlen(*),bp_ifa,
38     &  xxicn(3,*),flux(*)
39!
40      do i=nefa,nefb
41         indexf=ipnei(i)
42         do j=1,ipnei(i+1)-indexf
43!
44!     diffusion
45!
46            indexf=indexf+1
47            ifa=neifa(indexf)
48            iel=neiel(indexf)
49            if(iel.ne.0) then
50               coef=vfa(5,ifa)*area(ifa)*advfa(ifa)/
51     &              (xlet(indexf)*cosb(indexf))
52               ad(i)=ad(i)+coef
53               if(i.gt.iel) au(iau6(j,i))=au(iau6(j,i))-coef
54            else
55               if(ielfa(2,ifa).lt.0) then
56                  indexb=-ielfa(2,ifa)
57                  if(((ifabou(indexb+1).eq.0).or.
58     &                (ifabou(indexb+2).eq.0).or.
59     &                ( ifabou(indexb+3).eq.0)).and.
60     &               (ifabou(indexb+4).ne.0)) then
61!
62!                    pressure given (only if not all velocity
63!                    components are given)
64!
65                     coef=vfa(5,ifa)*area(ifa)*advfa(ifa)/
66     &                    (xle(indexf)*cosb(indexf))
67                     ad(i)=ad(i)+coef
68                  endif
69               endif
70            endif
71!
72!           right hand side: sum of the flux
73!
74            b(i)=b(i)-flux(indexf)
75         enddo
76      enddo
77!
78      return
79      end
80