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