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 patch(iterms,node,sti,scpav,mi,kon,ipkon, 20 & ipoints,members,linpatch,co,lakon,iavflag) 21! 22! computes the smoothed nodal stresses for an element patch 23! 24! author: Sascha Merz 25! 26 implicit none 27! 28 integer i,j,nope,mint3d,indexe,ielem,kon(*),ipkon(*), 29 & iterms,k,iflag,nrhs,info,node,ipnt, 30 & mi(*),irow,ipoints,members(*),linpatch,ielidx,iavflag 31! 32 real*8 xl(3,20),co(3,*),shp(4,20), 33 & pgauss(3),pol(20),pntdist,w,scpav(6,*),xsj, 34 & sti(6,mi(1),*),xi,et,ze,rv1(ipoints),pp(ipoints,iterms), 35 & pdat(ipoints,6),pfit(6),pwrk(iterms),pre(ipoints,iterms), 36 & z(ipoints,ipoints) 37! 38 real*8 tmpstr(6),gauss3d5e(3,4) 39! 40 character*8 lakon(*) 41! 42 logical matu,matv 43! 44 include 'gauss.f' 45! 46! initialize 47! 48! iavflag: if 1, then build average of integration point values 49! 50 if(iavflag.eq.0) then 51 do i=1,6 52 pfit(i)=0.d0 53 enddo 54 iflag=1 55! 56! irow: row number of the rectangular matrix of the overdetermined 57! system of equations 58! 59 irow=0 60! 61! loop over patch elements 62! linpatch: number of elements in patch 63! 64 do ielidx=1,linpatch 65! 66 ielem=members(ielidx) 67 if(lakon(ielem)(1:5).eq.'C3D20') then 68! 69! nope: nodes per element 70! 71 nope=20 72 if(lakon(ielem)(6:6).eq.'R') then 73 mint3d=8 74 else 75 mint3d=27 76 endif 77 elseif(lakon(ielem)(1:4).eq.'C3D8') then 78 nope=8 79 if(lakon(ielem)(5:5).eq.'R') then 80 mint3d=1 81 else 82 mint3d=8 83 endif 84 elseif(lakon(ielem)(1:5).eq.'C3D10') then 85 nope=10 86 mint3d=4 87 endif 88! 89 indexe=ipkon(ielem) 90! 91! coordinates of the nodes belonging to the element 92! 93 do j=1,nope 94 do k=1,3 95 xl(k,j)=co(k,kon(indexe+j)) 96 enddo 97 enddo 98! 99! loop over the integration points in one element 100! 101 do ipnt=1,mint3d 102 irow=irow+1 103 if(nope.eq.10) then 104 xi=gauss3d5(1,ipnt) 105 et=gauss3d5(2,ipnt) 106 ze=gauss3d5(3,ipnt) 107 call shape10tet(xi,et,ze,xl,xsj,shp,iflag) 108 elseif(nope.eq.20) then 109 if(mint3d.eq.8) then 110 xi=gauss3d2(1,ipnt) 111 et=gauss3d2(2,ipnt) 112 ze=gauss3d2(3,ipnt) 113 elseif(mint3d.eq.27) then 114 xi=gauss3d3(1,ipnt) 115 et=gauss3d3(2,ipnt) 116 ze=gauss3d3(3,ipnt) 117 endif 118 call shape20h(xi,et,ze,xl,xsj,shp,iflag) 119 elseif(nope.eq.8) then 120 if(mint3d.eq.1) then 121 xi=gauss3d1(1,ipnt) 122 et=gauss3d1(2,ipnt) 123 ze=gauss3d1(3,ipnt) 124 elseif(mint3d.eq.8) then 125 xi=gauss3d2(1,ipnt) 126 et=gauss3d2(2,ipnt) 127 ze=gauss3d2(3,ipnt) 128 endif 129 call shape8h(xi,et,ze,xl,xsj,shp,iflag) 130 endif 131! 132! in pgauss the global coordinates of the integration 133! point are saved. 134! 135 do j=1,3 136 pgauss(j)=0.d0 137 do k=1,nope 138 pgauss(j)=pgauss(j)+shp(4,k)*xl(j,k) 139 enddo 140! 141! the origin of the coordinate system is moved to the 142! evaluated node for higher numerical stability 143! 144 pgauss(j)=pgauss(j)-co(j,node) 145 enddo 146! 147! evaluate patch polynomial for the integration point 148! 149 pol(1)=1.d0 150 pol(2)=pgauss(1) 151 pol(3)=pgauss(2) 152 pol(4)=pgauss(3) 153 pol(5)=pgauss(1)*pgauss(2) 154 pol(6)=pgauss(1)*pgauss(3) 155 pol(7)=pgauss(2)*pgauss(3) 156 pol(8)=pgauss(1)*pgauss(1) 157 pol(9)=pgauss(2)*pgauss(2) 158 pol(10)=pgauss(3)*pgauss(3) 159 pol(11)=pgauss(1)*pgauss(2)*pgauss(3) 160 if(iterms.gt.11) then 161 pol(12)=pgauss(1)*pgauss(1)*pgauss(2) 162 pol(13)=pgauss(1)*pgauss(2)*pgauss(2) 163 pol(14)=pgauss(1)*pgauss(1)*pgauss(3) 164 pol(15)=pgauss(1)*pgauss(3)*pgauss(3) 165 pol(16)=pgauss(2)*pgauss(2)*pgauss(3) 166 pol(17)=pgauss(2)*pgauss(3)*pgauss(3) 167 if(iterms.gt.17) then 168 pol(18)=pgauss(1)*pgauss(1)*pgauss(1) 169 pol(19)=pgauss(2)*pgauss(2)*pgauss(2) 170 pol(20)=pgauss(3)*pgauss(3)*pgauss(3) 171 endif 172 endif 173! 174! weighting for integration point 175! 176 pntdist=dsqrt(pgauss(1)*pgauss(1) 177 & +pgauss(2)*pgauss(2) 178 & +pgauss(3)*pgauss(3)) 179 w=pntdist**(-1.5d0) 180! 181 do j=1,6 182 pdat(irow,j)=sti(j,ipnt,ielem)*w 183 enddo 184! 185 do j=1,iterms 186 pp(irow,j)=pol(j)*w 187 enddo 188 enddo 189 enddo 190! 191! using singular value decomposition for the least squares fit 192! 193 matu=.false. 194 matv=.true. 195 nrhs=6 196! 197 call hybsvd(ipoints,ipoints,ipoints,ipoints,ipoints, 198 & ipoints,iterms,pp,pwrk,matu,pp,matv, 199 & pre,z,pdat,nrhs,info,rv1) 200 if(info.ne.0) then 201 write(*,*) '*ERROR in patch: Bad conditioned matrix,', 202 & ' using average of sampling point values.' 203 iavflag=1 204 endif 205 endif 206! 207! matrix multiplication. only the first value of the 208! solution vector is needed. the singular values are manipulated 209! to increase the numerical stability 210! 211 if(iavflag.eq.0) then 212 do j=1,iterms 213 if(pwrk(j).lt.1.d-22) then 214 pwrk(j)=0.d0 215 else 216 pwrk(j)=1.d0/pwrk(j) 217 endif 218 enddo 219 do j=1,iterms 220 pre(1,j)=pre(1,j)*pwrk(j) 221 enddo 222 do j=1,nrhs 223 do k=1,iterms 224 pfit(j)=pfit(j)+pre(1,k)*pdat(k,j) 225 enddo 226 enddo 227! 228! pfit is an array containing the coefficients for the polynom 229! for the six stress components 230! 231! solution in the node 232! 233 do j=1,6 234 scpav(j,node)=scpav(j,node)+pfit(j) 235 enddo 236 endif 237! 238! if there are not enough elements to fit a polynomial, 239! build average value of the sampling point values 240! 241 if(iavflag.eq.1) then 242 do j=1,6 243 tmpstr(j)=0.d0 244 enddo 245 irow=0 246 do ielidx=1,linpatch 247 ielem=members(ielidx) 248 if(lakon(ielem)(1:5).eq.'C3D20') then 249 if(lakon(ielem)(6:6).eq.'R') then 250 mint3d=8 251 else 252 mint3d=27 253 endif 254 elseif(lakon(ielem)(1:5).eq.'C3D10') then 255 mint3d=4 256 elseif(lakon(ielem)(1:4).eq.'C3D8') then 257 if(lakon(ielem)(5:5).eq.'R') then 258 mint3d=1 259 else 260 mint3d=8 261 endif 262 endif 263 do ipnt=1,mint3d 264 irow=irow+1 265 do j=1,6 266 tmpstr(j)=tmpstr(j)+sti(j,ipnt,ielem) 267 enddo 268 enddo 269 enddo 270 do j=1,6 271 scpav(j,node)=scpav(j,node)+tmpstr(j)/irow 272 enddo 273 endif 274! 275 return 276 end 277