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 mafillprhs(nk,kon,ipkon,lakon,ipompc,nodempc,
20     &     coefmpc,nmpc,b1,nactdoh,mi,v,theta1,nea,neb,dtimef,ipvar,var,
21     &     compressible)
22!
23!     filling the rhs b of the pressure equations (step 2)
24!
25      implicit none
26!
27      character*8 lakon(*)
28!
29      integer kon(*),ipompc(*),nodempc(3,*),nactdoh(*),compressible,
30     &     mi(*),ipkon(*),nea,neb,ipvar(*),nk,nmpc,i,j,
31     &     id,ist,index,jdof1,node,indexe,nope
32!
33      real*8 coefmpc(*),b1(nk,0:mi(2)),v(nk,0:mi(2)),ff(8),theta1,
34     &     var(*),dtimef
35!
36      do i=nea,neb
37!
38        indexe=ipkon(i)
39        if(lakon(i)(4:4).eq.'8') then
40          nope=8
41        elseif(lakon(i)(4:4).eq.'4') then
42          nope=4
43        elseif(lakon(i)(4:4).eq.'6') then
44          nope=6
45        else
46          cycle
47        endif
48!
49        call e_c3d_prhs(nk,kon(indexe+1),lakon(i),ff,i,v,dtimef,mi(1),
50     &       theta1,ipvar,var)
51!
52        if(compressible.eq.1) then
53!
54!         compressible fluid: dof = node number
55!
56          do j=1,nope
57            node=kon(indexe+j)
58            b1(node,4)=b1(node,4)+ff(j)
59          enddo
60        else
61!
62!     compressible fluid: dofs may have been removed due to
63!     SPC's and/or MPC's
64!
65          do j=1,nope
66!
67            node=kon(indexe+j)
68            jdof1=nactdoh(node)
69!
70!     inclusion of ff
71!
72            if(jdof1.le.0) then
73              if(nmpc.ne.0) then
74                if(jdof1.ne.2*(jdof1/2)) then
75                  id=(-jdof1+1)/2
76                  ist=ipompc(id)
77                  index=nodempc(3,ist)
78                  if(index.eq.0) cycle
79                  do
80                    jdof1=nactdoh(nodempc(1,index))
81                    if(jdof1.gt.0) then
82                      b1(jdof1,4)=b1(jdof1,4)
83     &                     -coefmpc(index)*ff(j)
84     &                     /coefmpc(ist)
85                    endif
86                    index=nodempc(3,index)
87                    if(index.eq.0) exit
88                  enddo
89                endif
90              endif
91              cycle
92            endif
93            b1(jdof1,4)=b1(jdof1,4)+ff(j)
94!
95          enddo
96        endif
97      enddo
98c      write(*,*) 'mafillprhs '
99c      do i=1,nk
100c        write(*,*) i,b1(i,4)
101c      enddo
102!
103      return
104      end
105
106