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 presgradient(iponoel,inoel,sa,nk,shockcoef,
20     &     dtimef,ipkon,kon,lakon,vold,mi,nmethod,
21     &     nactdoh,neqp)
22!
23!     determining measure for the pressure gradient
24!
25!     Ref: The Finite Element Method for Fluid Dynamics,
26!     O.C. Zienkiewicz, R.L. Taylor & P. Nithiarasu
27!     6th edition (2006) ISBN 0 7506 6322 7
28!     p. 61
29!
30      implicit none
31!
32      character*8 lakon(*)
33!
34      integer iponoel(*),inoel(2,*),nk,i,j,k,index,indexe,nope,neqp,
35     &     ipkon(*),kon(*),node,ielem,mi(*),nmethod,
36     &     nactdoh(*)
37!
38      real*8 sa(*),shockcoef,dtimef,ca,sum,pa,
39     &     vold(0:mi(2),*),sumabs,contribution
40!
41      do i=1,nk
42        if(nactdoh(i).le.0) cycle
43        if(iponoel(i).le.0) cycle
44        j=nactdoh(i)
45!
46        sum=0.d0
47        sumabs=0.d0
48        pa=vold(4,i)
49        index=iponoel(i)
50!
51        do
52          ielem=inoel(1,index)
53          if(ipkon(ielem).lt.0) cycle
54          if(lakon(ielem)(1:1).ne.'F') cycle
55          if(lakon(ielem)(4:4).eq.'8') then
56            nope=8
57          elseif(lakon(ielem)(4:4).eq.'4') then
58            nope=4
59          elseif(lakon(ielem)(4:4).eq.'6') then
60            nope=6
61          endif
62          indexe=ipkon(ielem)
63          do k=1,nope
64            node=kon(indexe+k)
65            if(node.eq.i) cycle
66!
67!           connecting vector between node i and "node"
68!
69!
70            contribution=(pa-vold(4,node))
71            sum=sum+contribution
72            sumabs=sumabs+dabs(contribution)
73!
74          enddo
75          index=inoel(2,index)
76          if(index.eq.0) exit
77        enddo
78        if(sumabs.lt.1.d-10) then
79          sum=0.d0
80          sumabs=1.d0
81        endif
82        sa(j)=dabs(sum)/(sumabs*dtimef)
83      enddo
84!
85      if(nmethod.eq.4) then
86!
87!     transient compressible
88!
89        ca=shockcoef*dtimef
90        do i=1,neqp
91          sa(i)=ca*sa(i)
92        enddo
93      else
94!
95!     steady state compressible
96!
97        ca=shockcoef*dtimef
98        do i=1,neqp
99          sa(i)=ca*sa(i)
100        enddo
101      endif
102!
103      return
104      end
105