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 printoutfacefem(co,ntmat_,vold,shcon, 20 & nshcon,compressible,istartset,iendset,ipkon, 21 & lakon,kon,ialset,prset,ttime,nset,set,nprint,prlab,ielmat,mi, 22 & nelnew) 23! 24! calculation and printout of the lift and drag forces 25! 26 implicit none 27! 28 integer compressible 29! 30 character*8 lakonl,lakon(*) 31 character*6 prlab(*) 32 character*80 faset 33 character*81 set(*),prset(*) 34! 35 integer konl(8),ifaceq(8,6),nelem,ii,nprint,i,j,i1,i2,j1,id, 36 & k1,jj,ig,nshcon(*),ntmat_,nope,nopes,nelef,nelnew(*), 37 & imat,mint2d,ifacet(6,4),ifacew(8,5),iflag,indexe,jface, 38 & istartset(*),iendset(*),ipkon(*),kon(*),iset,ialset(*),nset, 39 & ipos,mi(*),ielmat(mi(3),*) 40! 41 real*8 co(3,*),xl(3,8),shp(4,8),xs2(3,7),dvi,f(3), 42 & vkl(3,3),t(3,3),div,shcon(0:3,ntmat_,*), 43 & voldl(0:mi(2),8),xl2(3,8),xsj2(3), 44 & shp2(7,8),vold(0:mi(2),*),xi,et,xsj,temp,xi3d,et3d,ze3d, 45 & weight,xlocal20(3,9,6),xlocal4(3,1,4),xlocal10(3,3,4), 46 & xlocal6(3,1,5),xlocal15(3,4,5),xlocal8(3,4,6), 47 & xlocal8r(3,1,6),ttime,pres,tf(3),tn,tt,dd,coords(3) 48! 49 include "gauss.f" 50 include "xlocal.f" 51! 52 data ifaceq /4,3,2,1,11,10,9,12, 53 & 5,6,7,8,13,14,15,16, 54 & 1,2,6,5,9,18,13,17, 55 & 2,3,7,6,10,19,14,18, 56 & 3,4,8,7,11,20,15,19, 57 & 4,1,5,8,12,17,16,20/ 58 data ifacet /1,3,2,7,6,5, 59 & 1,2,4,5,9,8, 60 & 2,3,4,6,10,9, 61 & 1,4,3,8,10,7/ 62 data ifacew /1,3,2,9,8,7,0,0, 63 & 4,5,6,10,11,12,0,0, 64 & 1,2,5,4,7,14,10,13, 65 & 2,3,6,5,8,15,11,14, 66 & 4,6,3,1,12,15,9,13/ 67 data iflag /3/ 68! 69! initialisierung forces 70! 71 do i=1,3 72 f(i)=0.d0 73 enddo 74! 75 do ii=1,nprint 76! 77! total drag 78! 79 if(prlab(ii)(1:4).eq.'DRAG') then 80! 81 ipos=index(prset(ii),' ') 82 faset=' ' 83 faset(1:ipos-1)=prset(ii)(1:ipos-1) 84! 85! printing the header 86! 87 write(5,*) 88 write(5,120) faset(1:ipos-2),ttime 89 120 format(' surface stress vector (tx,ty,tz), normal stress, sh 90 &ear stress and coordinates for set ',A,' and time ',e14.7) 91 write(5,*) 92! 93! printing the data 94! 95c do iset=1,nset 96c if(set(iset).eq.prset(ii)) exit 97c enddo 98 call cident81(set,prset(ii),nset,id) 99 iset=nset+1 100 if(id.gt.0) then 101 if(prset(ii).eq.set(id)) then 102 iset=id 103 endif 104 endif 105! 106 do jj=istartset(iset),iendset(iset) 107! 108 jface=ialset(jj) 109! 110 nelem=int(jface/10.d0) 111 ig=jface-10*nelem 112! 113 nelef=nelnew(nelem) 114! 115 lakonl=lakon(nelef) 116 indexe=ipkon(nelef) 117 imat=ielmat(1,nelef) 118! 119 if(lakonl(4:4).eq.'8') then 120 nope=8 121 nopes=4 122 elseif(lakonl(4:4).eq.'4') then 123 nope=4 124 nopes=3 125 elseif(lakonl(4:4).eq.'6') then 126 nope=6 127 endif 128! 129 if(lakonl(4:5).eq.'8R') then 130 mint2d=1 131 elseif(lakonl(4:4).eq.'8') then 132 mint2d=4 133 elseif(lakonl(4:4).eq.'4') then 134 mint2d=1 135 endif 136! 137! local topology 138! 139 do i=1,nope 140 konl(i)=kon(indexe+i) 141 enddo 142! 143! computation of the coordinates of the local nodes 144! 145 do i=1,nope 146 do j=1,3 147 xl(j,i)=co(j,konl(i)) 148 enddo 149 enddo 150! 151! temperature, velocity and auxiliary variables 152! (rho*energy density, rho*velocity and rho) 153! 154 do i1=1,nope 155 do i2=0,4 156 voldl(i2,i1)=vold(i2,konl(i1)) 157 enddo 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! 171 if(nope.eq.8) then 172 do i=1,nopes 173 do j=1,3 174 xl2(j,i)=co(j,konl(ifaceq(i,ig))) 175 enddo 176 enddo 177 elseif(nope.eq.4) then 178 do i=1,nopes 179 do j=1,3 180 xl2(j,i)=co(j,konl(ifacet(i,ig))) 181 enddo 182 enddo 183 else 184 do i=1,nopes 185 do j=1,3 186 xl2(j,i)=co(j,konl(ifacew(i,ig))) 187 enddo 188 enddo 189 endif 190! 191 do i=1,mint2d 192! 193! local coordinates of the surface integration 194! point within the surface local coordinate system 195! 196 if((lakonl(4:5).eq.'8R').or. 197 & ((lakonl(4:4).eq.'6').and.(nopes.eq.4))) then 198 xi=gauss2d1(1,i) 199 et=gauss2d1(2,i) 200 weight=weight2d1(i) 201 elseif(lakonl(4:4).eq.'8') then 202 xi=gauss2d2(1,i) 203 et=gauss2d2(2,i) 204 weight=weight2d2(i) 205 elseif((lakonl(4:4).eq.'4').or. 206 & ((lakonl(4:4).eq.'6').and.(nopes.eq.3))) then 207 xi=gauss2d4(1,i) 208 et=gauss2d4(2,i) 209 weight=weight2d4(i) 210 endif 211! 212! local surface normal 213! 214 if(nopes.eq.4) then 215 call shape4q(xi,et,xl2,xsj2,xs2,shp2,iflag) 216 else 217 call shape3tri(xi,et,xl2,xsj2,xs2,shp2,iflag) 218 endif 219! 220! global coordinates of the integration point 221! 222 do j1=1,3 223 coords(j1)=0.d0 224 do i1=1,nopes 225 coords(j1)=coords(j1)+shp2(4,i1)*xl2(j1,i1) 226 enddo 227 enddo 228! 229! local coordinates of the surface integration 230! point within the element local coordinate system 231! 232 if(lakonl(4:5).eq.'8R') then 233 xi3d=xlocal8r(1,i,ig) 234 et3d=xlocal8r(2,i,ig) 235 ze3d=xlocal8r(3,i,ig) 236 call shape8h(xi3d,et3d,ze3d,xl,xsj,shp,iflag) 237 elseif(lakonl(4:4).eq.'8') then 238 xi3d=xlocal8(1,i,ig) 239 et3d=xlocal8(2,i,ig) 240 ze3d=xlocal8(3,i,ig) 241 call shape8h(xi3d,et3d,ze3d,xl,xsj,shp,iflag) 242 elseif(lakonl(4:4).eq.'4') then 243 xi3d=xlocal4(1,i,ig) 244 et3d=xlocal4(2,i,ig) 245 ze3d=xlocal4(3,i,ig) 246 call shape4tet(xi3d,et3d,ze3d,xl,xsj,shp,iflag) 247 elseif(lakonl(4:4).eq.'6') then 248 xi3d=xlocal6(1,i,ig) 249 et3d=xlocal6(2,i,ig) 250 ze3d=xlocal6(3,i,ig) 251 call shape6w(xi3d,et3d,ze3d,xl,xsj,shp,iflag) 252 endif 253! 254! calculating of 255! the temperature temp 256! the static pressure pres 257! the velocity gradient vkl 258! in the integration point 259! 260 temp=0.d0 261 pres=0.d0 262 do i1=1,3 263 do j1=1,3 264 vkl(i1,j1)=0.d0 265 enddo 266 enddo 267 do i1=1,nope 268 temp=temp+shp(4,i1)*voldl(0,i1) 269 pres=pres+shp(4,i1)*voldl(4,i1) 270 do j1=1,3 271 do k1=1,3 272 vkl(j1,k1)=vkl(j1,k1)+shp(k1,i1)*voldl(j1,i1) 273 enddo 274 enddo 275 enddo 276 if(compressible.eq.1) div=vkl(1,1)+vkl(2,2)+vkl(3,3) 277! 278! material data (density, dynamic viscosity, heat capacity and 279! conductivity) 280! 281 call materialdata_dvifem(imat,ntmat_,temp,shcon,nshcon, 282 & dvi) 283! 284! determining the stress 285! 286 do i1=1,3 287 do j1=1,3 288 t(i1,j1)=vkl(i1,j1)+vkl(j1,i1) 289 enddo 290 if(compressible.eq.1) 291 & t(i1,i1)=t(i1,i1)-2.d0*div/3.d0 292 enddo 293! 294 dd=dsqrt(xsj2(1)*xsj2(1)+xsj2(2)*xsj2(2)+ 295 & xsj2(3)*xsj2(3)) 296 do i1=1,3 297 tf(i1)=dvi*(t(i1,1)*xsj2(1)+t(i1,2)*xsj2(2)+ 298 & t(i1,3)*xsj2(3))-pres*xsj2(i1) 299 f(i1)=f(i1)+tf(i1)*weight 300 tf(i1)=tf(i1)/dd 301 enddo 302 tn=(tf(1)*xsj2(1)+tf(2)*xsj2(2)+tf(3)*xsj2(3))/dd 303 tt=dsqrt((tf(1)-tn*xsj2(1)/dd)**2+ 304 & (tf(2)-tn*xsj2(2)/dd)**2+ 305 & (tf(3)-tn*xsj2(3)/dd)**2) 306 write(5,'(i6,1x,i3,1x,i3,1p,8(1x,e11.4))') nelem,ig,i, 307 & (tf(i1),i1=1,3),tn,tt,(coords(i1),i1=1,3) 308! 309 enddo 310 enddo 311! 312 write(5,*) 313 write(5,121) faset(1:ipos-2),ttime 314 121 format(' total surface force (fx,fy,fz) for set ',A, 315 & ' and time ',e14.7) 316 write(5,*) 317 write(5,'(1p,3(1x,e11.4))') (f(j),j=1,3) 318! 319 endif 320 enddo 321! 322 return 323 end 324 325 326