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 shape9q(xi,et,xl,xsj,xs,shp,iflag) 20! 21! shape functions and derivatives for a 9-node quadratic 22! isoparametric quadrilateral element. -1<=xi,et<=1 23! 24! iflag=1: calculate only the value of the shape functions 25! iflag=2: calculate the value of the shape functions, 26! their derivatives w.r.t. the local coordinates 27! and the Jacobian vector (local normal to the 28! surface) 29! iflag=3: calculate the value of the shape functions, the 30! value of their derivatives w.r.t. the global 31! coordinates and the Jacobian vector (local normal 32! to the surface) 33! iflag=4: calculate the value of the shape functions, the 34! value of their 1st and 2nd order derivatives 35! w.r.t. the local coordinates, the Jacobian vector 36! (local normal to the surface) 37! 38 implicit none 39! 40 integer i,j,k,iflag 41! 42 real*8 shp(7,9),xs(3,7),xsi(2,3),xl(3,9),sh(3),xsj(3),xi,et, 43 & fxi1,fxi2,fxi3,fet1,fet2,fet3,dfxi1,dfxi2,dfxi3,dfet1,dfet2, 44 & dfet3,ddfxi1,ddfxi2,ddfxi3,ddfet1,ddfet2,ddfet3 45! 46! 47! 48! shape functions and their glocal derivatives for an element 49! described with two local parameters and three global ones. 50! 51! shape functions in one dimension 52! 53 fxi1=xi*(xi-1.d0)/2.d0 54 fxi2=(1.d0-xi)*(1.d0+xi) 55 fxi3=xi*(xi+1.d0)/2.d0 56! 57 fet1=et*(et-1.d0)/2.d0 58 fet2=(1.d0-et)*(1.d0+et) 59 fet3=et*(et+1.d0)/2.d0 60! 61! shape functions 62! 63 shp(4,1)=fxi1*fet1 64 shp(4,2)=fxi3*fet1 65 shp(4,3)=fxi3*fet3 66 shp(4,4)=fxi1*fet3 67 shp(4,5)=fxi2*fet1 68 shp(4,6)=fxi3*fet2 69 shp(4,7)=fxi2*fet3 70 shp(4,8)=fxi1*fet2 71 shp(4,9)=fxi2*fet2 72! 73 if(iflag.eq.1) return 74! 75! derivative of the shape functions in one dimension 76! 77 dfxi1=(2.d0*xi-1.d0)/2.d0 78 dfxi2=-2.d0*xi 79 dfxi3=(2.d0*xi+1.d0)/2.d0 80! 81 dfet1=(2.d0*et-1.d0)/2.d0 82 dfet2=-2.d0*et 83 dfet3=(2.d0*et+1.d0)/2.d0 84! 85! local derivatives of the shape functions: xi-derivative 86! 87 shp(1,1)=dfxi1*fet1 88 shp(1,2)=dfxi3*fet1 89 shp(1,3)=dfxi3*fet3 90 shp(1,4)=dfxi1*fet3 91 shp(1,5)=dfxi2*fet1 92 shp(1,6)=dfxi3*fet2 93 shp(1,7)=dfxi2*fet3 94 shp(1,8)=dfxi1*fet2 95 shp(1,9)=dfxi2*fet2 96! 97! local derivatives of the shape functions: eta-derivative 98! 99 shp(2,1)=fxi1*dfet1 100 shp(2,2)=fxi3*dfet1 101 shp(2,3)=fxi3*dfet3 102 shp(2,4)=fxi1*dfet3 103 shp(2,5)=fxi2*dfet1 104 shp(2,6)=fxi3*dfet2 105 shp(2,7)=fxi2*dfet3 106 shp(2,8)=fxi1*dfet2 107 shp(2,9)=fxi2*dfet2 108! 109! computation of the local derivative of the global coordinates 110! (xs) 111! 112 do i=1,3 113 do j=1,2 114 xs(i,j)=0.d0 115 do k=1,9 116 xs(i,j)=xs(i,j)+xl(i,k)*shp(j,k) 117 enddo 118 enddo 119 enddo 120! 121! computation of the jacobian vector 122! 123 xsj(1)=xs(2,1)*xs(3,2)-xs(3,1)*xs(2,2) 124 xsj(2)=xs(1,2)*xs(3,1)-xs(3,2)*xs(1,1) 125 xsj(3)=xs(1,1)*xs(2,2)-xs(2,1)*xs(1,2) 126! 127 if(iflag.eq.3) then 128! 129! computation of the global derivative of the local coordinates 130! (xsi) (inversion of xs) 131! 132 if(dabs(xsj(3)).gt.1.d-10) then 133 xsi(1,1)=xs(2,2)/xsj(3) 134 xsi(2,2)=xs(1,1)/xsj(3) 135 xsi(1,2)=-xs(1,2)/xsj(3) 136 xsi(2,1)=-xs(2,1)/xsj(3) 137 if(dabs(xsj(2)).gt.1.d-10) then 138 xsi(2,3)=xs(1,1)/(-xsj(2)) 139 xsi(1,3)=-xs(1,2)/(-xsj(2)) 140 elseif(dabs(xsj(1)).gt.1.d-10) then 141 xsi(2,3)=xs(2,1)/xsj(1) 142 xsi(1,3)=-xs(2,2)/xsj(1) 143 else 144 xsi(2,3)=0.d0 145 xsi(1,3)=0.d0 146 endif 147 elseif(dabs(xsj(2)).gt.1.d-10) then 148 xsi(1,1)=xs(3,2)/(-xsj(2)) 149 xsi(2,3)=xs(1,1)/(-xsj(2)) 150 xsi(1,3)=-xs(1,2)/(-xsj(2)) 151 xsi(2,1)=-xs(3,1)/(-xsj(2)) 152 if(dabs(xsj(1)).gt.1.d-10) then 153 xsi(1,2)=xs(3,2)/xsj(1) 154 xsi(2,2)=-xs(3,1)/xsj(1) 155 else 156 xsi(1,2)=0.d0 157 xsi(2,2)=0.d0 158 endif 159 else 160 xsi(1,2)=xs(3,2)/xsj(1) 161 xsi(2,3)=xs(2,1)/xsj(1) 162 xsi(1,3)=-xs(2,2)/xsj(1) 163 xsi(2,2)=-xs(3,1)/xsj(1) 164 xsi(1,1)=0.d0 165 xsi(2,1)=0.d0 166 endif 167! 168! computation of the global derivatives of the shape functions 169! 170 do k=1,9 171 do j=1,3 172 sh(j)=shp(1,k)*xsi(1,j)+shp(2,k)*xsi(2,j) 173 enddo 174 do j=1,3 175 shp(j,k)=sh(j) 176 enddo 177 enddo 178! 179 elseif(iflag.eq.4) then 180! 181! second derivative of the shape functions in one dimension 182! 183 ddfxi1=1.d0 184 ddfxi2=-2.d0 185 ddfxi3=1.d0 186! 187 ddfet1=1.d0 188 ddfet2=-2.d0 189 ddfet3=1.d0 190! 191! local 2nd order derivatives of the shape functions: xi,xi-derivative 192! 193 shp(5,1)=ddfxi1*fet1 194 shp(5,2)=ddfxi3*fet1 195 shp(5,3)=ddfxi3*fet3 196 shp(5,4)=ddfxi1*fet3 197 shp(5,5)=ddfxi2*fet1 198 shp(5,6)=ddfxi3*fet2 199 shp(5,7)=ddfxi2*fet3 200 shp(5,8)=ddfxi1*fet2 201 shp(5,9)=ddfxi2*fet2 202! 203! local 2nd order derivatives of the shape functions: xi,eta-derivative 204! 205 shp(6,1)=dfxi1*dfet1 206 shp(6,2)=dfxi3*dfet1 207 shp(6,3)=dfxi3*dfet3 208 shp(6,4)=dfxi1*dfet3 209 shp(6,5)=dfxi2*dfet1 210 shp(6,6)=dfxi3*dfet2 211 shp(6,7)=dfxi2*dfet3 212 shp(6,8)=dfxi1*dfet2 213 shp(6,9)=dfxi2*dfet2 214! 215! local 2nd order derivatives of the shape functions: eta,eta-derivative 216! 217 shp(7,1)=fxi1*ddfet1 218 shp(7,2)=fxi3*ddfet1 219 shp(7,3)=fxi3*ddfet3 220 shp(7,4)=fxi1*ddfet3 221 shp(7,5)=fxi2*ddfet1 222 shp(7,6)=fxi3*ddfet2 223 shp(7,7)=fxi2*ddfet3 224 shp(7,8)=fxi1*ddfet2 225 shp(7,9)=fxi2*ddfet2 226! 227! computation of the local 2nd derivatives of the global coordinates 228! (xs) 229! 230 do i=1,3 231 do j=5,7 232 xs(i,j)=0.d0 233 do k=1,9 234 xs(i,j)=xs(i,j)+xl(i,k)*shp(j,k) 235 enddo 236 enddo 237 enddo 238 endif 239! 240 return 241 end 242