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 initialconditionss(inpc,textpart,set,istartset,iendset, 20 & ialset,nset,t0,t1,prestr,iprestr,ithermal,veold,inoelfree, 21 & nk_, 22 & mi,istep,istat,n,iline,ipol,inl,ipoinp,inp,lakon,kon,co,ne, 23 & ipkon,vold,ipoinpc,xstate,nstate_,nk,t0g,t1g,iaxial,ielprop, 24 & prop,ier,nuel_) 25! 26! reading the input deck: *INITIAL CONDITIONS 27! 28 implicit none 29! 30 logical user 31! 32 character*1 inpc(*) 33 character*8 lakon(*) 34 character*80 rebarn 35 character*81 set(*),noset 36 character*132 textpart(16) 37! 38 integer istartset(*),iendset(*),ialset(*),nset,iprestr, 39 & ithermal(*),m,id,nuel_, 40 & istep,istat,n,i,j,k,l,ii,key,idir,ipos,inoelfree,nk_,mi(*), 41 & iline,ipol,inl,ipoinp(2,*),inp(3,*),ij,jj,ntens,ncrds,layer, 42 & kspt,lrebar,iflag,i1,mint3d,nope,kon(*),konl(20),indexe, 43 & ipkon(*),ne,ipoinpc(0:*),nstate_,nk,jmax,ntot,numberoflines, 44 & iaxial,null,ielprop(*),ier 45! 46 real*8 t0(*),t1(*),beta(8),prestr(6,mi(1),*),veold(0:mi(2),*), 47 & temperature,velocity,tempgrad1,tempgrad2,pgauss(3), 48 & shp(4,20),xsj,xl(3,20),xi,et,ze,weight,co(3,*),pressure, 49 & vold(0:mi(2),*),xstate(nstate_,mi(1),*),dispvelo,totpres, 50 & xmassflow,t0g(2,*),t1g(2,*),prop(*),turbini(0:mi(2)) 51! 52 include "gauss.f" 53! 54 null=0 55! 56 if(istep.gt.0) then 57 write(*,*) 58 & '*ERROR reading *INITIAL CONDITIONS: *INITIAL CONDITIONS' 59 write(*,*) ' should be placed before all step definitions' 60 ier=1 61 return 62 endif 63! 64 do ij=2,n 65 if(textpart(ij)(1:16).eq.'TYPE=TEMPERATURE') then 66! 67 ithermal(1)=1 68 do 69 call getnewline(inpc,textpart,istat,n,key,iline,ipol,inl, 70 & ipoinp,inp,ipoinpc) 71 if((istat.lt.0).or.(key.eq.1)) return 72 read(textpart(2)(1:20),'(f20.0)',iostat=istat) 73 & temperature 74 if(istat.gt.0) then 75 call inputerror(inpc,ipoinpc,iline, 76 & "*INITIAL CONDITIONS%",ier) 77 return 78 endif 79 temperature=1.d-6*int(1.d6*temperature+0.5d0) 80! 81 if((inoelfree.ne.0).or.(nuel_.gt.0)) then 82 tempgrad1=0.d0 83 tempgrad2=0.d0 84 if(n.gt.2) then 85 read(textpart(3)(1:20),'(f20.0)',iostat=istat) 86 & tempgrad1 87 if(istat.gt.0) then 88 call inputerror(inpc,ipoinpc,iline, 89 & "*INITIAL CONDITIONS%",ier) 90 return 91 endif 92 endif 93 if(n.gt.3) then 94 read(textpart(4)(1:20),'(f20.0)',iostat=istat) 95 & tempgrad2 96 if(istat.gt.0) then 97 call inputerror(inpc,ipoinpc,iline, 98 & "*INITIAL CONDITIONS%",ier) 99 return 100 endif 101 endif 102 endif 103! 104 read(textpart(1)(1:10),'(i10)',iostat=istat) l 105 if(istat.eq.0) then 106 if(l.gt.nk) then 107 write(*,*) 108 & '*WARNING reading *INITIAL CONDITIONS: node ',l 109 write(*,*)' exceeds the largest defined ', 110 & 'node number' 111 cycle 112 endif 113 t0(l)=temperature 114 t1(l)=temperature 115 vold(0,l)=temperature 116 if((inoelfree.ne.0).or.(nuel_.gt.0)) then 117 t0g(1,l)=tempgrad1 118 t0g(2,l)=tempgrad2 119 t1g(1,l)=tempgrad1 120 t1g(2,l)=tempgrad2 121 endif 122 else 123 read(textpart(1)(1:80),'(a80)',iostat=istat) noset 124 noset(81:81)=' ' 125 ipos=index(noset,' ') 126 noset(ipos:ipos)='N' 127c do ii=1,nset 128c if(set(ii).eq.noset) exit 129c enddo 130 call cident81(set,noset,nset,id) 131 ii=nset+1 132 if(id.gt.0) then 133 if(noset.eq.set(id)) then 134 ii=id 135 endif 136 endif 137 if(ii.gt.nset) then 138 noset(ipos:ipos)=' ' 139 write(*,*) 140 & '*ERROR reading *INITIAL CONDITIONS: node set ' 141 & ,noset 142 write(*,*)' has not yet been defined. ' 143 call inputerror(inpc,ipoinpc,iline, 144 & "*INITIAL CONDITIONS%",ier) 145 return 146 endif 147 do j=istartset(ii),iendset(ii) 148 if(ialset(j).gt.0) then 149 t0(ialset(j))=temperature 150 t1(ialset(j))=temperature 151 vold(0,ialset(j))=temperature 152 if((inoelfree.ne.0).or.(nuel_.gt.0)) then 153 t0g(1,ialset(j))=tempgrad1 154 t0g(2,ialset(j))=tempgrad2 155 t1g(1,ialset(j))=tempgrad1 156 t1g(2,ialset(j))=tempgrad2 157 endif 158 else 159 k=ialset(j-2) 160 do 161 k=k-ialset(j) 162 if(k.ge.ialset(j-1)) exit 163 t0(k)=temperature 164 t1(k)=temperature 165 vold(0,k)=temperature 166 if((inoelfree.ne.0).or.(nuel_.gt.0)) then 167 t0g(1,k)=tempgrad1 168 t0g(2,k)=tempgrad2 169 t1g(1,k)=tempgrad1 170 t1g(2,k)=tempgrad2 171 endif 172 enddo 173 endif 174 enddo 175 endif 176 enddo 177 return 178 elseif(textpart(ij)(1:11).eq.'TYPE=STRESS') then 179! 180 iprestr=1 181 do jj=1,n 182 if(textpart(jj)(1:4).eq.'USER') then 183! 184! residual stresses are defined by user subroutine 185! sigini 186! 187 iflag=1 188 ntens=6 189 ncrds=3 190 lrebar=0 191 do i=1,ne 192 indexe=ipkon(i) 193 if(lakon(i)(4:4).eq.'2') then 194 nope=20 195 elseif(lakon(i)(4:4).eq.'8') then 196 nope=8 197 elseif(lakon(i)(4:5).eq.'10') then 198 nope=10 199 elseif(lakon(i)(4:4).eq.'4') then 200 nope=4 201 elseif(lakon(i)(4:5).eq.'15') then 202 nope=15 203 elseif(lakon(i)(4:4).eq.'6') then 204 nope=6 205 else 206 cycle 207 endif 208! 209 if(lakon(i)(4:5).eq.'8R') then 210 mint3d=1 211 elseif(lakon(i)(4:7).eq.'20RB') then 212 if((lakon(i)(8:8).eq.'R').or. 213 & (lakon(i)(8:8).eq.'C')) then 214 mint3d=50 215 else 216 call beamintscheme(lakon(i),mint3d, 217 & ielprop(i),prop, 218 & null,xi,et,ze,weight) 219 endif 220 elseif((lakon(i)(4:4).eq.'8').or. 221 & (lakon(i)(4:6).eq.'20R')) then 222 mint3d=8 223 elseif(lakon(i)(4:4).eq.'2') then 224 mint3d=27 225 elseif(lakon(i)(4:5).eq.'10') then 226 mint3d=4 227 elseif(lakon(i)(4:4).eq.'4') then 228 mint3d=1 229 elseif(lakon(i)(4:5).eq.'15') then 230 mint3d=9 231 elseif(lakon(i)(4:4).eq.'6') then 232 mint3d=2 233 endif 234! 235 do j=1,nope 236 konl(j)=kon(indexe+j) 237 do k=1,3 238 xl(k,j)=co(k,konl(j)) 239 enddo 240 enddo 241! 242 do j=1,mint3d 243 if(lakon(i)(4:5).eq.'8R') then 244 xi=gauss3d1(1,j) 245 et=gauss3d1(2,j) 246 ze=gauss3d1(3,j) 247 weight=weight3d1(j) 248 elseif(lakon(i)(4:7).eq.'20RB') then 249 if((lakon(i)(8:8).eq.'R').or. 250 & (lakon(i)(8:8).eq.'C')) then 251 xi=gauss3d13(1,j) 252 et=gauss3d13(2,j) 253 ze=gauss3d13(3,j) 254 weight=weight3d13(j) 255 else 256 call beamintscheme(lakon(i),mint3d, 257 & ielprop(i),prop, 258 & j,xi,et,ze,weight) 259 endif 260 elseif((lakon(i)(4:4).eq.'8').or. 261 & (lakon(i)(4:6).eq.'20R')) 262 & then 263 xi=gauss3d2(1,j) 264 et=gauss3d2(2,j) 265 ze=gauss3d2(3,j) 266 weight=weight3d2(j) 267 elseif(lakon(i)(4:4).eq.'2') then 268 xi=gauss3d3(1,j) 269 et=gauss3d3(2,j) 270 ze=gauss3d3(3,j) 271 weight=weight3d3(j) 272 elseif(lakon(i)(4:5).eq.'10') then 273 xi=gauss3d5(1,j) 274 et=gauss3d5(2,j) 275 ze=gauss3d5(3,j) 276 weight=weight3d5(j) 277 elseif(lakon(i)(4:4).eq.'4') then 278 xi=gauss3d4(1,j) 279 et=gauss3d4(2,j) 280 ze=gauss3d4(3,j) 281 weight=weight3d4(j) 282 elseif(lakon(i)(4:5).eq.'15') then 283 xi=gauss3d8(1,j) 284 et=gauss3d8(2,j) 285 ze=gauss3d8(3,j) 286 weight=weight3d8(j) 287 elseif(lakon(i)(4:4).eq.'6') then 288 xi=gauss3d7(1,j) 289 et=gauss3d7(2,j) 290 ze=gauss3d7(3,j) 291 weight=weight3d7(j) 292 endif 293! 294 if(nope.eq.20) then 295 call shape20h(xi,et,ze,xl,xsj,shp,iflag) 296 elseif(nope.eq.8) then 297 call shape8h(xi,et,ze,xl,xsj,shp,iflag) 298 elseif(nope.eq.10) then 299 call shape10tet(xi,et,ze,xl,xsj,shp,iflag) 300 elseif(nope.eq.4) then 301 call shape4tet(xi,et,ze,xl,xsj,shp,iflag) 302 elseif(nope.eq.15) then 303 call shape15w(xi,et,ze,xl,xsj,shp,iflag) 304 else 305 call shape6w(xi,et,ze,xl,xsj,shp,iflag) 306 endif 307! 308 do k=1,3 309 pgauss(k)=0.d0 310 do i1=1,nope 311 pgauss(k)=pgauss(k)+ 312 & shp(4,i1)*co(k,konl(i1)) 313 enddo 314 enddo 315! 316 call sigini(prestr(1,j,i),pgauss,ntens,ncrds, 317 & i,j,layer,kspt,lrebar,rebarn) 318! 319 enddo 320 enddo 321 call getnewline(inpc,textpart,istat,n,key,iline,ipol, 322 & inl,ipoinp,inp,ipoinpc) 323 return 324 endif 325 enddo 326! 327! residual stresses are written explicitly in the input deck 328! 329 do 330 call getnewline(inpc,textpart,istat,n,key,iline,ipol,inl, 331 & ipoinp,inp,ipoinpc) 332 if((istat.lt.0).or.(key.eq.1)) return 333 do j=1,6 334 read(textpart(j+2)(1:20),'(f20.0)',iostat=istat) 335 & beta(j) 336 if(istat.gt.0) then 337 call inputerror(inpc,ipoinpc,iline, 338 & "*INITIAL CONDITIONS%",ier) 339 return 340 endif 341 enddo 342 read(textpart(1)(1:10),'(i10)',iostat=istat) l 343 if(istat.ne.0) then 344 call inputerror(inpc,ipoinpc,iline, 345 & "*INITIAL CONDITIONS%",ier) 346 return 347 endif 348 if(l.gt.ne) then 349 write(*,*) 350 & '*WARNING reading *INITIAL CONDITIONS: element ',l 351 write(*,*)' exceeds the largest defined ', 352 & 'element number' 353 cycle 354 endif 355 read(textpart(2)(1:10),'(i10)',iostat=istat) k 356 if(istat.eq.0) then 357 do j=1,6 358 prestr(j,k,l)=beta(j) 359 enddo 360 else 361 call inputerror(inpc,ipoinpc,iline, 362 & "*INITIAL CONDITIONS%",ier) 363 return 364 endif 365 enddo 366 return 367 elseif(textpart(ij)(1:18).eq.'TYPE=PLASTICSTRAIN') then 368! 369 iprestr=2 370 do 371 call getnewline(inpc,textpart,istat,n,key,iline,ipol,inl, 372 & ipoinp,inp,ipoinpc) 373 if((istat.lt.0).or.(key.eq.1)) return 374 do j=1,6 375 read(textpart(j+2)(1:20),'(f20.0)',iostat=istat) 376 & beta(j) 377 if(istat.gt.0) then 378 call inputerror(inpc,ipoinpc,iline, 379 & "*INITIAL CONDITIONS%",ier) 380 return 381 endif 382 enddo 383 read(textpart(1)(1:10),'(i10)',iostat=istat) l 384 if(istat.ne.0) then 385 call inputerror(inpc,ipoinpc,iline, 386 & "*INITIAL CONDITIONS%",ier) 387 return 388 endif 389 if(l.gt.ne) then 390 write(*,*) 391 & '*WARNING reading *INITIAL CONDITIONS: element ',l 392 write(*,*)' exceeds the largest defined ', 393 & 'element number' 394 cycle 395 endif 396 read(textpart(2)(1:10),'(i10)',iostat=istat) k 397 if(istat.eq.0) then 398 do j=1,6 399 prestr(j,k,l)=beta(j) 400 enddo 401 else 402 call inputerror(inpc,ipoinpc,iline, 403 & "*INITIAL CONDITIONS%",ier) 404 return 405 endif 406 enddo 407 return 408 elseif((textpart(ij)(1:17).eq.'TYPE=DISPLACEMENT').or. 409 & (textpart(ij)(1:18).eq.'TYPE=FLUIDVELOCITY')) then 410! 411 do 412 call getnewline(inpc,textpart,istat,n,key,iline,ipol,inl, 413 & ipoinp,inp,ipoinpc) 414 if((istat.lt.0).or.(key.eq.1)) return 415 read(textpart(2)(1:10),'(i10)',iostat=istat) idir 416 if(istat.gt.0) then 417 call inputerror(inpc,ipoinpc,iline, 418 & "*INITIAL CONDITIONS%",ier) 419 return 420 endif 421 read(textpart(3)(1:20),'(f20.0)',iostat=istat) dispvelo 422 if(istat.gt.0) then 423 call inputerror(inpc,ipoinpc,iline, 424 & "*INITIAL CONDITIONS%",ier) 425 return 426 endif 427 read(textpart(1)(1:10),'(i10)',iostat=istat) l 428 if(istat.eq.0) then 429 if(l.gt.nk) then 430 write(*,*) 431 & '*WARNING reading *INITIAL CONDITIONS: node ',l 432 write(*,*)' exceeds the largest defined ', 433 & 'node number' 434 cycle 435 endif 436 vold(idir,l)=dispvelo 437 else 438 read(textpart(1)(1:80),'(a80)',iostat=istat) noset 439 noset(81:81)=' ' 440 ipos=index(noset,' ') 441 noset(ipos:ipos)='N' 442c do ii=1,nset 443c if(set(ii).eq.noset) exit 444c enddo 445 call cident81(set,noset,nset,id) 446 ii=nset+1 447 if(id.gt.0) then 448 if(noset.eq.set(id)) then 449 ii=id 450 endif 451 endif 452 if(ii.gt.nset) then 453 noset(ipos:ipos)=' ' 454 write(*,*) 455 & '*ERROR reading *INITIAL CONDITIONS: node set ' 456 & ,noset 457 write(*,*)' has not yet been defined. ' 458 call inputerror(inpc,ipoinpc,iline, 459 & "*INITIAL CONDITIONS%",ier) 460 return 461 endif 462 do j=istartset(ii),iendset(ii) 463 if(ialset(j).gt.0) then 464 vold(idir,ialset(j))=dispvelo 465 else 466 k=ialset(j-2) 467 do 468 k=k-ialset(j) 469 if(k.ge.ialset(j-1)) exit 470 vold(idir,k)=dispvelo 471 enddo 472 endif 473 enddo 474 endif 475 enddo 476 return 477 elseif(textpart(ij)(1:13).eq.'TYPE=VELOCITY') then 478! 479 do 480 call getnewline(inpc,textpart,istat,n,key,iline,ipol,inl, 481 & ipoinp,inp,ipoinpc) 482 if((istat.lt.0).or.(key.eq.1)) return 483 read(textpart(2)(1:10),'(i10)',iostat=istat) idir 484 if(istat.gt.0) then 485 call inputerror(inpc,ipoinpc,iline, 486 & "*INITIAL CONDITIONS%",ier) 487 return 488 endif 489 read(textpart(3)(1:20),'(f20.0)',iostat=istat) velocity 490 if(istat.gt.0) then 491 call inputerror(inpc,ipoinpc,iline, 492 & "*INITIAL CONDITIONS%",ier) 493 return 494 endif 495 read(textpart(1)(1:10),'(i10)',iostat=istat) l 496 if(istat.eq.0) then 497 if(l.gt.nk) then 498 write(*,*) 499 & '*WARNING reading *INITIAL CONDITIONS: node ',l 500 write(*,*)' exceeds the largest defined ', 501 & 'node number' 502 cycle 503 endif 504 veold(idir,l)=velocity 505 else 506 read(textpart(1)(1:80),'(a80)',iostat=istat) noset 507 noset(81:81)=' ' 508 ipos=index(noset,' ') 509 noset(ipos:ipos)='N' 510c do ii=1,nset 511c if(set(ii).eq.noset) exit 512c enddo 513 call cident81(set,noset,nset,id) 514 ii=nset+1 515 if(id.gt.0) then 516 if(noset.eq.set(id)) then 517 ii=id 518 endif 519 endif 520 if(ii.gt.nset) then 521 noset(ipos:ipos)=' ' 522 write(*,*) 523 & '*ERROR reading *INITIAL CONDITIONS: node set ' 524 & ,noset 525 write(*,*)' has not yet been defined. ' 526 call inputerror(inpc,ipoinpc,iline, 527 & "*INITIAL CONDITIONS%",ier) 528 return 529 endif 530 do j=istartset(ii),iendset(ii) 531 if(ialset(j).gt.0) then 532 veold(idir,ialset(j))=velocity 533 else 534 k=ialset(j-2) 535 do 536 k=k-ialset(j) 537 if(k.ge.ialset(j-1)) exit 538 veold(idir,k)=velocity 539 enddo 540 endif 541 enddo 542 endif 543 enddo 544 return 545 elseif(textpart(ij)(1:13).eq.'TYPE=PRESSURE') then 546! 547 do 548 call getnewline(inpc,textpart,istat,n,key,iline,ipol,inl, 549 & ipoinp,inp,ipoinpc) 550 if((istat.lt.0).or.(key.eq.1)) return 551 read(textpart(2)(1:20),'(f20.0)',iostat=istat) pressure 552 if(istat.gt.0) then 553 call inputerror(inpc,ipoinpc,iline, 554 & "*INITIAL CONDITIONS%",ier) 555 return 556 endif 557 read(textpart(1)(1:10),'(i10)',iostat=istat) l 558 if(istat.eq.0) then 559 if(l.gt.nk) then 560 write(*,*) 561 & '*WARNING reading *INITIAL CONDITIONS: node ',l 562 write(*,*)' exceeds the largest defined ', 563 & 'node number' 564 cycle 565 endif 566 vold(4,l)=pressure 567 else 568 read(textpart(1)(1:80),'(a80)',iostat=istat) noset 569 noset(81:81)=' ' 570 ipos=index(noset,' ') 571 noset(ipos:ipos)='N' 572c do ii=1,nset 573c if(set(ii).eq.noset) exit 574c enddo 575 call cident81(set,noset,nset,id) 576 ii=nset+1 577 if(id.gt.0) then 578 if(noset.eq.set(id)) then 579 ii=id 580 endif 581 endif 582 if(ii.gt.nset) then 583 noset(ipos:ipos)=' ' 584 write(*,*) 585 & '*ERROR reading *INITIAL CONDITIONS: node set ' 586 & ,noset 587 write(*,*)' has not yet been defined. ' 588 call inputerror(inpc,ipoinpc,iline, 589 & "*INITIAL CONDITIONS%",ier) 590 return 591 endif 592 do j=istartset(ii),iendset(ii) 593 if(ialset(j).gt.0) then 594 vold(4,ialset(j))=pressure 595 else 596 k=ialset(j-2) 597 do 598 k=k-ialset(j) 599 if(k.ge.ialset(j-1)) exit 600 vold(4,k)=pressure 601 enddo 602 endif 603 enddo 604 endif 605 enddo 606 return 607 elseif(textpart(ij)(1:15).eq.'TYPE=TURBULENCE') then 608! 609 do 610 call getnewline(inpc,textpart,istat,n,key,iline,ipol,inl, 611 & ipoinp,inp,ipoinpc) 612 if((istat.lt.0).or.(key.eq.1)) return 613 do m=5,mi(2) 614 read(textpart(m-3)(1:20),'(f20.0)',iostat=istat) 615 & turbini(m) 616 enddo 617 if(istat.gt.0) then 618 call inputerror(inpc,ipoinpc,iline, 619 & "*INITIAL CONDITIONS%",ier) 620 return 621 endif 622 read(textpart(1)(1:10),'(i10)',iostat=istat) l 623 if(istat.eq.0) then 624 if(l.gt.nk) then 625 write(*,*) 626 & '*WARNING reading *INITIAL CONDITIONS: node ',l 627 write(*,*)' exceeds the largest defined ', 628 & 'node number' 629 cycle 630 endif 631 do m=5,mi(2) 632 vold(m,l)=turbini(m) 633 enddo 634 else 635 read(textpart(1)(1:80),'(a80)',iostat=istat) noset 636 noset(81:81)=' ' 637 ipos=index(noset,' ') 638 noset(ipos:ipos)='N' 639c do ii=1,nset 640c if(set(ii).eq.noset) exit 641c enddo 642 call cident81(set,noset,nset,id) 643 ii=nset+1 644 if(id.gt.0) then 645 if(noset.eq.set(id)) then 646 ii=id 647 endif 648 endif 649 if(ii.gt.nset) then 650 noset(ipos:ipos)=' ' 651 write(*,*) 652 & '*ERROR reading *INITIAL CONDITIONS: node set ' 653 & ,noset 654 write(*,*)' has not yet been defined. ' 655 call inputerror(inpc,ipoinpc,iline, 656 & "*INITIAL CONDITIONS%",ier) 657 return 658 endif 659 do j=istartset(ii),iendset(ii) 660 if(ialset(j).gt.0) then 661 do m=5,mi(2) 662 vold(m,ialset(j))=turbini(m) 663 enddo 664 else 665 k=ialset(j-2) 666 do 667 k=k-ialset(j) 668 if(k.ge.ialset(j-1)) exit 669 do m=5,mi(2) 670 vold(m,k)=turbini(m) 671 enddo 672 enddo 673 endif 674 enddo 675 endif 676 enddo 677 return 678 elseif(textpart(ij)(1:18).eq.'TYPE=TOTALPRESSURE') then 679! 680 do 681 call getnewline(inpc,textpart,istat,n,key,iline,ipol,inl, 682 & ipoinp,inp,ipoinpc) 683 if((istat.lt.0).or.(key.eq.1)) return 684 read(textpart(2)(1:20),'(f20.0)',iostat=istat) totpres 685 if(istat.gt.0) then 686 call inputerror(inpc,ipoinpc,iline, 687 & "*INITIAL CONDITIONS%",ier) 688 return 689 endif 690 read(textpart(1)(1:10),'(i10)',iostat=istat) l 691 if(istat.eq.0) then 692 if(l.gt.nk) then 693 write(*,*) 694 & '*WARNING reading *INITIAL CONDITIONS: node ',l 695 write(*,*)' exceeds the largest defined ', 696 & 'node number' 697 cycle 698 endif 699 vold(2,l)=totpres 700 else 701 read(textpart(1)(1:80),'(a80)',iostat=istat) noset 702 noset(81:81)=' ' 703 ipos=index(noset,' ') 704 noset(ipos:ipos)='N' 705c do ii=1,nset 706c if(set(ii).eq.noset) exit 707c enddo 708 call cident81(set,noset,nset,id) 709 ii=nset+1 710 if(id.gt.0) then 711 if(noset.eq.set(id)) then 712 ii=id 713 endif 714 endif 715 if(ii.gt.nset) then 716 noset(ipos:ipos)=' ' 717 write(*,*) 718 & '*ERROR reading *INITIAL CONDITIONS: node set ' 719 & ,noset 720 write(*,*)' has not yet been defined. ' 721 call inputerror(inpc,ipoinpc,iline, 722 & "*INITIAL CONDITIONS%",ier) 723 return 724 endif 725 do j=istartset(ii),iendset(ii) 726 if(ialset(j).gt.0) then 727 vold(2,ialset(j))=totpres 728 else 729 k=ialset(j-2) 730 do 731 k=k-ialset(j) 732 if(k.ge.ialset(j-1)) exit 733 vold(2,k)=totpres 734 enddo 735 endif 736 enddo 737 endif 738 enddo 739 return 740 elseif(textpart(ij)(1:13).eq.'TYPE=MASSFLOW') then 741! 742 do 743 call getnewline(inpc,textpart,istat,n,key,iline,ipol,inl, 744 & ipoinp,inp,ipoinpc) 745 if((istat.lt.0).or.(key.eq.1)) return 746 read(textpart(2)(1:20),'(f20.0)',iostat=istat) xmassflow 747 if(iaxial.eq.180) xmassflow=xmassflow/iaxial 748 if(istat.gt.0) then 749 call inputerror(inpc,ipoinpc,iline, 750 & "*INITIAL CONDITIONS%",ier) 751 return 752 endif 753 read(textpart(1)(1:10),'(i10)',iostat=istat) l 754 if(istat.eq.0) then 755 if(l.gt.nk) then 756 write(*,*) 757 & '*WARNING reading *INITIAL CONDITIONS: node ',l 758 write(*,*)' exceeds the largest defined ', 759 & 'node number' 760 cycle 761 endif 762 vold(1,l)=xmassflow 763 else 764 read(textpart(1)(1:80),'(a80)',iostat=istat) noset 765 noset(81:81)=' ' 766 ipos=index(noset,' ') 767 noset(ipos:ipos)='N' 768c do ii=1,nset 769c if(set(ii).eq.noset) exit 770c enddo 771 call cident81(set,noset,nset,id) 772 ii=nset+1 773 if(id.gt.0) then 774 if(noset.eq.set(id)) then 775 ii=id 776 endif 777 endif 778 if(ii.gt.nset) then 779 noset(ipos:ipos)=' ' 780 write(*,*) 781 & '*ERROR reading *INITIAL CONDITIONS: node set ' 782 & ,noset 783 write(*,*)' has not yet been defined. ' 784 call inputerror(inpc,ipoinpc,iline, 785 & "*INITIAL CONDITIONS%",ier) 786 return 787 endif 788 do j=istartset(ii),iendset(ii) 789 if(ialset(j).gt.0) then 790 vold(1,ialset(j))=xmassflow 791 else 792 k=ialset(j-2) 793 do 794 k=k-ialset(j) 795 if(k.ge.ialset(j-1)) exit 796 vold(1,k)=xmassflow 797 enddo 798 endif 799 enddo 800 endif 801 enddo 802 return 803! 804 elseif(textpart(ij)(1:13).eq.'TYPE=SOLUTION') then 805 user=.false. 806 do j=2,n 807 if(textpart(j)(1:4).eq.'USER') user=.true. 808 enddo 809c if(.not.user) then 810c write(*,*) 811c & '*ERROR reading *INITIAL CONDITIONS: TYPE=SOLUTION' 812c write(*,*) ' can only be used in combination with' 813c write(*,*) ' USER' 814c ier=1 815c return 816c endif 817 if(user) then 818! 819! internal state variables are read in file sdvini.f 820! 821 iflag=1 822 ncrds=3 823 do i=1,ne 824 indexe=ipkon(i) 825 if(lakon(i)(4:4).eq.'2') then 826 nope=20 827 elseif(lakon(i)(4:4).eq.'8') then 828 nope=8 829 elseif(lakon(i)(4:5).eq.'10') then 830 nope=10 831 elseif(lakon(i)(4:4).eq.'4') then 832 nope=4 833 elseif(lakon(i)(4:5).eq.'15') then 834 nope=15 835 elseif(lakon(i)(4:4).eq.'6') then 836 nope=6 837 else 838 cycle 839 endif 840! 841 if(lakon(i)(4:5).eq.'8R') then 842 mint3d=1 843 elseif((lakon(i)(4:4).eq.'8').or. 844 & (lakon(i)(4:6).eq.'20R')) then 845 mint3d=8 846 elseif(lakon(i)(4:4).eq.'2') then 847 mint3d=27 848 elseif(lakon(i)(4:5).eq.'10') then 849 mint3d=4 850 elseif(lakon(i)(4:4).eq.'4') then 851 mint3d=1 852 elseif(lakon(i)(4:5).eq.'15') then 853 mint3d=9 854 elseif(lakon(i)(4:4).eq.'6') then 855 mint3d=2 856 endif 857! 858 do j=1,nope 859 konl(j)=kon(indexe+j) 860 do k=1,3 861 xl(k,j)=co(k,konl(j)) 862 enddo 863 enddo 864! 865 do j=1,mint3d 866 if(lakon(i)(4:5).eq.'8R') then 867 xi=gauss3d1(1,j) 868 et=gauss3d1(2,j) 869 ze=gauss3d1(3,j) 870 weight=weight3d1(j) 871 elseif((lakon(i)(4:4).eq.'8').or. 872 & (lakon(i)(4:6).eq.'20R')) 873 & then 874 xi=gauss3d2(1,j) 875 et=gauss3d2(2,j) 876 ze=gauss3d2(3,j) 877 weight=weight3d2(j) 878 elseif(lakon(i)(4:4).eq.'2') then 879 xi=gauss3d3(1,j) 880 et=gauss3d3(2,j) 881 ze=gauss3d3(3,j) 882 weight=weight3d3(j) 883 elseif(lakon(i)(4:5).eq.'10') then 884 xi=gauss3d5(1,j) 885 et=gauss3d5(2,j) 886 ze=gauss3d5(3,j) 887 weight=weight3d5(j) 888 elseif(lakon(i)(4:4).eq.'4') then 889 xi=gauss3d4(1,j) 890 et=gauss3d4(2,j) 891 ze=gauss3d4(3,j) 892 weight=weight3d4(j) 893 elseif(lakon(i)(4:5).eq.'15') then 894 xi=gauss3d8(1,j) 895 et=gauss3d8(2,j) 896 ze=gauss3d8(3,j) 897 weight=weight3d8(j) 898 elseif(lakon(i)(4:4).eq.'6') then 899 xi=gauss3d7(1,j) 900 et=gauss3d7(2,j) 901 ze=gauss3d7(3,j) 902 weight=weight3d7(j) 903 endif 904! 905 if(nope.eq.20) then 906 call shape20h(xi,et,ze,xl,xsj,shp,iflag) 907 elseif(nope.eq.8) then 908 call shape8h(xi,et,ze,xl,xsj,shp,iflag) 909 elseif(nope.eq.10) then 910 call shape10tet(xi,et,ze,xl,xsj,shp,iflag) 911 elseif(nope.eq.4) then 912 call shape4tet(xi,et,ze,xl,xsj,shp,iflag) 913 elseif(nope.eq.15) then 914 call shape15w(xi,et,ze,xl,xsj,shp,iflag) 915 else 916 call shape6w(xi,et,ze,xl,xsj,shp,iflag) 917 endif 918! 919 do k=1,3 920 pgauss(k)=0.d0 921 do i1=1,nope 922 pgauss(k)=pgauss(k)+ 923 & shp(4,i1)*co(k,konl(i1)) 924 enddo 925 enddo 926! 927 call sdvini(xstate(1,j,i),pgauss,nstate_,ncrds, 928 & i,j,layer,kspt) 929! 930 enddo 931 enddo 932 call getnewline(inpc,textpart,istat,n,key,iline,ipol, 933 & inl,ipoinp,inp,ipoinpc) 934 return 935 else 936! 937! internal variables are written explicitly in the input deck 938! 939 do 940 call getnewline(inpc,textpart,istat,n,key,iline,ipol, 941 & inl,ipoinp,inp,ipoinpc) 942 if((istat.lt.0).or.(key.eq.1)) return 943! 944 if(nstate_.lt.6) then 945 ntot=nstate_ 946 else 947 ntot=6 948 endif 949! 950 do j=1,ntot 951 read(textpart(j+2)(1:20),'(f20.0)',iostat=istat) 952 & beta(j) 953 if(istat.gt.0) then 954 call inputerror(inpc,ipoinpc,iline, 955 & "*INITIAL CONDITIONS%",ier) 956 return 957 endif 958 enddo 959 read(textpart(1)(1:10),'(i10)',iostat=istat) l 960 if(istat.ne.0) then 961 call inputerror(inpc,ipoinpc,iline, 962 & "*INITIAL CONDITIONS%",ier) 963 return 964 endif 965 if(l.gt.ne) then 966 write(*,*) 967 & '*WARNING reading *INITIAL CONDITIONS: element ',l 968 write(*,*)' exceeds the largest defined ', 969 & 'element number' 970 cycle 971 endif 972 read(textpart(2)(1:10),'(i10)',iostat=istat) k 973 if(istat.eq.0) then 974 do j=1,ntot 975 xstate(j,k,l)=beta(j) 976 enddo 977 else 978 call inputerror(inpc,ipoinpc,iline, 979 & "*INITIAL CONDITIONS%",ier) 980 return 981 endif 982! 983 if(nstate_.gt.6) then 984 numberoflines=(nstate_-7)/8+1 985 do ii=1,numberoflines 986 if(ii.lt.numberoflines) then 987 jmax=8 988 else 989 jmax=nstate_-ntot 990 endif 991 call getnewline(inpc,textpart,istat,n,key,iline, 992 & ipol,inl,ipoinp,inp,ipoinpc) 993 if((istat.lt.0).or.(key.eq.1)) return 994 do j=1,jmax 995 read(textpart(j+2)(1:20),'(f20.0)', 996 & iostat=istat) beta(j) 997 if(istat.gt.0) 998 & call inputerror(inpc,ipoinpc,iline, 999 & "*INITIAL CONDITIONS%",ier) 1000 return 1001 xstate(ntot+j,k,l)=beta(j) 1002 enddo 1003 ntot=ntot+jmax 1004 enddo 1005 endif 1006! 1007 enddo 1008 return 1009 endif 1010! 1011 else 1012 write(*,*) 1013 &'*WARNING reading *INITIAL CONDITIONS: parameter not recognized:' 1014 write(*,*) ' ', 1015 & textpart(ij)(1:index(textpart(ij),' ')-1) 1016 call inputwarning(inpc,ipoinpc,iline, 1017 & "*INITIAL CONDITIONS%") 1018 endif 1019 enddo 1020! 1021 write(*,*) '*ERROR reading *INITIAL CONDITIONS: unknown type' 1022 write(*,*) ' ' 1023 call inputerror(inpc,ipoinpc,iline, 1024 & "*INITIAL CONDITIONS%",ier) 1025! 1026 return 1027 end 1028 1029