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 solidsections(inpc,textpart,set,istartset,iendset, 20 & ialset,nset,ielmat,matname,nmat,ielorien,orname,norien, 21 & lakon,thicke,kon,ipkon,irstrt,istep,istat,n,iline,ipol,inl, 22 & ipoinp,inp,cs,mcs,iaxial,ipoinpc,mi,co,ixfree,xnor,iponor, 23 & ier,orab) 24! 25! reading the input deck: *SOLID SECTION 26! 27 implicit none 28! 29 logical nodalthickness 30! 31 character*1 inpc(*) 32 character*8 lakon(*) 33 character*80 matname(*),orname(*),material,orientation 34 character*81 set(*),elset 35 character*132 textpart(16) 36! 37 integer mi(*),istartset(*),iendset(*),ialset(*),ielmat(mi(3),*), 38 & ielorien(mi(3),*),kon(*),ipkon(*),indexe,irstrt(*),nset,nmat, 39 & norien,ielem,node1,node2,m,indexx,ixfree,iponor(2,*),id, 40 & istep,istat,n,key,i,j,k,l,imaterial,iorientation,ipos, 41 & iline,ipol,inl,ipoinp(2,*),inp(3,*),mcs,iaxial,ipoinpc(0:*), 42 & ier,numnod 43! 44 real*8 thicke(mi(3),*),thickness,pi,cs(17,*),xn(3),co(3,*),p(3), 45 & dd,xnor(*),orab(7,*) 46! 47 if((istep.gt.0).and.(irstrt(1).ge.0)) then 48 write(*,*) '*ERROR reading *SOLID SECTION: *SOLID SECTION' 49 write(*,*)' should be placed before all step definitions' 50 ier=1 51 return 52 endif 53! 54 nodalthickness=.false. 55 pi=4.d0*datan(1.d0) 56! 57 orientation=' 58 & ' 59 elset=' 60 & ' 61 ipos=0 62! 63 do i=2,n 64 if(textpart(i)(1:9).eq.'MATERIAL=') then 65 material=textpart(i)(10:89) 66 elseif(textpart(i)(1:12).eq.'ORIENTATION=') then 67 orientation=textpart(i)(13:92) 68 elseif(textpart(i)(1:6).eq.'ELSET=') then 69 elset=textpart(i)(7:86) 70 elset(81:81)=' ' 71 ipos=index(elset,' ') 72 elset(ipos:ipos)='E' 73 elseif(textpart(i)(1:14).eq.'NODALTHICKNESS') then 74 nodalthickness=.true. 75 else 76 write(*,*) 77 & '*WARNING reading *SOLID SECTION: parameter not recognized:' 78 write(*,*) ' ', 79 & textpart(i)(1:index(textpart(i),' ')-1) 80 call inputwarning(inpc,ipoinpc,iline, 81 &"*SOLID SECTION%") 82 endif 83 enddo 84! 85! check for the existence of the material 86! 87 do i=1,nmat 88 if(matname(i).eq.material) exit 89 enddo 90 if(i.gt.nmat) then 91 do i=1,nmat 92 if(matname(i)(1:11).eq.'ANISO_CREEP') then 93 if(matname(i)(12:20).eq.material(1:9)) exit 94 elseif(matname(i)(1:10).eq.'ANISO_PLAS') then 95 if(matname(i)(11:20).eq.material(1:10)) exit 96 endif 97 enddo 98 endif 99 if(i.gt.nmat) then 100 write(*,*) 101 & '*ERROR reading *SOLID SECTION: nonexistent material' 102 write(*,*) ' ' 103 call inputerror(inpc,ipoinpc,iline, 104 & "*SOLID SECTION%",ier) 105 return 106 endif 107 imaterial=i 108! 109! check for the existence of the orientation 110! 111 if(orientation.eq.' ') then 112 iorientation=0 113 else 114 do i=1,norien 115 if(orname(i).eq.orientation) exit 116 enddo 117 if(i.gt.norien) then 118 write(*,*) 119 & '*ERROR reading *SOLID SECTION: nonexistent orientation' 120 write(*,*) ' ' 121 call inputerror(inpc,ipoinpc,iline, 122 & "*SOLID SECTION%",ier) 123 return 124 endif 125 iorientation=i 126 endif 127! 128! check for the existence of the set 129! 130 if(ipos.eq.0) then 131 write(*,*) '*ERROR reading *SOLID SECTION: no element set ', 132 & elset 133 write(*,*) ' was been defined. ' 134 call inputerror(inpc,ipoinpc,iline, 135 & "*SOLID SECTION%",ier) 136 return 137 endif 138c do i=1,nset 139c if(set(i).eq.elset) exit 140c enddo 141 call cident81(set,elset,nset,id) 142 i=nset+1 143 if(id.gt.0) then 144 if(elset.eq.set(id)) then 145 i=id 146 endif 147 endif 148 if(i.gt.nset) then 149 elset(ipos:ipos)=' ' 150 write(*,*) '*ERROR reading *SOLID SECTION: element set ',elset 151 write(*,*) ' has not yet been defined. ' 152 call inputerror(inpc,ipoinpc,iline, 153 & "*SOLID SECTION%",ier) 154 return 155 endif 156! 157! assigning the elements of the set the appropriate material 158! and orientation number 159! 160 do j=istartset(i),iendset(i) 161 if(ialset(j).gt.0) then 162 k=ialset(j) 163 if((lakon(k)(1:1).eq.'B').or. 164 & (lakon(k)(1:1).eq.'S')) then 165 write(*,*) 166 & '*ERROR reading *SOLID SECTION: *SOLID SECTION can' 167 write(*,*) ' not be used for beam or shell elements 168 &' 169 write(*,*) ' Faulty element: ',k 170 ier=1 171 return 172 endif 173 ielmat(1,k)=imaterial 174 if(ielorien(1,k).lt.0) then 175! 176! an orientation based on a distribution has the 177! same name as the distribution. However, to a 178! distribution which is not used in any orientation 179! definition no local coordinate system has been 180! assigned (i.e. orab(7,..)=0.d0) 181! 182 if((orname(-ielorien(1,k)).eq.orientation).and. 183 & (orab(7,-ielorien(1,k)).ne.0.d0)) then 184 ielorien(1,k)=-ielorien(1,k) 185 else 186 ielorien(1,k)=iorientation 187 endif 188 else 189 ielorien(1,k)=iorientation 190 endif 191 else 192 k=ialset(j-2) 193 do 194 k=k-ialset(j) 195 if(k.ge.ialset(j-1)) exit 196 if((lakon(k)(1:1).eq.'B').or. 197 & (lakon(k)(1:1).eq.'S')) then 198 write(*,*) '*ERROR reading *SOLID SECTION: *SOLID SECT 199 &ION can' 200 write(*,*) ' not be used for beam or shell eleme 201 &nts' 202 write(*,*) ' Faulty element: ',k 203 ier=1 204 return 205 endif 206 ielmat(1,k)=imaterial 207 if(ielorien(1,k).lt.0) then 208! 209! an orientation based on a distribution has the 210! same name as the distribution. However, to a 211! distribution which is not used in any orientation 212! definition no local coordinate system has been 213! assigned (i.e. orab(7,..)=0.d0) 214! 215 if((orname(-ielorien(1,k)).eq.orientation).and. 216 & (orab(7,-ielorien(1,k)).ne.0.d0)) then 217 ielorien(1,k)=-ielorien(1,k) 218 else 219 ielorien(1,k)=iorientation 220 endif 221 else 222 ielorien(1,k)=iorientation 223 endif 224 enddo 225 endif 226 enddo 227! 228 call getnewline(inpc,textpart,istat,n,key,iline,ipol,inl, 229 & ipoinp,inp,ipoinpc) 230 if(istat.lt.0) return 231! 232! assigning a thickness to plane stress/strain elements and an angle to 233! axisymmetric elements 234! 235 if(key.eq.0) then 236! 237! second line with thickness given 238! 239 read(textpart(1)(1:20),'(f20.0)',iostat=istat) thickness 240 if(istat.gt.0) then 241 call inputerror(inpc,ipoinpc,iline, 242 & "*SOLID SECTION%",ier) 243 return 244 endif 245! 246! for axial symmetric structures: 247! thickness for axial symmetric elements: 2 degrees 248! thickness for plane stress elements: reduced by 180 249! thickness for plane strain elements: reduced by 180 250! 251 if(.not.nodalthickness) then 252 if(iaxial.eq.180) then 253 if(lakon(ialset(istartset(i)))(1:2).eq.'CA') then 254 thickness=datan(1.d0)*8.d0/iaxial 255 elseif(lakon(ialset(istartset(i)))(1:3).eq.'CPS') then 256 thickness=thickness/iaxial 257 elseif(lakon(ialset(istartset(i)))(1:3).eq.'CPE') then 258 thickness=thickness/iaxial 259 endif 260 endif 261 else 262! 263! for those elements for which nodal thickness is activated 264! the thickness is set to -1.d0 265! 266 thickness=-1.d0 267 endif 268! 269! assigning the thickness to each node of the corresponding 270! elements (thickness specified) 271! 272 do j=istartset(i),iendset(i) 273 if(ialset(j).gt.0) then 274 if((lakon(ialset(j))(1:2).eq.'CP').or. 275 & (lakon(ialset(j))(1:2).eq.'CA')) then 276! 277! plane stress/strain or axisymmetric elements 278! 279 indexe=ipkon(ialset(j)) 280 read(lakon(ialset(j))(4:4),'(i1)') numnod 281 do l=1,numnod 282 thicke(1,indexe+l)=thickness 283 enddo 284 elseif(lakon(ialset(j))(1:1).eq.'T') then 285 ielem=ialset(j) 286! 287! default cross section for trusses is the 288! rectangular cross section 289! 290 lakon(ielem)(8:8)='R' 291 indexe=ipkon(ielem) 292 node1=kon(indexe+1) 293 if(lakon(ielem)(4:4).eq.'2') then 294 node2=kon(indexe+2) 295 else 296 node2=kon(indexe+3) 297 endif 298! 299! determining a vector orthogonal to the truss 300! element 301! 302 do l=1,3 303 xn(l)=co(l,node2)-co(l,node1) 304 enddo 305 if(dabs(xn(1)).gt.0.d0) then 306 p(1)=-xn(3) 307 p(2)=0.d0 308 p(3)=xn(1) 309 elseif(dabs(xn(2)).gt.0.d0) then 310 p(1)=xn(2) 311 p(2)=-xn(1) 312 p(3)=0.d0 313 else 314 p(1)=0.d0 315 p(2)=xn(3) 316 p(3)=-xn(2) 317 endif 318 dd=dsqrt(p(1)*p(1)+p(2)*p(2)+p(3)*p(3)) 319 if(dd.lt.1.d-10) then 320 write(*,*) 321 & '*ERROR reading *SOLID SECTION: normal' 322 write(*,*) ' in direction 1 has zero size' 323 ier=1 324 return 325 endif 326 do l=1,3 327 p(l)=p(l)/dd 328 enddo 329 do l=1,3 330 thicke(1,indexe+l)=dsqrt(thickness) 331 thicke(2,indexe+l)=dsqrt(thickness) 332 indexx=ixfree 333 do m=1,3 334 xnor(indexx+m)=p(m) 335 enddo 336 ixfree=ixfree+6 337 iponor(1,indexe+l)=indexx 338 enddo 339 endif 340 else 341 k=ialset(j-2) 342 do 343 k=k-ialset(j) 344 if(k.ge.ialset(j-1)) exit 345 if((lakon(k)(1:2).eq.'CP').or. 346 & (lakon(k)(1:2).eq.'CA')) then 347 indexe=ipkon(k) 348 read(lakon(k)(4:4),'(i1)') numnod 349 do l=1,numnod 350 thicke(1,indexe+l)=thickness 351 enddo 352 elseif(lakon(k)(1:1).eq.'T') then 353! 354! default cross section for trusses is the 355! rectangular cross section 356! 357 lakon(k)(8:8)='R' 358 indexe=ipkon(k) 359 node1=kon(indexe+1) 360 if(lakon(k)(4:4).eq.'2') then 361 node2=kon(indexe+2) 362 else 363 node2=kon(indexe+3) 364 endif 365! 366! determining a vector orthogonal to the truss 367! element 368! 369 do l=1,3 370 xn(l)=co(l,node2)-co(l,node1) 371 enddo 372 if(dabs(xn(1)).gt.0.d0) then 373 p(1)=-xn(3) 374 p(2)=0.d0 375 p(3)=xn(1) 376 elseif(dabs(xn(2)).gt.0.d0) then 377 p(1)=xn(2) 378 p(2)=-xn(1) 379 p(3)=0.d0 380 else 381 p(1)=0.d0 382 p(2)=xn(3) 383 p(3)=-xn(2) 384 endif 385 dd=dsqrt(p(1)*p(1)+p(2)*p(2)+p(3)*p(3)) 386 if(dd.lt.1.d-10) then 387 write(*,*) 388 & '*ERROR reading *SOLID SECTION: normal' 389 write(*,*) ' in direction 1 has zero size' 390 ier=1 391 return 392 endif 393 do l=1,3 394 p(l)=p(l)/dd 395 enddo 396 do l=1,3 397 thicke(1,indexe+l)=dsqrt(thickness) 398 thicke(2,indexe+l)=dsqrt(thickness) 399 indexx=ixfree 400 do m=1,3 401 xnor(indexx+m)=p(m) 402 enddo 403 ixfree=ixfree+6 404 iponor(1,indexe+l)=indexx 405 enddo 406 endif 407 enddo 408 endif 409 enddo 410 else 411! 412! no second line (no thickness given) 413! 414! assigning the thickness to each node of the corresponding 415! elements (thickness not specified: only axisymmetric elements 416! or plane stress/strain elements with nodal thickness) 417! 418 thickness=datan(1.d0)*8.d0/iaxial 419 do j=istartset(i),iendset(i) 420 if(ialset(j).gt.0) then 421 if(lakon(ialset(j))(1:2).eq.'CA') then 422! 423! axisymmetric elements 424! 425 indexe=ipkon(ialset(j)) 426 read(lakon(ialset(j))(4:4),'(i1)') numnod 427 do l=1,numnod 428 thicke(1,indexe+l)=thickness 429 enddo 430 elseif((lakon(ialset(j))(1:3).eq.'CPS').or. 431 & (lakon(ialset(j))(1:3).eq.'CPE')) then 432 if(nodalthickness) then 433 indexe=ipkon(ialset(j)) 434 read(lakon(ialset(j))(4:4),'(i1)') numnod 435 do l=1,numnod 436 thicke(1,indexe+l)=-1.d0 437 enddo 438 endif 439 endif 440 else 441 k=ialset(j-2) 442 do 443 k=k-ialset(j) 444 if(k.ge.ialset(j-1)) exit 445 if(lakon(k)(1:2).eq.'CA') then 446! 447! axisymmetric elements 448! 449 indexe=ipkon(k) 450 read(lakon(k)(4:4),'(i1)') numnod 451 do l=1,numnod 452 thicke(1,indexe+l)=thickness 453 enddo 454 elseif((lakon(k)(1:3).eq.'CPS').or. 455 & (lakon(k)(1:3).eq.'CPE')) then 456 if(nodalthickness) then 457 indexe=ipkon(k) 458 read(lakon(k)(4:4),'(i1)') numnod 459 do l=1,numnod 460 thicke(1,indexe+l)=-1.d0 461 enddo 462 endif 463 endif 464 enddo 465 endif 466 enddo 467 endif 468! 469! defining cyclic symmetric conditions for axisymmetric 470! elements (needed for cavity radiation) 471! 472 do j=istartset(i),iendset(i) 473 if(ialset(j).gt.0) then 474 if(lakon(ialset(j))(1:2).eq.'CA') then 475 if(mcs.gt.1) then 476 write(*,*) '*ERROR reading *SOLID SECTION: ' 477 write(*,*) ' axisymmetric elements cannot be 478 &combined with cyclic symmetry' 479 ier=1 480 return 481 elseif(mcs.eq.1) then 482 if(int(cs(1,1)).ne.int(2.d0*pi/thickness+0.5d0)) 483 & then 484 write(*,*) '*ERROR reading *SOLID SECTION: ' 485 write(*,*) ' it is not allowed to define t 486 &wo different' 487 write(*,*) ' angles for an axisymmetric st 488 &ructure' 489 ier=1 490 return 491 else 492 exit 493 endif 494 endif 495 mcs=1 496 cs(1,1)=2.d0*pi/thickness+0.5d0 497 cs(2,1)=-0.5d0 498 cs(3,1)=-0.5d0 499 cs(5,1)=1.5d0 500 do k=6,9 501 cs(k,1)=0.d0 502 enddo 503 cs(10,1)=1.d0 504 cs(11,1)=0.d0 505 cs(12,1)=-1.d0 506 cs(14,1)=0.5d0 507 cs(15,1)=dcos(thickness) 508 cs(16,1)=dsin(thickness) 509 exit 510 endif 511 else 512 k=ialset(j-2) 513 do 514 k=k-ialset(j) 515 if(k.ge.ialset(j-1)) exit 516 if(lakon(k)(1:2).eq.'CA') then 517 if(mcs.gt.1) then 518 write(*,*) '*ERROR reading *SOLID SECTION: ' 519 write(*,*) ' axisymmetric elements cannot 520 &be combined with cyclic symmetry' 521 ier=1 522 return 523 elseif(mcs.eq.1) then 524 if(int(cs(1,1)).ne.int(2.d0*pi/thickness+0.5d0)) 525 & then 526 write(*,*) '*ERROR reading *SOLID SECTION: ' 527 write(*,*) ' it is not allowed to defin 528 &e two different' 529 write(*,*) ' angles for an axisymmetric 530 &structure' 531 ier=1 532 return 533 else 534 exit 535 endif 536 endif 537 mcs=1 538 cs(1,1)=2.d0*pi/thickness+0.5d0 539 cs(2,1)=-0.5d0 540 cs(3,1)=-0.5d0 541 cs(5,1)=1.5d0 542 do k=6,9 543 cs(k,1)=0.d0 544 enddo 545 cs(10,1)=1.d0 546 cs(11,1)=0.d0 547 cs(12,1)=-1.d0 548 cs(14,1)=0.5d0 549 cs(15,1)=dcos(thickness) 550 cs(16,1)=dsin(thickness) 551 exit 552 endif 553 enddo 554 endif 555 enddo 556! 557 if(key.eq.0) then 558 call getnewline(inpc,textpart,istat,n,key,iline,ipol,inl, 559 & ipoinp,inp,ipoinpc) 560 endif 561! 562 return 563 end 564 565