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 shape20h_ax(xi,et,ze,xl,xsj,shp,iflag) 20! 21! shape functions and derivatives for a 20-node quadratic 22! isoparametric brick element. -1<=xi,et,ze<=1 23! special case: axisymmetric elements 24! 25! iflag=1: calculate only the value of the shape functions 26! iflag=2: calculate the value of the shape functions and 27! the Jacobian determinant 28! iflag=3: calculate the value of the shape functions, the 29! value of their derivatives w.r.t. the global 30! coordinates and the Jacobian determinant 31! 32 implicit none 33! 34 integer j,k,iflag 35! 36 real*8 shp(4,20),xs(3,3),xsi(3,3),xl(3,20),shpe(4,20),dd, 37 & dd1,dd2,dd3 38! 39 real*8 xi,et,ze,xsj,omg,omh,omr,opg,oph,opr, 40 & tpgphpr,tmgphpr,tmgmhpr,tpgmhpr,tpgphmr,tmgphmr,tmgmhmr,tpgmhmr, 41 & omgopg,omhoph,omropr,omgmopg,omhmoph,omrmopr 42! 43! 44! 45! shape functions and their glocal derivatives 46! 47 omg=1.d0-xi 48 omh=1.d0-et 49 omr=1.d0-ze 50 opg=1.d0+xi 51 oph=1.d0+et 52 opr=1.d0+ze 53 tpgphpr=opg+oph+ze 54 tmgphpr=omg+oph+ze 55 tmgmhpr=omg+omh+ze 56 tpgmhpr=opg+omh+ze 57 tpgphmr=opg+oph-ze 58 tmgphmr=omg+oph-ze 59 tmgmhmr=omg+omh-ze 60 tpgmhmr=opg+omh-ze 61 omgopg=omg*opg/4.d0 62 omhoph=omh*oph/4.d0 63 omropr=omr*opr/4.d0 64 omgmopg=(omg-opg)/4.d0 65 omhmoph=(omh-oph)/4.d0 66 omrmopr=(omr-opr)/4.d0 67! 68! shape functions 69! 70 shp(4, 1)=-omg*omh*omr*tpgphpr/8.d0 71 shp(4, 2)=-opg*omh*omr*tmgphpr/8.d0 72 shp(4, 3)=-opg*oph*omr*tmgmhpr/8.d0 73 shp(4, 4)=-omg*oph*omr*tpgmhpr/8.d0 74 shp(4, 5)=-omg*omh*opr*tpgphmr/8.d0 75 shp(4, 6)=-opg*omh*opr*tmgphmr/8.d0 76 shp(4, 7)=-opg*oph*opr*tmgmhmr/8.d0 77 shp(4, 8)=-omg*oph*opr*tpgmhmr/8.d0 78 shp(4, 9)=omgopg*omh*omr 79 shp(4,10)=omhoph*opg*omr 80 shp(4,11)=omgopg*oph*omr 81 shp(4,12)=omhoph*omg*omr 82 shp(4,13)=omgopg*omh*opr 83 shp(4,14)=omhoph*opg*opr 84 shp(4,15)=omgopg*oph*opr 85 shp(4,16)=omhoph*omg*opr 86 shp(4,17)=omropr*omg*omh 87 shp(4,18)=omropr*opg*omh 88 shp(4,19)=omropr*opg*oph 89 shp(4,20)=omropr*omg*oph 90! 91 if(iflag.eq.1) return 92! 93! local derivatives of the shape functions: xi-derivative 94! 95 shpe(1, 1)=omh*omr*(tpgphpr-omg)/8.d0 96 shpe(1, 2)=(opg-tmgphpr)*omh*omr/8.d0 97 shpe(1, 3)=(opg-tmgmhpr)*oph*omr/8.d0 98 shpe(1, 4)=oph*omr*(tpgmhpr-omg)/8.d0 99 shpe(1, 5)=omh*opr*(tpgphmr-omg)/8.d0 100 shpe(1, 6)=(opg-tmgphmr)*omh*opr/8.d0 101 shpe(1, 7)=(opg-tmgmhmr)*oph*opr/8.d0 102 shpe(1, 8)=oph*opr*(tpgmhmr-omg)/8.d0 103 shpe(1, 9)=omgmopg*omh*omr 104 shpe(1,10)=omhoph*omr 105 shpe(1,11)=omgmopg*oph*omr 106 shpe(1,12)=-omhoph*omr 107 shpe(1,13)=omgmopg*omh*opr 108 shpe(1,14)=omhoph*opr 109 shpe(1,15)=omgmopg*oph*opr 110 shpe(1,16)=-omhoph*opr 111 shpe(1,17)=-omropr*omh 112 shpe(1,18)=omropr*omh 113 shpe(1,19)=omropr*oph 114 shpe(1,20)=-omropr*oph 115! 116! local derivatives of the shape functions: eta-derivative 117! 118 shpe(2, 1)=omg*omr*(tpgphpr-omh)/8.d0 119 shpe(2, 2)=opg*omr*(tmgphpr-omh)/8.d0 120 shpe(2, 3)=opg*(oph-tmgmhpr)*omr/8.d0 121 shpe(2, 4)=omg*(oph-tpgmhpr)*omr/8.d0 122 shpe(2, 5)=omg*opr*(tpgphmr-omh)/8.d0 123 shpe(2, 6)=opg*opr*(tmgphmr-omh)/8.d0 124 shpe(2, 7)=opg*(oph-tmgmhmr)*opr/8.d0 125 shpe(2, 8)=omg*(oph-tpgmhmr)*opr/8.d0 126 shpe(2, 9)=-omgopg*omr 127 shpe(2,10)=omhmoph*opg*omr 128 shpe(2,11)=omgopg*omr 129 shpe(2,12)=omhmoph*omg*omr 130 shpe(2,13)=-omgopg*opr 131 shpe(2,14)=omhmoph*opg*opr 132 shpe(2,15)=omgopg*opr 133 shpe(2,16)=omhmoph*omg*opr 134 shpe(2,17)=-omropr*omg 135 shpe(2,18)=-omropr*opg 136 shpe(2,19)=omropr*opg 137 shpe(2,20)=omropr*omg 138! 139! local derivatives of the shape functions: zeta-derivative 140! 141 shpe(3, 1)=omg*omh*(tpgphpr-omr)/8.d0 142 shpe(3, 2)=opg*omh*(tmgphpr-omr)/8.d0 143 shpe(3, 3)=opg*oph*(tmgmhpr-omr)/8.d0 144 shpe(3, 4)=omg*oph*(tpgmhpr-omr)/8.d0 145 shpe(3, 5)=omg*omh*(opr-tpgphmr)/8.d0 146 shpe(3, 6)=opg*omh*(opr-tmgphmr)/8.d0 147 shpe(3, 7)=opg*oph*(opr-tmgmhmr)/8.d0 148 shpe(3, 8)=omg*oph*(opr-tpgmhmr)/8.d0 149 shpe(3, 9)=-omgopg*omh 150 shpe(3,10)=-omhoph*opg 151 shpe(3,11)=-omgopg*oph 152 shpe(3,12)=-omhoph*omg 153 shpe(3,13)=omgopg*omh 154 shpe(3,14)=omhoph*opg 155 shpe(3,15)=omgopg*oph 156 shpe(3,16)=omhoph*omg 157 shpe(3,17)=omrmopr*omg*omh 158 shpe(3,18)=omrmopr*opg*omh 159 shpe(3,19)=omrmopr*opg*oph 160 shpe(3,20)=omrmopr*omg*oph 161! 162! computation of the local derivative of the global coordinates 163! (xs) 164! 165c do i=1,3 166c do j=1,3 167c xs(i,j)=0.d0 168c do k=1,20 169c xs(i,j)=xs(i,j)+xl(i,k)*shpe(j,k) 170c enddo 171c enddo 172c enddo 173 do j=1,3 174 xs(1,j)=xl(1,1)*(shpe(j,1)+shpe(j,5)) 175 & +xl(1,2)*(shpe(j,2)+shpe(j,6)) 176 & +xl(1,3)*(shpe(j,3)+shpe(j,7)) 177 & +xl(1,4)*(shpe(j,4)+shpe(j,8)) 178 & +xl(1,9)*(shpe(j,9)+shpe(j,13)) 179 & +xl(1,10)*(shpe(j,10)+shpe(j,14)) 180 & +xl(1,11)*(shpe(j,11)+shpe(j,15)) 181 & +xl(1,12)*(shpe(j,12)+shpe(j,16)) 182 & +xl(1,17)*shpe(j,17)+xl(1,18)*shpe(j,18) 183 & +xl(1,19)*shpe(j,19)+xl(1,20)*shpe(j,20) 184 xs(2,j)=xl(2,1)*(shpe(j,1)+shpe(j,5)+shpe(j,17)) 185 & +xl(2,2)*(shpe(j,2)+shpe(j,6)+shpe(j,18)) 186 & +xl(2,3)*(shpe(j,3)+shpe(j,7)+shpe(j,19)) 187 & +xl(2,4)*(shpe(j,4)+shpe(j,8)+shpe(j,20)) 188 & +xl(2,9)*(shpe(j,9)+shpe(j,13)) 189 & +xl(2,10)*(shpe(j,10)+shpe(j,14)) 190 & +xl(2,11)*(shpe(j,11)+shpe(j,15)) 191 & +xl(2,12)*(shpe(j,12)+shpe(j,16)) 192 xs(3,j)=xl(3,1)*(shpe(j,1)-shpe(j,5)) 193 & +xl(3,2)*(shpe(j,2)-shpe(j,6)) 194 & +xl(3,3)*(shpe(j,3)-shpe(j,7)) 195 & +xl(3,4)*(shpe(j,4)-shpe(j,8)) 196 & +xl(3,9)*(shpe(j,9)-shpe(j,13)) 197 & +xl(3,10)*(shpe(j,10)-shpe(j,14)) 198 & +xl(3,11)*(shpe(j,11)-shpe(j,15)) 199 & +xl(3,12)*(shpe(j,12)-shpe(j,16)) 200 enddo 201! 202! computation of the jacobian determinant 203! 204 dd1=xs(2,2)*xs(3,3)-xs(2,3)*xs(3,2) 205 dd2=xs(2,3)*xs(3,1)-xs(2,1)*xs(3,3) 206 dd3=xs(2,1)*xs(3,2)-xs(2,2)*xs(3,1) 207 xsj=xs(1,1)*dd1+xs(1,2)*dd2+xs(1,3)*dd3 208c xsj=xs(1,1)*(xs(2,2)*xs(3,3)-xs(2,3)*xs(3,2)) 209c & -xs(1,2)*(xs(2,1)*xs(3,3)-xs(2,3)*xs(3,1)) 210c & +xs(1,3)*(xs(2,1)*xs(3,2)-xs(2,2)*xs(3,1)) 211! 212 if(iflag.eq.2) return 213! 214 dd=1.d0/xsj 215! 216! computation of the global derivative of the local coordinates 217! (xsi) (inversion of xs) 218! 219 xsi(1,1)=dd1*dd 220 xsi(1,2)=(xs(1,3)*xs(3,2)-xs(1,2)*xs(3,3))*dd 221 xsi(1,3)=(xs(1,2)*xs(2,3)-xs(2,2)*xs(1,3))*dd 222 xsi(2,1)=dd2*dd 223 xsi(2,2)=(xs(1,1)*xs(3,3)-xs(3,1)*xs(1,3))*dd 224 xsi(2,3)=(xs(1,3)*xs(2,1)-xs(1,1)*xs(2,3))*dd 225 xsi(3,1)=dd3*dd 226 xsi(3,2)=(xs(1,2)*xs(3,1)-xs(1,1)*xs(3,2))*dd 227 xsi(3,3)=(xs(1,1)*xs(2,2)-xs(2,1)*xs(1,2))*dd 228c xsi(1,1)=(xs(2,2)*xs(3,3)-xs(3,2)*xs(2,3))*dd 229c xsi(1,2)=(xs(1,3)*xs(3,2)-xs(1,2)*xs(3,3))*dd 230c xsi(1,3)=(xs(1,2)*xs(2,3)-xs(2,2)*xs(1,3))*dd 231c xsi(2,1)=(xs(2,3)*xs(3,1)-xs(2,1)*xs(3,3))*dd 232c xsi(2,2)=(xs(1,1)*xs(3,3)-xs(3,1)*xs(1,3))*dd 233c xsi(2,3)=(xs(1,3)*xs(2,1)-xs(1,1)*xs(2,3))*dd 234c xsi(3,1)=(xs(2,1)*xs(3,2)-xs(3,1)*xs(2,2))*dd 235c xsi(3,2)=(xs(1,2)*xs(3,1)-xs(1,1)*xs(3,2))*dd 236c xsi(3,3)=(xs(1,1)*xs(2,2)-xs(2,1)*xs(1,2))*dd 237! 238! computation of the global derivatives of the shape functions 239! 240 do k=1,20 241 do j=1,3 242 shp(j,k)=shpe(1,k)*xsi(1,j)+shpe(2,k)*xsi(2,j) 243 & +shpe(3,k)*xsi(3,j) 244 enddo 245 enddo 246! 247 return 248 end 249