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 dloads(inpc,textpart,set,istartset,iendset, 20 & ialset,nset,nelemload,sideload,xload,nload,nload_, 21 & ielmat,iamload,amname,nam,lakon,ne,dload_flag,istep, 22 & istat,n,iline,ipol,inl,ipoinp,inp,cbody,ibody,xbody,nbody, 23 & nbody_,xbodyold,iperturb,physcon,nam_,namtot_,namta,amta, 24 & nmethod,ipoinpc,maxsectors,mi,idefload,idefbody,ipkon, 25 & thicke,iamplitudedefault,namtot,ier) 26! 27! reading the input deck: *DLOAD 28! 29 implicit none 30! 31 logical dload_flag,submodel,edgeload,surface 32! 33 character*1 inpc(*) 34 character*8 lakon(*) 35 character*20 sideload(*),label 36 character*80 amname(*),amplitude 37 character*81 set(*),elset,cbody(*) 38 character*132 textpart(16) 39! 40 integer istartset(*),iendset(*),ialset(*),nelemload(2,*),mi(*), 41 & ielmat(mi(3),*),nset,nload,nload_,istep,istat,n,i,j,l,key, 42 & iamload(2,*),nam,iamplitude,ipos,ne,iline,ipol,iperturb(*), 43 & inl,ipoinp(2,*),inp(3,*),ibody(3,*),nbody,nbody_,nam_,namtot, 44 & namtot_,namta(3,*),idelay,nmethod,lc,isector,node,id, 45 & ipoinpc(0:*),maxsectors,jsector,iglobstep,idefload(*), 46 % idefbody(*),ipkon(*),k,indexe,iamplitudedefault,ier 47! 48 real*8 xload(2,*),xbody(7,*),xmagnitude,dd,p1(3),p2(3),bodyf(3), 49 & xbodyold(7,*),physcon(*),amta(2,*),xxmagnitude, 50 & thicke(mi(3),*),thickness 51! 52 iamplitude=iamplitudedefault 53 idelay=0 54 lc=1 55 isector=0 56 submodel=.false. 57 iglobstep=0 58 edgeload=.false. 59 surface=.false. 60! 61 if(istep.lt.1) then 62 write(*,*) '*ERROR reading *DLOAD: *DLOAD should only be used' 63 write(*,*) ' within a STEP' 64 ier=1 65 return 66 endif 67! 68 do i=2,n 69 if((textpart(i)(1:6).eq.'OP=NEW').and.(.not.dload_flag)) then 70 do j=1,nload 71 if(sideload(j)(1:1).eq.'P') then 72 sideload(j)(3:4)=' ' 73 xload(1,j)=0.d0 74 endif 75 enddo 76 do j=1,nbody 77 xbody(1,j)=0.d0 78 enddo 79 elseif(textpart(i)(1:10).eq.'AMPLITUDE=') then 80 read(textpart(i)(11:90),'(a80)') amplitude 81 do j=1,nam 82 if(amname(j).eq.amplitude) then 83 iamplitude=j 84 exit 85 endif 86 enddo 87 if(j.gt.nam) then 88 write(*,*)'*ERROR reading *DLOAD: nonexistent amplitude' 89 write(*,*) ' ' 90 call inputerror(inpc,ipoinpc,iline, 91 & "*DLOAD%",ier) 92 return 93 endif 94 iamplitude=j 95 elseif(textpart(i)(1:10).eq.'TIMEDELAY=') THEN 96 if(idelay.ne.0) then 97 write(*,*) 98 & '*ERROR reading *DLOAD: the parameter TIME DELAY' 99 write(*,*) ' is used twice in the same keyword' 100 write(*,*) ' ' 101 call inputerror(inpc,ipoinpc,iline, 102 & "*DLOAD%",ier) 103 return 104 else 105 idelay=1 106 endif 107 nam=nam+1 108 if(nam.gt.nam_) then 109 write(*,*) '*ERROR reading *DLOAD: increase nam_' 110 ier=1 111 return 112 endif 113 amname(nam)=' 114 &' 115 if(iamplitude.eq.0) then 116 write(*,*) '*ERROR reading *DLOAD: time delay must be' 117 write(*,*) ' preceded by the amplitude parameter' 118 ier=1 119 return 120 endif 121 namta(3,nam)=sign(iamplitude,namta(3,iamplitude)) 122 iamplitude=nam 123c if(nam.eq.1) then 124c namtot=0 125c else 126c namtot=namta(2,nam-1) 127c endif 128 namtot=namtot+1 129 if(namtot.gt.namtot_) then 130 write(*,*) '*ERROR dloads: increase namtot_' 131 ier=1 132 return 133 endif 134 namta(1,nam)=namtot 135 namta(2,nam)=namtot 136c call reorderampl(amname,namta,nam) 137 read(textpart(i)(11:30),'(f20.0)',iostat=istat) 138 & amta(1,namtot) 139 if(istat.gt.0) then 140 call inputerror(inpc,ipoinpc,iline, 141 & "*DLOAD%",ier) 142 return 143 endif 144 elseif(textpart(i)(1:9).eq.'LOADCASE=') then 145 read(textpart(i)(10:19),'(i10)',iostat=istat) lc 146 if(istat.gt.0) then 147 call inputerror(inpc,ipoinpc,iline, 148 & "*DLOAD%",ier) 149 return 150 endif 151 if(nmethod.ne.5) then 152 write(*,*) 153 & '*ERROR reading *DLOAD: the parameter LOAD CASE' 154 write(*,*) ' is only allowed in STEADY STATE' 155 write(*,*) ' DYNAMICS calculations' 156 ier=1 157 return 158 endif 159 elseif(textpart(i)(1:7).eq.'SECTOR=') then 160 read(textpart(i)(8:17),'(i10)',iostat=istat) isector 161 if(istat.gt.0) then 162 call inputerror(inpc,ipoinpc,iline, 163 & "*DLOAD%",ier) 164 return 165 endif 166 if((nmethod.le.3).or.(iperturb(1).gt.1)) then 167 write(*,*) '*ERROR reading *DLOAD: the parameter SECTOR' 168 write(*,*) ' is only allowed in MODAL DYNAMICS or' 169 write(*,*) ' STEADY STATE DYNAMICS calculations' 170 ier=1 171 return 172 endif 173 if(isector.gt.maxsectors) then 174 write(*,*) '*ERROR reading *DLOAD: sector ',isector 175 write(*,*) ' exceeds number of sectors' 176 ier=1 177 return 178 endif 179 isector=isector-1 180 elseif(textpart(i)(1:8).eq.'SUBMODEL') then 181 submodel=.true. 182 elseif(textpart(i)(1:5).eq.'STEP=') then 183 read(textpart(i)(6:15),'(i10)',iostat=istat) iglobstep 184 if(istat.gt.0) then 185 call inputerror(inpc,ipoinpc,iline, 186 & "*DLOAD%",ier) 187 return 188 endif 189 elseif(textpart(i)(1:8).eq.'DATASET=') then 190 read(textpart(i)(9:18),'(i10)',iostat=istat) iglobstep 191 if(istat.gt.0) then 192 call inputerror(inpc,ipoinpc,iline, 193 & "*DLOAD%",ier) 194 return 195 endif 196! 197! the mode number for submodels 198! is stored as a negative global step 199! 200 iglobstep=-iglobstep 201 else 202 write(*,*) 203 & '*WARNING reading *DLOAD: parameter not recognized:' 204 write(*,*) ' ', 205 & textpart(i)(1:index(textpart(i),' ')-1) 206 call inputwarning(inpc,ipoinpc,iline, 207 & "*DLOAD%") 208 endif 209 enddo 210! 211! check for the presence of an amplitude in submodel cases 212! 213 if(submodel) then 214 if(iamplitude.ne.0) then 215 write(*,*) '*WARNING reading *DSLOAD:' 216 write(*,*) ' no amplitude definition is allowed' 217 write(*,*) ' in combination with a submodel' 218 endif 219 endif 220! 221! check whether global step was specified for submodel 222! 223 if((submodel).and.(iglobstep.eq.0)) then 224 write(*,*) '*ERROR reading *DLOAD: no global step' 225 write(*,*) ' step specified for the submodel' 226 call inputerror(inpc,ipoinpc,iline, 227 & "*DLOAD%",ier) 228 return 229 endif 230! 231 do 232 call getnewline(inpc,textpart,istat,n,key,iline,ipol,inl, 233 & ipoinp,inp,ipoinpc) 234 if((istat.lt.0).or.(key.eq.1)) return 235! 236 read(textpart(2)(1:20),'(a20)',iostat=istat) label 237! 238! for submodels the load label is modified and the global 239! step is stored in iamload(1,*) 240! 241 if(submodel) then 242 label(3:4)='SM' 243 iamplitude=iglobstep 244 endif 245! 246 if(label(3:4).ne.'NP') then 247 read(textpart(3)(1:20),'(f20.0)',iostat=istat) xmagnitude 248 else 249 read(textpart(3)(1:10),'(i10)',iostat=istat) node 250 endif 251 if(istat.gt.0) then 252 call inputerror(inpc,ipoinpc,iline, 253 & "*DLOAD%",ier) 254 return 255 endif 256 if(label(1:7).eq.'CENTRIF') then 257 do i=1,3 258 read(textpart(i+3)(1:20),'(f20.0)',iostat=istat) p1(i) 259 if(istat.gt.0) then 260 call inputerror(inpc,ipoinpc,iline, 261 & "*DLOAD%",ier) 262 return 263 endif 264 enddo 265 do i=1,3 266 read(textpart(i+6)(1:20),'(f20.0)',iostat=istat) p2(i) 267 if(istat.gt.0) then 268 call inputerror(inpc,ipoinpc,iline, 269 & "*DLOAD%",ier) 270 return 271 endif 272 enddo 273 dd=dsqrt(p2(1)**2+p2(2)**2+p2(3)**2) 274 do i=1,3 275 p2(i)=p2(i)/dd 276 enddo 277 elseif(label(1:4).eq.'GRAV') then 278 do i=1,3 279 read(textpart(i+3)(1:20),'(f20.0)',iostat=istat) bodyf(i) 280 if(istat.gt.0) then 281 call inputerror(inpc,ipoinpc,iline, 282 & "*DLOAD%",ier) 283 return 284 endif 285 enddo 286 elseif(label(1:6).eq.'NEWTON') then 287 if(iperturb(1).le.1) then 288 write(*,*) '*ERROR reading *DLOAD: NEWTON gravity force' 289 write(*,*) ' can only be used in a nonlinear' 290 write(*,*) ' procedure' 291 ier=1 292 return 293 endif 294 if(physcon(3).le.0.d0) then 295 write(*,*) '*ERROR reading *DLOAD: NEWTON gravity force' 296 write(*,*) ' requires the definition of a' 297 write(*,*) ' positive gravity constant with' 298 write(*,*) ' a *PHYSICAL CONSTANTS card' 299 ier=1 300 return 301 endif 302 elseif(((label(1:2).ne.'P1').and.(label(1:2).ne.'P2').and. 303 & (label(1:2).ne.'P3').and.(label(1:2).ne.'P4').and. 304 & (label(1:2).ne.'P5').and.(label(1:2).ne.'P6').and. 305 & (label(1:2).ne.'P ').and.(label(1:2).ne.'BX').and. 306 & (label(1:2).ne.'BY').and.(label(1:2).ne.'BZ').and. 307c BernhardiStart 308 & (label(1:2).ne.'ED')).or. 309 & ((label(3:6).ne.'NOR1').and.(label(3:6).ne.'NOR2').and. 310 & (label(3:6).ne.'NOR3').and.(label(3:6).ne.'NOR4')).and. 311c BernhardiEnd 312 & ((label(3:4).ne.' ').and.(label(3:4).ne.'NU').and. 313 & (label(3:4).ne.'NP').and.(label(3:4).ne.'SM'))) then 314 call inputerror(inpc,ipoinpc,iline, 315 & "*DLOAD%",ier) 316 return 317 endif 318! 319 read(textpart(1)(1:10),'(i10)',iostat=istat) l 320 if(istat.eq.0) then 321 if(l.gt.ne) then 322 write(*,*) '*ERROR reading *DLOAD: element ',l 323 write(*,*) ' is not defined' 324 ier=1 325 return 326 endif 327 if((label(1:7).eq.'CENTRIF').or.(label(1:4).eq.'GRAV').or. 328 & (label(1:6).eq.'NEWTON')) then 329 elset(1:80)=textpart(1)(1:80) 330 elset(81:81)=' ' 331 call bodyadd(cbody,ibody,xbody,nbody,nbody_,elset,label, 332 & iamplitude,xmagnitude,p1,p2,bodyf,xbodyold,lc,idefbody) 333 else 334 xxmagnitude=xmagnitude 335 if((lakon(l)(1:2).eq.'CP').or. 336 & (lakon(l)(2:2).eq.'A').or. 337 & (lakon(l)(7:7).eq.'E').or. 338 & (lakon(l)(7:7).eq.'S').or. 339 & (lakon(l)(7:7).eq.'A')) then 340 if(label(1:2).eq.'P1') then 341 label(1:2)='P3' 342 elseif(label(1:2).eq.'P2') then 343 label(1:2)='P4' 344 elseif(label(1:2).eq.'P3') then 345 label(1:2)='P5' 346 elseif(label(1:2).eq.'P4') then 347 label(1:2)='P6' 348 endif 349 elseif((lakon(l)(1:1).eq.'B').or. 350 & (lakon(l)(7:7).eq.'B')) then 351 if(label(1:2).eq.'P2') label(1:2)='P5' 352 elseif((lakon(l)(1:1).eq.'S').or. 353 & (lakon(l)(7:7).eq.'L')) then 354c BernhardiStart 355 if(label(1:6).eq.'EDNOR1') then 356 label(1:2)='P3' 357 edgeload=.true. 358 elseif(label(1:6).eq.'EDNOR2') then 359 label(1:2)='P4' 360 edgeload=.true. 361 elseif(label(1:6).eq.'EDNOR3') then 362 label(1:2)='P5' 363 edgeload=.true. 364 elseif(label(1:6).eq.'EDNOR4') then 365 label(1:2)='P6' 366 edgeload=.true. 367 else 368 label(1:2)='P1' 369 endif 370! 371! EDNOR is an edge load 372! 373 if(edgeload) then 374 indexe=ipkon(l) 375 thickness=0.d0 376 do k=1,mi(3) 377 if(ielmat(k,l).ne.0) then 378 thickness=thickness+thicke(k,indexe+1) 379 else 380 exit 381 endif 382 enddo 383 xxmagnitude=xmagnitude/thickness 384 endif 385c BernhardiEnd 386 endif 387 if(lc.ne.1) then 388 jsector=isector+maxsectors 389 else 390 jsector=isector 391 endif 392 if(label(3:4).ne.'NP') then 393 call loadadd(l,label,xxmagnitude,nelemload,sideload, 394 & xload,nload,nload_,iamload,iamplitude, 395 & nam,jsector,idefload) 396 else 397 call loadaddp(l,label,nelemload,sideload, 398 & xload,nload,nload_,iamload,iamplitude, 399 & nam,node) 400 endif 401 endif 402 else 403 read(textpart(1)(1:80),'(a80)',iostat=istat) elset 404 elset(81:81)=' ' 405 ipos=index(elset,' ') 406! 407! check for element set 408! 409 elset(ipos:ipos)='E' 410c do i=1,nset 411c if(set(i).eq.elset) exit 412c enddo 413 call cident81(set,elset,nset,id) 414 i=nset+1 415 if(id.gt.0) then 416 if(elset.eq.set(id)) then 417 i=id 418 endif 419 endif 420 if(i.gt.nset) then 421! 422! check for facial surface 423! 424 surface=.true. 425 elset(ipos:ipos)='T' 426c do i=1,nset 427c if(set(i).eq.elset) exit 428c enddo 429 call cident81(set,elset,nset,id) 430 i=nset+1 431 if(id.gt.0) then 432 if(elset.eq.set(id)) then 433 i=id 434 endif 435 endif 436 if(i.gt.nset) then 437 elset(ipos:ipos)=' ' 438 write(*,*) '*ERROR reading *DLOAD: element set ' 439 write(*,*) ' or facial surface ',elset 440 write(*,*) ' has not yet been defined. ' 441 call inputerror(inpc,ipoinpc,iline, 442 & "*DLOAD%",ier) 443 return 444 endif 445 endif 446! 447 if((label(1:7).eq.'CENTRIF').or.(label(1:4).eq.'GRAV').or. 448 & (label(1:6).eq.'NEWTON')) then 449 call bodyadd(cbody,ibody,xbody,nbody,nbody_,elset,label, 450 & iamplitude,xmagnitude,p1,p2,bodyf,xbodyold,lc,idefbody) 451 else 452 l=ialset(istartset(i)) 453 if(surface) then 454 write(label(2:2),'(i1)') l-10*(l/10) 455 l=l/10 456 endif 457 if((lakon(l)(1:2).eq.'CP').or. 458 & (lakon(l)(2:2).eq.'A').or. 459 & (lakon(l)(7:7).eq.'E').or. 460 & (lakon(l)(7:7).eq.'S').or. 461 & (lakon(l)(7:7).eq.'A')) then 462 if(label(1:2).eq.'P1') then 463 label(1:2)='P3' 464 elseif(label(1:2).eq.'P2') then 465 label(1:2)='P4' 466 elseif(label(1:2).eq.'P3') then 467 label(1:2)='P5' 468 elseif(label(1:2).eq.'P4') then 469 label(1:2)='P6' 470 endif 471 elseif((lakon(l)(1:1).eq.'B').or. 472 & (lakon(l)(7:7).eq.'B')) then 473 if(label(1:2).eq.'P2') label(1:2)='P5' 474 elseif((lakon(l)(1:1).eq.'S').or. 475 & (lakon(l)(7:7).eq.'L')) then 476c BernhardiStart 477 if(label(1:6).eq.'EDNOR1') then 478 label(1:2)='P3' 479 edgeload=.true. 480 elseif(label(1:6).eq.'EDNOR2') then 481 label(1:2)='P4' 482 edgeload=.true. 483 elseif(label(1:6).eq.'EDNOR3') then 484 label(1:2)='P5' 485 edgeload=.true. 486 elseif(label(1:6).eq.'EDNOR4') then 487 label(1:2)='P6' 488 edgeload=.true. 489 else 490 label(1:2)='P1' 491 endif 492c BernhardiEnd 493 endif 494! 495 do j=istartset(i),iendset(i) 496 if(ialset(j).gt.0) then 497 l=ialset(j) 498 if(surface) then 499 write(label(2:2),'(i1)') l-10*(l/10) 500 l=l/10 501 endif 502 xxmagnitude=xmagnitude 503! 504! EDNOR is an edge load 505! 506 if(edgeload) then 507 indexe=ipkon(l) 508 thickness=0.d0 509 do k=1,mi(3) 510 if(ielmat(k,l).ne.0) then 511 thickness=thickness+thicke(k,indexe+1) 512 else 513 exit 514 endif 515 enddo 516 xxmagnitude=xmagnitude/thickness 517 endif 518! 519 if(lc.ne.1) then 520 jsector=isector+maxsectors 521 else 522 jsector=isector 523 endif 524 if(label(3:4).ne.'NP') then 525 call loadadd(l,label,xxmagnitude,nelemload, 526 & sideload,xload,nload,nload_,iamload, 527 & iamplitude,nam,jsector,idefload) 528 else 529 call loadaddp(l,label,nelemload, 530 & sideload,xload,nload,nload_,iamload, 531 & iamplitude,nam,node) 532 endif 533 else 534 l=ialset(j-2) 535 do 536 l=l-ialset(j) 537 if(l.ge.ialset(j-1)) exit 538 xxmagnitude=xmagnitude 539! 540! EDNOR is an edge load 541! 542 if(edgeload) then 543 indexe=ipkon(l) 544 thickness=0.d0 545 do k=1,mi(3) 546 if(ielmat(k,l).ne.0) then 547 thickness=thickness+thicke(k,indexe+1) 548 else 549 exit 550 endif 551 enddo 552 xxmagnitude=xmagnitude/thickness 553 endif 554! 555 if(lc.ne.1) then 556 jsector=isector+maxsectors 557 else 558 jsector=isector 559 endif 560 if(label(3:4).ne.'NP') then 561 call loadadd(l,label,xxmagnitude,nelemload, 562 & sideload,xload,nload,nload_, 563 & iamload,iamplitude,nam,jsector,idefload) 564 else 565 call loadaddp(l,label,nelemload, 566 & sideload,xload,nload,nload_, 567 & iamload,iamplitude,nam,node) 568 endif 569 enddo 570 endif 571 enddo 572 endif 573 endif 574 enddo 575! 576 return 577 end 578