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 printoutface(co,rhcon,nrhcon,ntmat_,vold,shcon,nshcon, 20 & cocon,ncocon,compressible,istartset,iendset,ipkonf,lakonf,konf, 21 & ialset,prset,ttime,nset,set,nprint,prlab,ielmat,mi, 22 & ithermal,nactdoh,icfd,time,stn) 23! 24! calculation and printout of the lift and drag forces 25! 26 implicit none 27! 28 integer compressible 29! 30 character*8 lakonl,lakonf(*) 31 character*6 prlab(*) 32 character*80 faset 33 character*81 set(*),prset(*) 34! 35 integer konl(20),ifaceq(8,6),nelem,ii,nprint,i,j,i1,i2,j1, 36 & ncocon(2,*),k1,jj,ig,nrhcon(*),nshcon(*),ntmat_,nope,nopes,imat, 37 & mint2d,ifacet(6,4),ifacew(8,5),iflag,indexe,jface,istartset(*), 38 & iendset(*),ipkonf(*),konf(*),iset,ialset(*),nset,ipos,id, 39 & mi(*),ielmat(mi(3),*),nactdoh(*),icfd,nelemcfd,ithermal(*) 40! 41 real*8 co(3,*),xl(3,20),shp(4,20),xs2(3,7),dvi,f(0:3),time, 42 & vkl(0:3,3),rhcon(0:1,ntmat_,*),t(3,3),div,shcon(0:3,ntmat_,*), 43 & voldl(0:mi(2),20),cocon(0:6,ntmat_,*),xl2(3,8),xsj2(3), 44 & shp2(7,8),vold(0:mi(2),*),xi,et,xsj,temp,xi3d,et3d,ze3d,weight, 45 & xlocal20(3,9,6),xlocal4(3,1,4),xlocal10(3,3,4),xlocal6(3,1,5), 46 & xlocal15(3,4,5),xlocal8(3,4,6),xlocal8r(3,1,6),ttime,pres, 47 & tf(0:3),tn,tt,dd,coords(3),cond,stn(6,*),xm(3),df(3),cg(3), 48 & area,xn(3),xnormforc,shearforc 49! 50! 51 include "gauss.f" 52 include "xlocal.f" 53! 54 data ifaceq /4,3,2,1,11,10,9,12, 55 & 5,6,7,8,13,14,15,16, 56 & 1,2,6,5,9,18,13,17, 57 & 2,3,7,6,10,19,14,18, 58 & 3,4,8,7,11,20,15,19, 59 & 4,1,5,8,12,17,16,20/ 60 data ifacet /1,3,2,7,6,5, 61 & 1,2,4,5,9,8, 62 & 2,3,4,6,10,9, 63 & 1,4,3,8,10,7/ 64 data ifacew /1,3,2,9,8,7,0,0, 65 & 4,5,6,10,11,12,0,0, 66 & 1,2,5,4,7,14,10,13, 67 & 2,3,6,5,8,15,11,14, 68 & 4,6,3,1,12,15,9,13/ 69 data iflag /3/ 70! 71! flag to check whether sof and/or som output was already 72! printed 73! 74 do ii=1,nprint 75! 76! check whether there are facial print requests 77! 78! DRAG: drag forces in cfd-calculations 79! FLUX: heat flux 80! SOF: forces and moments on a section 81! (= an internal or external surface) 82! 83! DRAG removed on 13 Dec 2020: DRAG only accessible for FEM-CBS 84! through printoutfacefem.f 85! 86 if((prlab(ii)(1:4).eq.'FLUX').or. 87 & (prlab(ii)(1:3).eq.'SOF')) 88c if((prlab(ii)(1:4).eq.'DRAG').or.(prlab(ii)(1:4).eq.'FLUX').or. 89c & (prlab(ii)(1:3).eq.'SOF')) 90 & then 91! 92 ipos=index(prset(ii),' ') 93 do i=1,80 94 faset(i:i)=' ' 95 enddo 96c faset=' ' 97 faset(1:ipos-1)=prset(ii)(1:ipos-1) 98! 99! printing the header 100! 101 write(5,*) 102 if(prlab(ii)(1:4).eq.'DRAG') then 103! 104! initialisierung forces 105! 106 do i=1,3 107 f(i)=0.d0 108 enddo 109! 110 write(5,120) faset(1:ipos-2),ttime+time 111 120 format( 112 & ' surface stress at the integration points for set ', 113 & A,' and time ',e14.7) 114 write(5,*) 115 write(5,124) 116 124 format(' el fa int tx ty 117 & tz normal stress shear stress and coordinates 118 &') 119 elseif(prlab(ii)(1:4).eq.'FLUX') then 120! 121! initialisierung of the flux 122! 123 f(0)=0.d0 124 write(5,121) faset(1:ipos-2),ttime+time 125 121 format(' heat flux at the integration points for set ', 126 & A,' and time ',e14.7) 127 write(5,*) 128 write(5,125) 129 125 format(' el fa int heat flux q and c 130 &oordinates') 131 else 132! 133! initialize total force and total moment 134! 135 do i=1,3 136 f(i)=0.d0 137 xm(i)=0.d0 138 cg(i)=0.d0 139 xn(i)=0.d0 140 enddo 141 area=0.d0 142 endif 143 write(5,*) 144! 145! printing the data 146! 147c do iset=1,nset 148c if(set(iset).eq.prset(ii)) exit 149c enddo 150 call cident81(set,prset(ii),nset,id) 151 iset=nset+1 152 if(id.gt.0) then 153 if(prset(ii).eq.set(id)) then 154 iset=id 155 endif 156 endif 157! 158 do jj=istartset(iset),iendset(iset) 159! 160 jface=ialset(jj) 161! 162 nelem=int(jface/10.d0) 163 ig=jface-10*nelem 164c write(*,*) 'printoutface elem ig ',nelem,ig 165! 166! for CFD calculations the elements were renumbered 167! 168 if(icfd.ne.0) then 169 nelemcfd=nactdoh(nelem) 170 indexe=ipkonf(nelemcfd) 171 lakonl=lakonf(nelemcfd) 172 imat=ielmat(1,nelemcfd) 173 else 174 indexe=ipkonf(nelem) 175 lakonl=lakonf(nelem) 176 imat=ielmat(1,nelem) 177 endif 178! 179 if(lakonl(4:4).eq.'2') then 180 nope=20 181 nopes=8 182 elseif(lakonl(4:4).eq.'8') then 183 nope=8 184 nopes=4 185 elseif(lakonl(4:5).eq.'10') then 186 nope=10 187 nopes=6 188 elseif(lakonl(4:4).eq.'4') then 189 nope=4 190 nopes=3 191 elseif(lakonl(4:5).eq.'15') then 192 nope=15 193 elseif(lakonl(4:4).eq.'6') then 194 nope=6 195 endif 196! 197 if(lakonl(4:5).eq.'8R') then 198 mint2d=1 199 elseif((lakonl(4:4).eq.'8').or.(lakonl(4:6).eq.'20R')) 200 & then 201 if((lakonl(7:7).eq.'A').or.(lakonl(7:7).eq.'S').or. 202 & (lakonl(7:7).eq.'E')) then 203 mint2d=2 204 else 205 mint2d=4 206 endif 207 elseif(lakonl(4:4).eq.'2') then 208 mint2d=9 209 elseif(lakonl(4:5).eq.'10') then 210 mint2d=3 211 elseif(lakonl(4:4).eq.'4') then 212 mint2d=1 213 endif 214! 215! local topology 216! 217 do i=1,nope 218 konl(i)=konf(indexe+i) 219 enddo 220! 221! computation of the coordinates of the local nodes 222! 223 do i=1,nope 224 do j=1,3 225 xl(j,i)=co(j,konl(i)) 226 enddo 227 enddo 228! 229! temperature, displacement (structures) or 230! temperature, velocities and pressure (fluids) 231! 232 do i1=1,nope 233 do i2=0,mi(2) 234 voldl(i2,i1)=vold(i2,konl(i1)) 235 enddo 236 enddo 237! 238! for structural calculations: adding the displacements to the 239! coordinates 240! 241 if(icfd.eq.0) then 242 do i=1,nope 243 do j=1,3 244 xl(j,i)=xl(j,i)+voldl(j,i) 245 enddo 246 enddo 247 endif 248! 249! treatment of wedge faces 250! 251 if(lakonl(4:4).eq.'6') then 252 mint2d=1 253 if(ig.le.2) then 254 nopes=3 255 else 256 nopes=4 257 endif 258 endif 259 if(lakonl(4:5).eq.'15') then 260 if(ig.le.2) then 261 mint2d=3 262 nopes=6 263 else 264 mint2d=4 265 nopes=8 266 endif 267 endif 268! 269 if(icfd.ne.0) then 270! 271! CFD: undeformed structure 272! 273 if((nope.eq.20).or.(nope.eq.8)) then 274 do i=1,nopes 275 do j=1,3 276 xl2(j,i)=co(j,konl(ifaceq(i,ig))) 277 enddo 278 enddo 279 elseif((nope.eq.10).or.(nope.eq.4)) then 280 do i=1,nopes 281 do j=1,3 282 xl2(j,i)=co(j,konl(ifacet(i,ig))) 283 enddo 284 enddo 285 else 286 do i=1,nopes 287 do j=1,3 288 xl2(j,i)=co(j,konl(ifacew(i,ig))) 289 enddo 290 enddo 291 endif 292 else 293! 294! no CFD: deformed structure 295! 296 if((nope.eq.20).or.(nope.eq.8)) then 297 do i=1,nopes 298 do j=1,3 299 xl2(j,i)=co(j,konl(ifaceq(i,ig))) 300 & +vold(j,konl(ifaceq(i,ig))) 301 enddo 302 enddo 303 elseif((nope.eq.10).or.(nope.eq.4)) then 304 do i=1,nopes 305 do j=1,3 306 xl2(j,i)=co(j,konl(ifacet(i,ig))) 307 & +vold(j,konl(ifacet(i,ig))) 308 enddo 309 enddo 310 else 311 do i=1,nopes 312 do j=1,3 313 xl2(j,i)=co(j,konl(ifacew(i,ig)))+ 314 & vold(j,konl(ifacew(i,ig))) 315 enddo 316 enddo 317 endif 318 endif 319! 320 do i=1,mint2d 321! 322! local coordinates of the surface integration 323! point within the surface local coordinate system 324! 325 if((lakonl(4:5).eq.'8R').or. 326 & ((lakonl(4:4).eq.'6').and.(nopes.eq.4))) then 327 xi=gauss2d1(1,i) 328 et=gauss2d1(2,i) 329 weight=weight2d1(i) 330 elseif((lakonl(4:4).eq.'8').or. 331 & (lakonl(4:6).eq.'20R').or. 332 & ((lakonl(4:5).eq.'15').and.(nopes.eq.8))) then 333 xi=gauss2d2(1,i) 334 et=gauss2d2(2,i) 335 weight=weight2d2(i) 336 elseif(lakonl(4:4).eq.'2') then 337 xi=gauss2d3(1,i) 338 et=gauss2d3(2,i) 339 weight=weight2d3(i) 340 elseif((lakonl(4:5).eq.'10').or. 341 & ((lakonl(4:5).eq.'15').and.(nopes.eq.6))) then 342 xi=gauss2d5(1,i) 343 et=gauss2d5(2,i) 344 weight=weight2d5(i) 345 elseif((lakonl(4:4).eq.'4').or. 346 & ((lakonl(4:4).eq.'6').and.(nopes.eq.3))) then 347 xi=gauss2d4(1,i) 348 et=gauss2d4(2,i) 349 weight=weight2d4(i) 350 endif 351! 352! local surface normal 353! 354 if(nopes.eq.8) then 355 call shape8q(xi,et,xl2,xsj2,xs2,shp2,iflag) 356 elseif(nopes.eq.4) then 357 call shape4q(xi,et,xl2,xsj2,xs2,shp2,iflag) 358 elseif(nopes.eq.6) then 359 call shape6tri(xi,et,xl2,xsj2,xs2,shp2,iflag) 360 else 361 call shape3tri(xi,et,xl2,xsj2,xs2,shp2,iflag) 362 endif 363! 364! global coordinates of the integration point 365! 366 do j1=1,3 367 coords(j1)=0.d0 368 do i1=1,nopes 369 coords(j1)=coords(j1)+shp2(4,i1)*xl2(j1,i1) 370 enddo 371 enddo 372! 373! local coordinates of the surface integration 374! point within the element local coordinate system 375! 376 if((prlab(ii)(1:4).eq.'DRAG').or. 377 & (prlab(ii)(1:4).eq.'FLUX')) then 378! 379! deformation gradient is only needed for 380! DRAG and FLUX applications 381! 382 if(lakonl(4:5).eq.'8R') then 383 xi3d=xlocal8r(1,i,ig) 384 et3d=xlocal8r(2,i,ig) 385 ze3d=xlocal8r(3,i,ig) 386 call shape8h(xi3d,et3d,ze3d,xl,xsj,shp,iflag) 387 elseif(lakonl(4:4).eq.'8') then 388 xi3d=xlocal8(1,i,ig) 389 et3d=xlocal8(2,i,ig) 390 ze3d=xlocal8(3,i,ig) 391 call shape8h(xi3d,et3d,ze3d,xl,xsj,shp,iflag) 392 elseif(lakonl(4:6).eq.'20R') then 393 xi3d=xlocal8(1,i,ig) 394 et3d=xlocal8(2,i,ig) 395 ze3d=xlocal8(3,i,ig) 396 call shape20h(xi3d,et3d,ze3d,xl,xsj,shp,iflag) 397 elseif(lakonl(4:4).eq.'2') then 398 xi3d=xlocal20(1,i,ig) 399 et3d=xlocal20(2,i,ig) 400 ze3d=xlocal20(3,i,ig) 401 call shape20h(xi3d,et3d,ze3d,xl,xsj,shp,iflag) 402 elseif(lakonl(4:5).eq.'10') then 403 xi3d=xlocal10(1,i,ig) 404 et3d=xlocal10(2,i,ig) 405 ze3d=xlocal10(3,i,ig) 406 call shape10tet(xi3d,et3d,ze3d,xl,xsj,shp,iflag) 407 elseif(lakonl(4:4).eq.'4') then 408 xi3d=xlocal4(1,i,ig) 409 et3d=xlocal4(2,i,ig) 410 ze3d=xlocal4(3,i,ig) 411 call shape4tet(xi3d,et3d,ze3d,xl,xsj,shp,iflag) 412 elseif(lakonl(4:5).eq.'15') then 413 xi3d=xlocal15(1,i,ig) 414 et3d=xlocal15(2,i,ig) 415 ze3d=xlocal15(3,i,ig) 416 call shape15w(xi3d,et3d,ze3d,xl,xsj,shp,iflag) 417 elseif(lakonl(4:4).eq.'6') then 418 xi3d=xlocal6(1,i,ig) 419 et3d=xlocal6(2,i,ig) 420 ze3d=xlocal6(3,i,ig) 421 call shape6w(xi3d,et3d,ze3d,xl,xsj,shp,iflag) 422 endif 423 endif 424! 425! calculating of 426! the temperature temp 427! the static pressure pres 428! the velocity gradient vkl 429! in the integration point 430! 431 if(prlab(ii)(1:4).eq.'DRAG') then 432 temp=0.d0 433 pres=0.d0 434 do i1=1,3 435 do j1=1,3 436 vkl(i1,j1)=0.d0 437 enddo 438 enddo 439 do i1=1,nope 440 temp=temp+shp(4,i1)*voldl(0,i1) 441 pres=pres+shp(4,i1)*voldl(4,i1) 442 do j1=1,3 443 do k1=1,3 444 vkl(j1,k1)=vkl(j1,k1)+ 445 & shp(k1,i1)*voldl(j1,i1) 446 enddo 447 enddo 448 enddo 449 if(compressible.eq.1) 450 & div=vkl(1,1)+vkl(2,2)+vkl(3,3) 451! 452! material data (density, dynamic viscosity, heat capacity and 453! conductivity) 454! 455 call materialdata_dvi(shcon,nshcon,imat,dvi,temp, 456 & ntmat_,ithermal) 457! 458! determining the stress 459! 460 do i1=1,3 461 do j1=1,3 462 t(i1,j1)=vkl(i1,j1)+vkl(j1,i1) 463 enddo 464 if(compressible.eq.1) 465 & t(i1,i1)=t(i1,i1)-2.d0*div/3.d0 466 enddo 467! 468 dd=dsqrt(xsj2(1)*xsj2(1)+xsj2(2)*xsj2(2)+ 469 & xsj2(3)*xsj2(3)) 470 do i1=1,3 471 tf(i1)=dvi*(t(i1,1)*xsj2(1)+t(i1,2)*xsj2(2)+ 472 & t(i1,3)*xsj2(3))-pres*xsj2(i1) 473 f(i1)=f(i1)+tf(i1)*weight 474 tf(i1)=tf(i1)/dd 475 enddo 476 tn=(tf(1)*xsj2(1)+tf(2)*xsj2(2)+tf(3)*xsj2(3))/dd 477 tt=dsqrt((tf(1)-tn*xsj2(1)/dd)**2+ 478 & (tf(2)-tn*xsj2(2)/dd)**2+ 479 & (tf(3)-tn*xsj2(3)/dd)**2) 480 write(5,'(i10,1x,i3,1x,i3,1p,8(1x,e13.6))')nelem, 481 & ig,i,(tf(i1),i1=1,3),tn,tt,(coords(i1),i1=1,3) 482! 483 elseif(prlab(ii)(1:4).eq.'FLUX') then 484 temp=0.d0 485 do j1=1,3 486 vkl(0,j1)=0.d0 487 enddo 488 do i1=1,nope 489 temp=temp+shp(4,i1)*voldl(0,i1) 490 do k1=1,3 491 vkl(0,k1)=vkl(0,k1)+shp(k1,i1)*voldl(0,i1) 492 enddo 493 enddo 494! 495! material data (conductivity) 496! 497 call materialdata_cond(imat,ntmat_,temp,cocon, 498 & ncocon,cond) 499! 500! determining the stress 501! 502 dd=dsqrt(xsj2(1)*xsj2(1)+xsj2(2)*xsj2(2)+ 503 & xsj2(3)*xsj2(3)) 504! 505 tf(0)=-cond*(vkl(0,1)*xsj2(1)+ 506 & vkl(0,2)*xsj2(2)+ 507 & vkl(0,3)*xsj2(3)) 508 f(0)=f(0)+tf(0)*weight 509 tf(0)=tf(0)/dd 510 write(5,'(i10,1x,i3,1x,i3,1p,4(1x,e13.6))')nelem, 511 & ig,i,tf(0),(coords(i1),i1=1,3) 512! 513 else 514! 515! forces and moments in sections 516! These are obtained by integrating the stress tensor; 517! the stresses at the integration points are interpolated 518! from the values at the nodes of the face; these values 519! were extrapolated from the integration point values within 520! the element. This approach is different from the one 521! under "DRAG" because here the stress-strain relationship 522! can be highly nonlinear (plasticity..); this is not the 523! case in cfd. 524! 525! interpolation of the stress tensor at the integration 526! point 527! 528 do i1=1,3 529 do j1=i1,3 530 t(i1,j1)=0.d0 531 enddo 532 enddo 533! 534 if((nope.eq.20).or.(nope.eq.8)) then 535 do i1=1,nopes 536 t(1,1)=t(1,1)+ 537 & shp2(4,i1)*stn(1,konl(ifaceq(i1,ig))) 538 t(2,2)=t(2,2)+ 539 & shp2(4,i1)*stn(2,konl(ifaceq(i1,ig))) 540 t(3,3)=t(3,3)+ 541 & shp2(4,i1)*stn(3,konl(ifaceq(i1,ig))) 542 t(1,2)=t(1,2)+ 543 & shp2(4,i1)*stn(4,konl(ifaceq(i1,ig))) 544 t(1,3)=t(1,3)+ 545 & shp2(4,i1)*stn(5,konl(ifaceq(i1,ig))) 546 t(2,3)=t(2,3)+ 547 & shp2(4,i1)*stn(6,konl(ifaceq(i1,ig))) 548 enddo 549 elseif((nope.eq.10).or.(nope.eq.4)) then 550 do i1=1,nopes 551 t(1,1)=t(1,1)+ 552 & shp2(4,i1)*stn(1,konl(ifacet(i1,ig))) 553 t(2,2)=t(2,2)+ 554 & shp2(4,i1)*stn(2,konl(ifacet(i1,ig))) 555 t(3,3)=t(3,3)+ 556 & shp2(4,i1)*stn(3,konl(ifacet(i1,ig))) 557 t(1,2)=t(1,2)+ 558 & shp2(4,i1)*stn(4,konl(ifacet(i1,ig))) 559 t(1,3)=t(1,3)+ 560 & shp2(4,i1)*stn(5,konl(ifacet(i1,ig))) 561 t(2,3)=t(2,3)+ 562 & shp2(4,i1)*stn(6,konl(ifacet(i1,ig))) 563 enddo 564 else 565 do i1=1,nopes 566 t(1,1)=t(1,1)+ 567 & shp2(4,i1)*stn(1,konl(ifacew(i1,ig))) 568 t(2,2)=t(2,2)+ 569 & shp2(4,i1)*stn(2,konl(ifacew(i1,ig))) 570 t(3,3)=t(3,3)+ 571 & shp2(4,i1)*stn(3,konl(ifacew(i1,ig))) 572 t(1,2)=t(1,2)+ 573 & shp2(4,i1)*stn(4,konl(ifacew(i1,ig))) 574 t(1,3)=t(1,3)+ 575 & shp2(4,i1)*stn(5,konl(ifacew(i1,ig))) 576 t(2,3)=t(2,3)+ 577 & shp2(4,i1)*stn(6,konl(ifacew(i1,ig))) 578 enddo 579 endif 580 t(2,1)=t(1,2) 581 t(3,1)=t(1,3) 582 t(3,2)=t(2,3) 583! 584! calculating the force contribution 585! 586 do i1=1,3 587 df(i1)=(t(i1,1)*xsj2(1)+ 588 & t(i1,2)*xsj2(2)+ 589 & t(i1,3)*xsj2(3))*weight 590 enddo 591! 592! update total force and total moment 593! 594 do i1=1,3 595 f(i1)=f(i1)+df(i1) 596 enddo 597 xm(1)=xm(1)+coords(2)*df(3)-coords(3)*df(2) 598 xm(2)=xm(2)+coords(3)*df(1)-coords(1)*df(3) 599 xm(3)=xm(3)+coords(1)*df(2)-coords(2)*df(1) 600! 601! update area, center of gravity and mean normal 602! 603 dd=dsqrt(xsj2(1)*xsj2(1)+xsj2(2)*xsj2(2)+ 604 & xsj2(3)*xsj2(3)) 605 area=area+dd*weight 606 do i1=1,3 607 cg(i1)=cg(i1)+coords(i1)*dd*weight 608 xn(i1)=xn(i1)+xsj2(i1)*weight 609 enddo 610 endif 611 enddo 612 enddo 613! 614 if(prlab(ii)(1:4).eq.'DRAG') then 615 write(5,*) 616 write(5,122) faset(1:ipos-2),ttime+time 617 122 format(' total surface force (fx,fy,fz) for set ',A, 618 & ' and time ',e14.7) 619 write(5,*) 620 write(5,'(1p,3(1x,e13.6))') (f(j),j=1,3) 621 elseif(prlab(ii)(1:4).eq.'FLUX') then 622 write(5,*) 623 write(5,123) faset(1:ipos-2),ttime+time 624 123 format(' total surface flux (q) for set ',A, 625 & ' and time ',e14.7) 626 write(5,*) 627 write(5,'(1p,1x,e13.6)') f(0) 628 else 629! 630! surface set statistics 631! 632 write(5,*) 633 write(5,130) faset(1:ipos-2),ttime+time 634 130 format(' statistics for surface set ',A, 635 & ' and time ',e14.7) 636! 637! total force and moment about the origin 638! 639 write(5,*) 640 write(5,126) 641! 642 126 format(' total surface force (fx,fy,fz) ', 643 & 'and moment about the origin(mx,my,mz)') 644 write(5,*) 645 write(5,'(2x,1p,6(1x,e13.6))') (f(j),j=1,3),(xm(j),j=1,3) 646! 647! center of gravity and mean normal 648! 649 do i1=1,3 650 cg(i1)=cg(i1)/area 651 xn(i1)=xn(i1)/area 652 enddo 653 write(5,*) 654 write(5,127) 655 127 format(' center of gravity and mean normal') 656 write(5,*) 657 write(5,'(2x,1p,6(1x,e13.6))')(cg(j),j=1,3),(xn(j),j=1,3) 658! 659! moment about the center of gravity 660! 661 write(5,*) 662 write(5,129) 663 129 format( 664 & ' moment about the center of gravity(mx,my,mz)') 665 write(5,*) 666 write(5,'(2x,1p,6(1x,e13.6))') 667 & xm(1)-cg(2)*f(3)+cg(3)*f(2), 668 & xm(2)-cg(3)*f(1)+cg(1)*f(3), 669 & xm(3)-cg(1)*f(2)+cg(2)*f(1) 670! 671! area, normal and shear force 672! 673 xnormforc=f(1)*xn(1)+f(2)*xn(2)+f(3)*xn(3) 674 shearforc=sqrt((f(1)-xnormforc*xn(1))**2+ 675 & (f(2)-xnormforc*xn(2))**2+ 676 & (f(3)-xnormforc*xn(3))**2) 677 write(5,*) 678 write(5,128) 679 128 format(' area, ', 680 & ' normal force (+ = tension) and shear force (size)') 681 write(5,*) 682 write(5,'(2x,1p,3(1x,e13.6))') area,xnormforc,shearforc 683 endif 684! 685 endif 686 enddo 687! 688 return 689 end 690