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 printoutintfluidfem(prlab,ipkon,lakon,stx, 20 & mi,ii,nelem,qfx,orab,ielorien,norien,co,kon, 21 & ielmat,thicke,ielprop,prop,nelel,ithermal,orname,nelnew) 22! 23! stores integration point results for element "nelem" in the .dat file 24! 25! nelem is the element number from the input deck 26! nelel is a local number, created e.g. by renumbering; thus 27! far only applicable for FVM-CFD; for all other applications 28! nelem=nelel 29! 30 implicit none 31! 32 character*6 prlab(*) 33 character*8 lakon(*) 34 character*80 orname(*) 35! 36 integer ipkon(*),mi(*),nelem,l,ii,mint3d,j,k,nope, 37 & ielorien(mi(3),*),norien,kon(*),konl,indexe,m,iorien,iflag, 38 & ielmat(mi(3),*),nopes,mint2d,kk,ki,kl,nlayer,ilayer, 39 & null,ielprop(*),nelel,ithermal(*),nelef,nelnew(*) 40! 41 real*8 stx(6,mi(1),*), 42 & qfx(3,mi(1),*),xi,et,ze,xl(3,20),xsj,shp(4,20), 43 & coords(3,27),weight,orab(7,*),co(3,*),a(3,3),b(3,3),c(3,3), 44 & qfxl(3),thicke(mi(3),*),xsj2(3),shp2(7,8),xl2(3,8),xs2(3,7), 45 & thickness,tlayer(4),dlayer(4),xlayer(mi(3),4), 46 & prop(*) 47! 48 include "gauss.f" 49! 50 data iflag /1/ 51 data null /0/ 52! 53 write(*,*) 'printoutintfluidfem ',nelel 54 write(*,*) 'printoutintfluidfem new ',nelnew(nelel) 55 nelef=nelnew(nelel) 56! 57 if(ipkon(nelef).lt.0) then 58 return 59 else 60 indexe=ipkon(nelef) 61 endif 62! 63! check whether transformation is necessary (if orientation 64! is applied and output in local system is requested) 65! 66 if(lakon(nelef)(7:8).ne.'LC') then 67 if((norien.eq.0).or.(prlab(ii)(6:6).eq.'G')) then 68 iorien=0 69 else 70 iorien=max(0,ielorien(1,nelef)) 71 endif 72 elseif(lakon(nelef)(4:5).eq.'20') then 73! 74! composite materials 75! 76 if((norien.eq.0).or.(prlab(ii)(6:6).eq.'G')) then 77 iorien=0 78 else 79! 80! check whether at least one layer has a transformation 81! 82 iorien=0 83 do k=1,mi(3) 84 if(ielorien(k,nelef).ne.0) then 85 iorien=max(0,ielorien(k,nelef)) 86 exit 87 endif 88 enddo 89 endif 90! 91! determining the number of layers 92! 93 mint2d=4 94 nopes=8 95! 96 nlayer=0 97 do k=1,mi(3) 98 if(ielmat(k,nelef).ne.0) then 99 nlayer=nlayer+1 100 endif 101 enddo 102! 103! determining the layer thickness and global thickness 104! at the shell integration points 105! 106 iflag=1 107 do kk=1,mint2d 108 xi=gauss3d2(1,kk) 109 et=gauss3d2(2,kk) 110 call shape8q(xi,et,xl2,xsj2,xs2,shp2,iflag) 111 tlayer(kk)=0.d0 112 do k=1,nlayer 113 thickness=0.d0 114 do j=1,nopes 115 thickness=thickness+thicke(k,indexe+j)*shp2(4,j) 116 enddo 117 tlayer(kk)=tlayer(kk)+thickness 118 xlayer(k,kk)=thickness 119 enddo 120 enddo 121 iflag=3 122! 123 ilayer=0 124 do k=1,4 125 dlayer(k)=0.d0 126 enddo 127! 128! 129 elseif(lakon(nelef)(4:5).eq.'15') then 130! 131! composite materials 132! 133 if((norien.eq.0).or.(prlab(ii)(6:6).eq.'G')) then 134 iorien=0 135 else 136! 137! check whether at least one layer has a transformation 138! 139 iorien=0 140 do k=1,mi(3) 141 if(ielorien(k,nelef).ne.0) then 142 iorien=max(0,ielorien(k,nelef)) 143 exit 144 endif 145 enddo 146 endif 147! 148! determining the number of layers 149! 150 mint2d=3 151 nopes=6 152! 153 nlayer=0 154 do k=1,mi(3) 155 if(ielmat(k,nelef).ne.0) then 156 nlayer=nlayer+1 157 endif 158 enddo 159! 160! determining the layer thickness and global thickness 161! at the shell integration points 162! 163 iflag=1 164 do kk=1,mint2d 165 xi=gauss3d10(1,kk) 166 et=gauss3d10(2,kk) 167 call shape6tri(xi,et,xl2,xsj2,xs2,shp2,iflag) 168 tlayer(kk)=0.d0 169 do k=1,nlayer 170 thickness=0.d0 171 do j=1,nopes 172 thickness=thickness+thicke(k,indexe+j)*shp2(4,j) 173 enddo 174 tlayer(kk)=tlayer(kk)+thickness 175 xlayer(k,kk)=thickness 176 enddo 177 enddo 178 iflag=3 179! 180 ilayer=0 181 do k=1,3 182 dlayer(k)=0.d0 183 enddo 184! 185 endif 186! 187! number of integration points 188! 189 if(lakon(nelef)(4:5).eq.'8R') then 190 mint3d=1 191 elseif(lakon(nelef)(4:7).eq.'20RB') then 192 if((lakon(nelef)(8:8).eq.'R').or. 193 & (lakon(nelef)(8:8).eq.'C')) then 194 mint3d=50 195 else 196 call beamintscheme(lakon(nelef),mint3d,ielprop(nelel),prop, 197 & null,xi,et,ze,weight) 198 endif 199 elseif((lakon(nelef)(4:4).eq.'8').or. 200 & (lakon(nelef)(4:6).eq.'20R')) then 201 if(lakon(nelef)(7:8).eq.'LC') then 202 mint3d=8*nlayer 203 else 204 mint3d=8 205 endif 206 elseif(lakon(nelef)(4:4).eq.'2') then 207 mint3d=27 208 elseif(lakon(nelef)(4:5).eq.'10') then 209 mint3d=4 210 elseif(lakon(nelef)(4:4).eq.'4') then 211 mint3d=1 212 elseif(lakon(nelef)(4:5).eq.'15') then 213 if(lakon(nelef)(7:8).eq.'LC') then 214 mint3d=6*nlayer 215 else 216 mint3d=9 217 endif 218 elseif(lakon(nelef)(4:4).eq.'6') then 219 mint3d=2 220 elseif(lakon(nelef)(1:1).eq.'U') then 221 mint3d=ichar(lakon(nelef)(6:6)) 222 else 223 return 224 endif 225! 226! calculation of the integration point coordinates for 227! output in the local system (if needed) 228! 229 if((iorien.ne.0).or.(prlab(ii)(1:4).eq.'COOR')) then 230 if(lakon(nelef)(4:4).eq.'2') then 231 nope=20 232 elseif(lakon(nelef)(4:4).eq.'8') then 233 nope=8 234 elseif(lakon(nelef)(4:5).eq.'10') then 235 nope=10 236 elseif(lakon(nelef)(4:4).eq.'4') then 237 nope=4 238 elseif(lakon(nelef)(4:5).eq.'15') then 239 nope=15 240 elseif(lakon(nelef)(4:4).eq.'6') then 241 nope=6 242 endif 243! 244 do j=1,nope 245 konl=kon(indexe+j) 246 do k=1,3 247 xl(k,j)=co(k,konl) 248 enddo 249 enddo 250! 251 do j=1,mint3d 252 if((lakon(nelef)(4:5).eq.'8R').or. 253 & (lakon(nelef)(1:4).eq.'F3D8')) then 254 xi=gauss3d1(1,j) 255 et=gauss3d1(2,j) 256 ze=gauss3d1(3,j) 257 weight=weight3d1(j) 258 elseif(lakon(nelef)(4:8).eq.'20RB') then 259 if((lakon(nelef)(8:8).eq.'B').or. 260 & (lakon(nelef)(8:8).eq.'C')) then 261 xi=gauss3d13(1,j) 262 et=gauss3d13(2,j) 263 ze=gauss3d13(3,j) 264 weight=weight3d13(j) 265 else 266 call beamintscheme(lakon(nelef),mint3d,ielprop(nelel), 267 & prop,j,xi,et,ze,weight) 268 endif 269 elseif((lakon(nelef)(4:4).eq.'8').or. 270 & (lakon(nelef)(4:6).eq.'20R')) 271 & then 272 if(lakon(nelef)(7:8).ne.'LC') then 273 xi=gauss3d2(1,j) 274 et=gauss3d2(2,j) 275 ze=gauss3d2(3,j) 276 weight=weight3d2(j) 277 else 278 kl=mod(j,8) 279 if(kl.eq.0) kl=8 280! 281 xi=gauss3d2(1,kl) 282 et=gauss3d2(2,kl) 283 ze=gauss3d2(3,kl) 284 weight=weight3d2(kl) 285! 286 ki=mod(j,4) 287 if(ki.eq.0) ki=4 288! 289 if(kl.eq.1) then 290 ilayer=ilayer+1 291 if(ilayer.gt.1) then 292 do k=1,4 293 dlayer(k)=dlayer(k)+xlayer(ilayer-1,k) 294 enddo 295 endif 296 endif 297 ze=2.d0*(dlayer(ki)+(ze+1.d0)/2.d0*xlayer(ilayer,ki))/ 298 & tlayer(ki)-1.d0 299 weight=weight*xlayer(ilayer,ki)/tlayer(ki) 300 endif 301 elseif(lakon(nelef)(4:4).eq.'2') then 302 xi=gauss3d3(1,j) 303 et=gauss3d3(2,j) 304 ze=gauss3d3(3,j) 305 weight=weight3d3(j) 306 elseif(lakon(nelef)(4:5).eq.'10') then 307 xi=gauss3d5(1,j) 308 et=gauss3d5(2,j) 309 ze=gauss3d5(3,j) 310 weight=weight3d5(j) 311 elseif(lakon(nelef)(4:4).eq.'4') then 312 xi=gauss3d4(1,j) 313 et=gauss3d4(2,j) 314 ze=gauss3d4(3,j) 315 weight=weight3d4(j) 316 elseif(lakon(nelef)(4:5).eq.'15') then 317 if(lakon(nelef)(7:8).ne.'LC') then 318 xi=gauss3d8(1,j) 319 et=gauss3d8(2,j) 320 ze=gauss3d8(3,j) 321 weight=weight3d8(j) 322 else 323 kl=mod(j,6) 324 if(kl.eq.0) kl=6 325! 326 xi=gauss3d10(1,kl) 327 et=gauss3d10(2,kl) 328 ze=gauss3d10(3,kl) 329 weight=weight3d10(kl) 330! 331 ki=mod(j,3) 332 if(ki.eq.0) ki=3 333! 334 if(kl.eq.1) then 335 ilayer=ilayer+1 336 if(ilayer.gt.1) then 337 do k=1,3 338 dlayer(k)=dlayer(k)+xlayer(ilayer-1,k) 339 enddo 340 endif 341 endif 342 ze=2.d0*(dlayer(ki)+(ze+1.d0)/2.d0*xlayer(ilayer,ki))/ 343 & tlayer(ki)-1.d0 344 weight=weight*xlayer(ilayer,ki)/tlayer(ki) 345 endif 346 elseif(lakon(nelef)(1:4).eq.'C3D6') then 347 xi=gauss3d7(1,j) 348 et=gauss3d7(2,j) 349 ze=gauss3d7(3,j) 350 weight=weight3d7(j) 351 elseif(lakon(nelef)(1:4).eq.'F3D6') then 352 xi=gauss3d14(1,j) 353 et=gauss3d14(2,j) 354 ze=gauss3d14(3,j) 355 weight=weight3d14(j) 356 endif 357! 358 if(nope.eq.20) then 359 call shape20h(xi,et,ze,xl,xsj,shp,iflag) 360 elseif(nope.eq.8) then 361 call shape8h(xi,et,ze,xl,xsj,shp,iflag) 362 elseif(nope.eq.10) then 363 call shape10tet(xi,et,ze,xl,xsj,shp,iflag) 364 elseif(nope.eq.4) then 365 call shape4tet(xi,et,ze,xl,xsj,shp,iflag) 366 elseif(nope.eq.15) then 367 call shape15w(xi,et,ze,xl,xsj,shp,iflag) 368 else 369 call shape6w(xi,et,ze,xl,xsj,shp,iflag) 370 endif 371! 372 do k=1,3 373 coords(k,j)=0.d0 374 do l=1,nope 375 coords(k,j)=coords(k,j)+xl(k,l)*shp(4,l) 376 enddo 377 enddo 378 enddo 379 endif 380! 381 if((prlab(ii)(1:4).eq.'S ').or.(prlab(ii)(1:4).eq.'SVF ')) then 382 do j=1,mint3d 383! 384! composite materials 385! 386 if(lakon(nelef)(7:8).eq.'LC') then 387 if((norien.eq.0).or.(prlab(ii)(6:6).eq.'G')) then 388 iorien=0 389 elseif(lakon(nelef)(4:5).eq.'20') then 390 kl=mod(j,8) 391 if(kl.eq.0) kl=8 392 ilayer=(j-kl)/8+1 393 iorien=max(0,ielorien(ilayer,nelef)) 394 elseif(lakon(nelef)(4:5).eq.'15') then 395 kl=mod(j,6) 396 if(kl.eq.0) kl=6 397 ilayer=(j-kl)/6+1 398 iorien=max(0,ielorien(ilayer,nelef)) 399 endif 400 endif 401! 402 if(iorien.eq.0) then 403 write(5,'(i10,1x,i3,1p,6(1x,e13.6))') nelem,j, 404 & (stx(k,j,nelef),k=1,6) 405 else 406 call transformatrix(orab(1,iorien),coords(1,j),a) 407 b(1,1)=stx(1,j,nelef) 408 b(2,2)=stx(2,j,nelef) 409 b(3,3)=stx(3,j,nelef) 410 b(1,2)=stx(4,j,nelef) 411 b(1,3)=stx(5,j,nelef) 412 b(2,3)=stx(6,j,nelef) 413 b(2,1)=b(1,2) 414 b(3,1)=b(1,3) 415 b(3,2)=b(2,3) 416 do k=1,3 417 do l=1,3 418 c(k,l)=0.d0 419 do m=1,3 420 c(k,l)=c(k,l)+b(k,m)*a(m,l) 421 enddo 422 enddo 423 enddo 424 do k=1,3 425 do l=k,3 426 b(k,l)=0.d0 427 do m=1,3 428 b(k,l)=b(k,l)+a(m,k)*c(m,l) 429 enddo 430 enddo 431 enddo 432 write(5,'(i10,1x,i3,1p,6(1x,e13.6),1x,a20)') nelem,j, 433 & b(1,1),b(2,2),b(3,3),b(1,2),b(1,3),b(2,3), 434 & orname(iorien)(1:20) 435 endif 436 enddo 437 elseif(((prlab(ii)(1:4).eq.'HFL ').or.(prlab(ii)(1:4).eq.'HFLF')) 438 & .and.(ithermal(1).gt.1)) then 439 do j=1,mint3d 440! 441! composite materials 442! 443 if(lakon(nelef)(7:8).eq.'LC') then 444 if((norien.eq.0).or.(prlab(ii)(6:6).eq.'G')) then 445 iorien=0 446 elseif(lakon(nelef)(4:5).eq.'20') then 447 kl=mod(j,8) 448 if(kl.eq.0) kl=8 449 ilayer=(j-kl)/8+1 450 iorien=max(0,ielorien(ilayer,nelef)) 451 elseif(lakon(nelef)(4:5).eq.'15') then 452 kl=mod(j,6) 453 if(kl.eq.0) kl=6 454 ilayer=(j-kl)/6+1 455 iorien=max(0,ielorien(ilayer,nelef)) 456 endif 457 endif 458! 459 if(iorien.eq.0) then 460 write(5,'(i10,1x,i3,1p,3(1x,e13.6))') nelem,j, 461 & (qfx(k,j,nelef),k=1,3) 462 else 463 do k=1,3 464 qfxl(k)=qfx(k,j,nelef) 465 enddo 466 call transformatrix(orab(1,iorien),coords(1,j),a) 467 write(5,'(i10,1x,i3,1p,3(1x,e13.6),1x,a20)') nelem,j, 468 & qfxl(1)*a(1,1)+qfxl(2)*a(2,1)+qfxl(3)*a(3,1), 469 & qfxl(1)*a(1,2)+qfxl(2)*a(2,2)+qfxl(3)*a(3,2), 470 & qfxl(1)*a(1,3)+qfxl(2)*a(2,3)+qfxl(3)*a(3,3), 471 & orname(iorien)(1:20) 472 endif 473 enddo 474 elseif(prlab(ii)(1:4).eq.'COOR') then 475 do j=1,mint3d 476 write(5,'(i10,1x,i3,1p,3(1x,e13.6))') nelem,j, 477 & (coords(k,j),k=1,3) 478 enddo 479 endif 480! 481 return 482 end 483