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 calcexternalwork(co,vold,istartset,iendset,ipkon, 20 & lakon,kon,ialset,nset,set,mi,externalfaces,nelemload, 21 & nload,xload,sideload,delexternalwork,voldprev) 22! 23! calculation and printout of the lift and drag forces 24! 25 implicit none 26! 27 character*8 lakonl,lakon(*) 28 character*20 sideload(*) 29 character*81 set(*),externalfaces 30! 31 integer konl(20),ifaceq(8,6),nelem,i,j,i1,j1, 32 & jj,ig,nope,nopes,igpres,id, 33 & mint2d,ifacet(6,4),ifacew(8,5),iflag,indexe,jface,istartset(*), 34 & iendset(*),ipkon(*),kon(*),iset,ialset(*),nset,ipos, 35 & mi(*),nelemload(2,*),nload,ipres 36! 37 real*8 co(3,*),xs2(3,7),pres,voldl2(3,8),xl2(3,8),xsj2(3), 38 & shp2(7,8),vold(0:mi(2),*),xi,et,weight,voldprev(0:mi(2),*), 39 & xlocal20(3,9,6),xlocal4(3,1,4),xlocal10(3,3,4),xlocal6(3,1,5), 40 & xlocal15(3,4,5),xlocal8(3,4,6),xlocal8r(3,1,6), 41 & disp(3),xload(2,*),delexternalwork 42! 43 include "gauss.f" 44 include "xlocal.f" 45! 46 data ifaceq /4,3,2,1,11,10,9,12, 47 & 5,6,7,8,13,14,15,16, 48 & 1,2,6,5,9,18,13,17, 49 & 2,3,7,6,10,19,14,18, 50 & 3,4,8,7,11,20,15,19, 51 & 4,1,5,8,12,17,16,20/ 52 data ifacet /1,3,2,7,6,5, 53 & 1,2,4,5,9,8, 54 & 2,3,4,6,10,9, 55 & 1,4,3,8,10,7/ 56 data ifacew /1,3,2,9,8,7,0,0, 57 & 4,5,6,10,11,12,0,0, 58 & 1,2,5,4,7,14,10,13, 59 & 2,3,6,5,8,15,11,14, 60 & 4,6,3,1,12,15,9,13/ 61 data iflag /3/ 62! 63! determining the set with the external faces 64! 65 ipos=index(externalfaces,' ') 66 externalfaces(ipos:ipos)='T' 67 do iset=1,nset 68 if(set(iset)(1:ipos).eq.externalfaces(1:ipos)) exit 69 enddo 70 externalfaces(ipos:ipos)=' ' 71! 72 delexternalwork=0.d0 73! 74 do jj=istartset(iset),iendset(iset) 75! 76 jface=ialset(jj) 77! 78 nelem=int(jface/10.d0) 79c** start change 20191219 80 indexe=ipkon(nelem) 81 if(indexe.lt.0) cycle 82c** end chenge 20191219 83 ig=jface-10*nelem 84! 85! determining the distributed load on the face 86! 87 ipres=0 88 call nident2(nelemload,nelem,nload,id) 89 do 90 if((id.eq.0).or.(nelemload(1,id).ne.nelem)) exit 91 if(sideload(id)(1:1).ne.'P') then 92 id=id-1 93 cycle 94 endif 95 igpres=ichar(sideload(id)(2:2))-48 96 if(ig.ne.igpres) then 97 id=id-1 98 cycle 99 else 100 ipres=1 101 pres=xload(1,id) 102 exit 103 endif 104 enddo 105c if(ipres.eq.0) then 106c write(*,*) '*ERROR in calcexternalwork: external face', 107c & ig,' of element ',nelem,' has no external load' 108c call exit(201) 109c endif 110! 111c indexe=ipkon(nelem) 112 lakonl=lakon(nelem) 113! 114! determining the number of nodes in the face 115! 116 if(lakonl(4:4).eq.'2') then 117 nope=20 118 nopes=8 119 elseif(lakonl(4:4).eq.'8') then 120 nope=8 121 nopes=4 122 elseif(lakonl(4:5).eq.'10') then 123 nope=10 124 nopes=6 125 elseif(lakonl(4:4).eq.'4') then 126 nope=4 127 nopes=3 128 elseif(lakonl(4:5).eq.'15') then 129 nope=15 130 elseif(lakonl(4:4).eq.'6') then 131 nope=6 132 endif 133! 134! determining the number of integration points in the face 135! 136 if(lakonl(4:5).eq.'8R') then 137 mint2d=1 138 elseif((lakonl(4:4).eq.'8').or.(lakonl(4:6).eq.'20R')) 139 & then 140 if((lakonl(7:7).eq.'A').or.(lakonl(7:7).eq.'S').or. 141 & (lakonl(7:7).eq.'E')) then 142 mint2d=2 143 else 144 mint2d=4 145 endif 146 elseif(lakonl(4:4).eq.'2') then 147 mint2d=9 148 elseif(lakonl(4:5).eq.'10') then 149 mint2d=3 150 elseif(lakonl(4:4).eq.'4') then 151 mint2d=1 152 endif 153! 154! local topology 155! 156 do i=1,nope 157 konl(i)=kon(indexe+i) 158 enddo 159! 160! treatment of wedge faces 161! 162 if(lakonl(4:4).eq.'6') then 163 mint2d=1 164 if(ig.le.2) then 165 nopes=3 166 else 167 nopes=4 168 endif 169 endif 170 if(lakonl(4:5).eq.'15') then 171 if(ig.le.2) then 172 mint2d=3 173 nopes=6 174 else 175 mint2d=4 176 nopes=8 177 endif 178 endif 179 180! 181! coordinates and differential displacements of the nodes belonging 182! to the face face 183! 184 if((nope.eq.20).or.(nope.eq.8)) then 185 do i=1,nopes 186 do j=1,3 187 voldl2(j,i)=vold(j,konl(ifaceq(i,ig))) 188 xl2(j,i)=co(j,konl(ifaceq(i,ig))) 189 & +voldl2(j,i) 190 voldl2(j,i)=voldl2(j,i)-voldprev(j,konl(ifaceq(i,ig))) 191 enddo 192 enddo 193 elseif((nope.eq.10).or.(nope.eq.4)) then 194 do i=1,nopes 195 do j=1,3 196 voldl2(j,i)=vold(j,konl(ifacet(i,ig))) 197 xl2(j,i)=co(j,konl(ifacet(i,ig))) 198 & +voldl2(j,i) 199 voldl2(j,i)=voldl2(j,i)-voldprev(j,konl(ifacet(i,ig))) 200 enddo 201 enddo 202 else 203 do i=1,nopes 204 do j=1,3 205 voldl2(j,i)=vold(j,konl(ifacew(i,ig))) 206 xl2(j,i)=co(j,konl(ifacew(i,ig))) 207 & +voldl2(j,i) 208 voldl2(j,i)=voldl2(j,i)-voldprev(j,konl(ifacew(i,ig))) 209 enddo 210 enddo 211 endif 212! 213 do i=1,mint2d 214! 215! local coordinates of the surface integration 216! point within the surface local coordinate system 217! 218 if((lakonl(4:5).eq.'8R').or. 219 & ((lakonl(4:4).eq.'6').and.(nopes.eq.4))) then 220 xi=gauss2d1(1,i) 221 et=gauss2d1(2,i) 222 weight=weight2d1(i) 223 elseif((lakonl(4:4).eq.'8').or. 224 & (lakonl(4:6).eq.'20R').or. 225 & ((lakonl(4:5).eq.'15').and.(nopes.eq.8))) then 226 xi=gauss2d2(1,i) 227 et=gauss2d2(2,i) 228 weight=weight2d2(i) 229 elseif(lakonl(4:4).eq.'2') then 230 xi=gauss2d3(1,i) 231 et=gauss2d3(2,i) 232 weight=weight2d3(i) 233 elseif((lakonl(4:5).eq.'10').or. 234 & ((lakonl(4:5).eq.'15').and.(nopes.eq.6))) then 235 xi=gauss2d5(1,i) 236 et=gauss2d5(2,i) 237 weight=weight2d5(i) 238 elseif((lakonl(4:4).eq.'4').or. 239 & ((lakonl(4:4).eq.'6').and.(nopes.eq.3))) then 240 xi=gauss2d4(1,i) 241 et=gauss2d4(2,i) 242 weight=weight2d4(i) 243 endif 244! 245! local surface normal 246! 247 if(nopes.eq.8) then 248 call shape8q(xi,et,xl2,xsj2,xs2,shp2,iflag) 249 elseif(nopes.eq.4) then 250 call shape4q(xi,et,xl2,xsj2,xs2,shp2,iflag) 251 elseif(nopes.eq.6) then 252 call shape6tri(xi,et,xl2,xsj2,xs2,shp2,iflag) 253 else 254 call shape3tri(xi,et,xl2,xsj2,xs2,shp2,iflag) 255 endif 256! 257! differential displacements at the integration point 258! 259 do j1=1,3 260 disp(j1)=0.d0 261 do i1=1,nopes 262 disp(j1)=disp(j1)+shp2(4,i1)*voldl2(j1,i1) 263 enddo 264 enddo 265! 266! total external work (pressure is a negative load) 267! 268 delexternalwork=delexternalwork-pres*weight* 269 & (disp(1)*xsj2(1)+disp(2)*xsj2(2)+disp(3)*xsj2(3)) 270 enddo 271 enddo 272! 273 return 274 end 275