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