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