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 e_c3d_v2rhs(konl,lakonl,bb,nelem,v,dtimef,mi, 20 & ipvar,var,nk) 21! 22! computation of the velocity element matrix and rhs for the element with 23! element with the topology in konl: step 3 (correction **) 24! 25! bb: rhs 26! 27 implicit none 28! 29 character*8 lakonl 30! 31 integer konl(8),nelem,i,j,k,nope,mint3d,mi(*),ipvar(*),index,nk 32! 33 real*8 shp(4,8),bb(3,8),xsjmod,vl(0:mi(2),8),v(nk,0:mi(2)), 34 & var(*),dtimef,xsj,ddpress(3) 35! 36 if(lakonl(4:4).eq.'4') then 37 nope=4 38 mint3d=1 39 elseif(lakonl(4:4).eq.'6') then 40 nope=6 41 mint3d=2 42 elseif(lakonl(4:5).eq.'8R') then 43 nope=8 44 mint3d=1 45 elseif(lakonl(4:4).eq.'8') then 46 nope=8 47 mint3d=8 48 endif 49! 50! initialisation of the rhs 51! 52 do i=1,nope 53 do j=1,3 54 bb(j,i)=0.d0 55 enddo 56 enddo 57! 58! change in pressure 59! 60 do i=1,nope 61 vl(4,i)=v(konl(i),4) 62 enddo 63! 64! computation of the matrix: loop over the Gauss points 65! 66 index=ipvar(nelem) 67 do k=1,mint3d 68! 69! copying the shape functions, their derivatives and the 70! Jacobian determinant from field var 71! 72 do j=1,nope 73 do i=1,4 74 index=index+1 75 shp(i,j)=var(index) 76 enddo 77 enddo 78 index=index+1 79 xsj=var(index) 80! 81 xsjmod=dtimef*xsj 82! 83 index=index+1 84! 85! only for the semi-implicit procedure: calculate ddpress; 86! 87 do j=1,3 88 ddpress(j)=0.d0 89 enddo 90 do i=1,nope 91 do j=1,3 92 ddpress(j)=ddpress(j)+shp(j,i)*vl(4,i) 93 enddo 94 enddo 95! 96! 3. contribution to the second part of the 97! momentum equation 98! 99 do j=1,nope 100 bb(1,j)=bb(1,j)-xsjmod*shp(4,j)*ddpress(1) 101 bb(2,j)=bb(2,j)-xsjmod*shp(4,j)*ddpress(2) 102 bb(3,j)=bb(3,j)-xsjmod*shp(4,j)*ddpress(3) 103 enddo 104! 105 enddo 106! 107c write(*,*) 'e_c3d_v2rhs' 108c do k=1,8 109c write(*,*) nelem,k,(bb(j,k),j=1,3) 110c enddo 111! 112 return 113 end 114