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 liquidpump(node1,node2,nodem,nelem,
20     &     nactdog,identity,ielprop,prop,iflag,v,xflow,f,
21     &     nodef,idirf,df,rho,g,co,numf,mi,ttime,time,
22     &     iaxial,iplausi)
23!
24!     pump for incompressible media
25!
26      implicit none
27!
28      logical identity
29!
30      integer nelem,nactdog(0:3,*),node1,node2,nodem,
31     &     ielprop(*),nodef(*),idirf(*),index,iflag,
32     &     inv,id,numf,npu,i,mi(*),iaxial,iplausi
33!
34      real*8 prop(*),v(0:mi(2),*),xflow,f,df(*),ttime,time,
35     &     p1,p2,rho,g(3),dg,z1,z2,co(3,*),
36     &     xpu(10),ypu(10),xxpu(10),yypu(10),dh
37!
38!
39!
40      numf=3
41!
42      if (iflag.eq.0) then
43         identity=.true.
44!
45         if(nactdog(2,node1).ne.0)then
46            identity=.false.
47         elseif(nactdog(2,node2).ne.0)then
48            identity=.false.
49         elseif(nactdog(1,nodem).ne.0)then
50            identity=.false.
51         endif
52!
53      elseif((iflag.eq.1).or.(iflag.eq.2)) then
54         if(iflag.eq.1) then
55            if(v(1,nodem).ne.0.d0) then
56               xflow=v(1,nodem)
57               return
58            endif
59         endif
60!
61         index=ielprop(nelem)
62!
63         npu=nint(prop(index+1))
64         do i=1,npu
65            xpu(i)=prop(index+2*i)
66            ypu(i)=prop(index+2*i+1)
67         enddo
68!
69         p1=v(2,node1)
70         p2=v(2,node2)
71!
72         z1=-g(1)*co(1,node1)-g(2)*co(2,node1)-g(3)*co(3,node1)
73         z2=-g(1)*co(1,node2)-g(2)*co(2,node2)-g(3)*co(3,node2)
74!
75         if(iflag.eq.2) then
76            xflow=v(1,nodem)*iaxial
77            if(xflow.ge.0.d0) then
78               inv=1
79            else
80               inv=-1
81            endif
82            nodef(1)=node1
83            nodef(2)=nodem
84            nodef(3)=node2
85            idirf(1)=2
86            idirf(2)=1
87            idirf(3)=2
88         endif
89!
90         dg=dsqrt(g(1)*g(1)+g(2)*g(2)+g(3)*g(3))
91!
92         if(iflag.eq.1) then
93            dh=(z2-z1+(p2-p1)/rho)/dg
94!
95!           reverting the order in xpu and ypu and storing the
96!           result in xxpu and yypu
97!
98            do i=1,npu
99               xxpu(i)=xpu(npu+1-i)
100               yypu(i)=ypu(npu+1-i)
101            enddo
102            call ident(yypu,dh,npu,id)
103            if(id.eq.0) then
104               xflow=xxpu(1)
105            elseif(id.eq.npu) then
106               xflow=0.d0
107            else
108               xflow=xxpu(id)+(xxpu(id+1)-xxpu(id))*(dh-yypu(id))/
109     &               (yypu(id+1)-yypu(id))
110            endif
111         else
112            df(1)=1.d0/rho
113            df(3)=-df(1)
114            xflow=xflow/rho
115            call ident(xpu,xflow,npu,id)
116            if(id.eq.0) then
117               if(xflow.ge.0.d0) then
118                  f=z1-z2+(p1-p2)/rho+dg*ypu(1)
119                  df(2)=0.d0
120               else
121                  df(2)=-1.d10
122                  f=z1-z2+(p1-p2)/rho+dg*(ypu(1)+xflow*df(2))
123                  df(2)=df(2)*dg/rho
124               endif
125            elseif(id.eq.npu) then
126               df(2)=-1.d10
127               f=z1-z2+(p1-p2)/rho+dg*(ypu(npu)+df(2)*(xflow-xpu(npu)))
128               df(2)=df(2)*dg/rho
129            else
130               df(2)=(ypu(id+1)-ypu(id))/(xpu(id+1)-xpu(id))
131               f=z1-z2+(p1-p2)/rho+dg*(ypu(id)+(xflow-xpu(id))*df(2))
132               df(2)=df(2)*dg/rho
133            endif
134         endif
135!
136      endif
137!
138      xflow=xflow/iaxial
139      df(2)=df(2)*iaxial
140!
141      return
142      end
143