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 envtemp(itg,ieg,ntg,ntr,sideload,nelemload, 20 & ipkon,kon,lakon,ielmat,ne,nload,kontri,ntri,nloadtr, 21 & nflow,ndirboun,nactdog,nodeboun,nacteq,nboun, 22 & ielprop,prop,nteq,v,network,physcon,shcon,ntmat_, 23 & co,vold,set,nshcon,rhcon,nrhcon,mi, 24 & nmpc,nodempc,ipompc,labmpc,ikboun,nasym,ttime,time,iaxial) 25! 26! determines the number of gas temperatures and radiation 27! temperatures 28! 29 implicit none 30! 31 logical identity,walltemp,temperaturebc,pressurebc,massflowbcall, 32 & pressurebcall 33! 34 character*8 lakon(*) 35 character*20 labmpc(*) 36 character*20 sideload(*) 37 character*81 set(*) 38! 39 integer itg(*),ntg,ntr,nelemload(2,*),ipkon(*),network,mi(*), 40 & kon(*),ielmat(mi(3),*),ne,i,j,k,l,index,id,node,nload, 41 & ifaceq(8,6),ider,nasym,indexe,iaxial,networkmpcs, 42 & ifacet(6,4),ifacew(8,5),kontri3(3,1),kontri4(3,2), 43 & kontri6(3,4),kontri8(3,6),kontri(4,*),ntri, 44 & konf(8),nloadtr(*),nelem,nope,nopes,ig,nflow,ieg(*), 45 & ndirboun(*),nactdog(0:3,*),nboun,nodeboun(*),ntmat_, 46 & idir,ntq,nteq,nacteq(0:3,*),node1,node2,nodem, 47 & ielprop(*),idirf(8),iflag,imat,numf,nrhcon(*),nshcon(*), 48 & nmpc,nodempc(3,*),ipompc(*),ikboun(*),iplausi 49! 50 real*8 prop(*),f,xflow,nodef(8),df(8),v(0:mi(2),*),g(3), 51 & cp,r,physcon(*),shcon(0:3,ntmat_,*),rho,ttime,time, 52 & co(3,*),dvi,vold(0:mi(2),*),rhcon(*) 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 kontri3 /1,2,3/ 70 data kontri4 /1,2,4,2,3,4/ 71 data kontri6 /1,4,6,4,5,6,4,2,5,6,5,3/ 72 data kontri8 /1,5,8,8,5,7,8,7,4,5,2,6,5,6,7,7,6,3/ 73! 74 ntg=0 75 ntr=0 76 ntri=0 77! 78 walltemp=.false. 79 temperaturebc=.false. 80 pressurebc=.false. 81 massflowbcall=.true. 82 pressurebcall=.true. 83! 84 networkmpcs=0 85! 86! ordering the gas temperature nodes and counting them 87! counting the radiation temperatures 88! 89 do i=1,nload 90 if(sideload(i)(3:4).eq.'FC') then 91 walltemp=.true. 92 call nident(itg,nelemload(2,i),ntg,id) 93 if(id.gt.0) then 94 if(itg(id).eq.nelemload(2,i)) then 95 nactdog(0,nelemload(2,i))=1 96 cycle 97 endif 98 endif 99 ntg=ntg+1 100 do j=ntg,id+2,-1 101 itg(j)=itg(j-1) 102 enddo 103 itg(id+1)=nelemload(2,i) 104 nactdog(0,nelemload(2,i))=1 105! 106 elseif(sideload(i)(3:4).eq.'NP') then 107 call nident(itg,nelemload(2,i),ntg,id) 108 if(id.gt.0) then 109 if(itg(id).eq.nelemload(2,i)) then 110 nactdog(2,nelemload(2,i))=1 111 cycle 112 endif 113 endif 114 ntg=ntg+1 115 do j=ntg,id+2,-1 116 itg(j)=itg(j-1) 117 enddo 118 itg(id+1)=nelemload(2,i) 119 nactdog(2,nelemload(2,i))=1 120! 121 elseif(sideload(i)(3:4).eq.'CR') then 122 ntr=ntr+1 123 nelem=nelemload(1,i) 124 read(sideload(i)(2:2),'(i1)') ig 125! 126! number of nodes in the face 127! 128 if(lakon(nelem)(4:4).eq.'2') then 129 nope=20 130 nopes=8 131 elseif(lakon(nelem)(4:4).eq.'8') then 132 nope=8 133 nopes=4 134 elseif(lakon(nelem)(4:5).eq.'10') then 135 nope=10 136 nopes=6 137 elseif(lakon(nelem)(4:4).eq.'4') then 138 nope=4 139 nopes=3 140 elseif(lakon(nelem)(4:4).eq.'6') then 141 nope=6 142 if(ig.le.2) then 143 nopes=3 144 else 145 nopes=4 146 endif 147 elseif(lakon(nelem)(4:5).eq.'15') then 148 nope=15 149 if(ig.le.2) then 150 nopes=6 151 else 152 nopes=8 153 endif 154 endif 155! 156! nodes in the face 157! 158 if((nope.eq.20).or.(nope.eq.8)) then 159 do k=1,nopes 160 konf(k)=kon(ipkon(nelem)+ifaceq(k,ig)) 161 enddo 162 elseif((nope.eq.10).or.(nope.eq.4)) then 163 do k=1,nopes 164 konf(k)=kon(ipkon(nelem)+ifacet(k,ig)) 165 enddo 166 else 167 do k=1,nopes 168 konf(k)=kon(ipkon(nelem)+ifacew(k,ig)) 169 enddo 170 endif 171! 172! triangulation of the face 173! 174 nloadtr(ntr)=i 175 if((lakon(nelem)(4:4).eq.'2').or. 176 & ((lakon(nelem)(4:5).eq.'15').and.(ig.gt.2))) then 177 do k=1,6 178 ntri=ntri+1 179 do l=1,3 180 kontri(l,ntri)=konf(kontri8(l,k)) 181 enddo 182 kontri(4,ntri)=ntr 183 enddo 184 elseif((lakon(nelem)(4:4).eq.'8').or. 185 & ((lakon(nelem)(4:4).eq.'6').and.(ig.gt.2))) then 186 do k=1,2 187 ntri=ntri+1 188 do l=1,3 189 kontri(l,ntri)=konf(kontri4(l,k)) 190 enddo 191 kontri(4,ntri)=ntr 192 enddo 193 elseif((lakon(nelem)(4:5).eq.'10').or. 194 & ((lakon(nelem)(4:5).eq.'15').and.(ig.le.2))) then 195 do k=1,4 196 ntri=ntri+1 197 do l=1,3 198 kontri(l,ntri)=konf(kontri6(l,k)) 199 enddo 200 kontri(4,ntri)=ntr 201 enddo 202 elseif((lakon(nelem)(4:4).eq.'4').or. 203 & ((lakon(nelem)(4:4).eq.'6').and.(ig.le.2))) then 204 do k=1,1 205 ntri=ntri+1 206 do l=1,3 207 kontri(l,ntri)=konf(kontri3(l,k)) 208 enddo 209 kontri(4,ntri)=ntr 210 enddo 211 endif 212 endif 213 enddo 214! 215! storing the gas elements in a dedicated array 216! 217 nflow=0 218! 219 do i=1,ne 220 if(lakon(i)(1:1).eq.'D') then 221 if((lakon(i)(2:2).ne.' ').and.(network.ne.1)) then 222 nflow=nflow+1 223 ieg(nflow)=i 224 else 225 nasym=1 226! 227! removing gas nodes belonging to 'D '-elements 228! in which a 'FC'-film condition was applied from the 229! itg vector 230! 231 indexe=ipkon(i) 232 do j=1,3,2 233 node=kon(indexe+j) 234 call nident(itg,node,ntg,id) 235 if(id.gt.0) then 236 if(itg(id).eq.node) then 237 ntg=ntg-1 238 do k=id,ntg 239 itg(k)=itg(k+1) 240 enddo 241 nactdog(0,node)=0 242 nactdog(2,node)=0 243 endif 244 endif 245 enddo 246! 247 endif 248 endif 249! 250! removing gas nodes belonging to advective elements 251! (last node in the element topology) 252! these are usually gas nodes not belonging to any 253! "D"-element 254! 255 if(lakon(i)(1:7).eq.'ESPRNGF') then 256 read(lakon(i)(8:8),'(i1)') nopes 257 nope=nopes+1 258 node=kon(ipkon(i)+nope) 259! 260 call nident(itg,node,ntg,id) 261 if(id.gt.0) then 262 if(itg(id).eq.node) then 263 ntg=ntg-1 264 do k=id,ntg 265 itg(k)=itg(k+1) 266 enddo 267 nactdog(0,node)=0 268 nactdog(2,node)=0 269 endif 270 endif 271 endif 272 enddo 273! 274! mass flux nodes are also taken as unknowns in the 275! gas temperature system; determining the active 276! degrees of freedom 277! 278! first node of the flow element 279! 280 do i=1,nflow 281 index=ipkon(ieg(i)) 282 node=kon(index+1) 283 if (node.eq.0) cycle 284 call nident(itg,node,ntg,id) 285 if(id.gt.0) then 286 if(itg(id).eq.node) then 287! 288! upstream depth of SO,WO and DO is known 289! 290 if((lakon(ieg(i))(4:7).eq.'CHSO').or. 291 & (lakon(ieg(i))(4:7).eq.'CHWO').or. 292 & (lakon(ieg(i))(4:7).eq.'CHDO')) cycle 293 nactdog(0,node)=1 294 nactdog(2,node)=1 295 cycle 296 endif 297 endif 298 ntg=ntg+1 299 do j=ntg,id+2,-1 300 itg(j)=itg(j-1) 301 enddo 302 itg(id+1)=node 303! 304! upstream depth of SO,WO and DO is known 305! 306 if((lakon(ieg(i))(4:7).eq.'CHSO').or. 307 & (lakon(ieg(i))(4:7).eq.'CHWO').or. 308 & (lakon(ieg(i))(4:7).eq.'CHDO')) cycle 309 nactdog(0,node)=1 310 nactdog(2,node)=1 311 enddo 312! 313! middle node of the flow element :flux 314! 315 do i=1,nflow 316 index=ipkon(ieg(i)) 317 node=kon(index+2) 318 call nident(itg,node,ntg,id) 319 if(id.gt.0) then 320 if(itg(id).eq.node) cycle 321 endif 322 ntg=ntg+1 323 do j=ntg,id+2,-1 324 itg(j)=itg(j-1) 325 enddo 326 itg(id+1)=node 327 nactdog(1,node)=1 328! 329! variable geometric property 330! 331 if(lakon(ieg(i))(6:7).eq.'GV') then 332 index=ielprop(ieg(i)) 333 if(prop(index+2).le.0.d0) nactdog(3,node)=1 334 elseif((lakon(ieg(i))(4:7).eq.'CHSG').or. 335 & (lakon(ieg(i))(4:7).eq.'CHWE').or. 336 & (lakon(ieg(i))(4:7).eq.'CHDS')) then 337 nactdog(3,node)=1 338 elseif(lakon(ieg(i))(2:7).eq.'ACCTUB') then 339! 340 index=ielprop(ieg(i)) 341 if(prop(index+1).eq.2) then 342! Interval factor unknown,setting a DOF for geometry 343 nactdog(3,node)=1 344 elseif(prop(index+1).eq.3) then 345! Hole diameter factor unknown,setting a DOF for geometry 346 nactdog(3,node)=1 347 endif 348 endif 349 enddo 350! 351! third node of the flow element 352! 353 do i=1,nflow 354 index=ipkon(ieg(i)) 355 node=kon(index+3) 356 if (node.eq.0) cycle 357 call nident(itg,node,ntg,id) 358 if(id.gt.0) then 359 if(itg(id).eq.node) then 360! 361! downstream depth of SG,WE and DS is known 362! 363 if((lakon(ieg(i))(4:7).eq.'CHSG').or. 364 & (lakon(ieg(i))(4:7).eq.'CHWE').or. 365 & (lakon(ieg(i))(4:7).eq.'CHDS')) cycle 366 nactdog(0,node)=1 367 nactdog(2,node)=1 368 cycle 369 endif 370 endif 371 ntg=ntg+1 372 do j=ntg,id+2,-1 373 itg(j)=itg(j-1) 374 enddo 375 itg(id+1)=node 376! 377! downstream depth of SG,WE and DS is known 378! 379 if((lakon(ieg(i))(4:7).eq.'CHSG').or. 380 & (lakon(ieg(i))(4:7).eq.'CHWE').or. 381 & (lakon(ieg(i))(4:7).eq.'CHDS')) cycle 382 nactdog(0,node)=1 383 nactdog(2,node)=1 384 enddo 385! 386! tagging the network MPC's 387! 388 do i=1,nmpc 389! 390! check whether network MPC 391! 392 index=ipompc(i) 393 do 394 node=nodempc(1,index) 395 call nident(itg,node,ntg,id) 396 if(id.gt.0) then 397 if(itg(id).eq.node) then 398 labmpc(i)(1:7)='NETWORK' 399 networkmpcs=1 400 exit 401 endif 402 endif 403 index=nodempc(3,index) 404 if(index.eq.0) exit 405 enddo 406 enddo 407! 408! taking the MPC network nodes into account 409! 410 if(networkmpcs.eq.1) then 411 do i=1,nmpc 412 if(labmpc(i)(1:7).ne.'NETWORK') cycle 413! 414 index=ipompc(i) 415 do 416 node=nodempc(1,index) 417 call nident(itg,node,ntg,id) 418 if(id.gt.0) then 419 if(itg(id).eq.node) then 420 nactdog(nodempc(2,index),node)=1 421 index=nodempc(3,index) 422 if(index.eq.0) then 423 exit 424 else 425 cycle 426 endif 427 endif 428 endif 429! 430! adding a node to itg 431! 432 ntg=ntg+1 433 do j=ntg,id+2,-1 434 itg(j)=itg(j-1) 435 enddo 436 itg(id+1)=node 437 nactdog(nodempc(2,index),node)=1 438 index=nodempc(3,index) 439 if(index.eq.0) exit 440 enddo 441 enddo 442 endif 443! 444! subtracting the SPC conditions 445! 446 do i=1,nboun 447 node=nodeboun(i) 448 call nident(itg,node,ntg,id) 449 if (id.gt.0) then 450 if (itg(id).eq.node) then 451 idir=ndirboun(i) 452 nactdog(idir,node)=0 453 if(idir.eq.0) then 454 temperaturebc=.true. 455 elseif(idir.eq.2) then 456 pressurebc=.true. 457 endif 458 endif 459 endif 460 enddo 461! 462! temporarily removing the dependent nodes of the MPC's 463! only for mass flow and pressure 464! 465! these are the only dofs for which the corresponding equations 466! (mass equilibrium in end node and element equation in middle 467! node) are not applied in the nodes where the dofs are lacking 468! (mass flow in middle node and pressure in end node) 469! 470 if(networkmpcs.eq.1) then 471 do i=1,nmpc 472 if(labmpc(i)(1:7).ne.'NETWORK') cycle 473 index=ipompc(i) 474 idir=nodempc(2,index) 475 if((idir.eq.1).or.(idir.eq.2)) then 476 node=nodempc(1,index) 477 call nident(itg,node,ntg,id) 478 if(id.gt.0) then 479 if(itg(id).eq.node) then 480 nactdog(idir,node)=0 481 endif 482 endif 483 endif 484 enddo 485 endif 486! 487! determining the active equations 488! 489! element contributions 490! 491 iflag=0 492 do i=1,nflow 493 nelem=ieg(i) 494 index=ipkon(nelem) 495 node1=kon(index+1) 496 nodem=kon(index+2) 497 node2=kon(index+3) 498! 499! "end of network" element ---X---O 500! 1 m 2 501! 502 if(node1.eq.0)then 503 if ((nactdog(1,nodem).ne.0).or.(nactdog(3,nodem).ne.0))then 504 nacteq(1,node2)=1 ! mass equation 505 endif 506 if (nactdog(0,node2).ne.0)then 507 nacteq(0,node2)=1 ! energy equation 508 endif 509! 510! "end of network" element node O---X--- 511! 1 m 2 512 elseif (node2.eq.0) then 513 if ((nactdog(1,nodem).ne.0).or.(nactdog(3,nodem).ne.0))then 514 nacteq(1,node1)=1 ! mass equation 515 endif 516 if (nactdog(0,node1).ne.0)then 517 nacteq(0,node1)=1 ! energy equation 518 endif 519! 520! "flow element" O---X---O 521! 1 m 2 522! 523 else 524 if((nactdog(1,nodem).ne.0).or.(nactdog(3,nodem).ne.0)) then 525 nacteq(1,node2)=1 ! mass equation 526 nacteq(1,node1)=1 ! mass equation 527 endif 528! 529 ider=0 530 call flux(node1,node2,nodem,nelem,lakon,kon,ipkon, 531 & nactdog,identity,ielprop,prop,iflag,v,xflow,f, 532 & nodef,idirf,df,cp,r,rho,physcon,g,co,dvi,numf, 533 & vold,set,shcon,nshcon,rhcon,nrhcon,ntmat_,mi,ider, 534 & ttime,time,iaxial,iplausi) 535! 536 if (.not.identity) then 537 nacteq(2,nodem)=1 ! momentum equation 538 endif 539! 540 if (nactdog(0,node1).ne.0)then 541 nacteq(0,node1)=1 ! energy equation 542 endif 543! 544 if (nactdog(0,node2).ne.0)then 545 nacteq(0,node2)=1 ! energy equation 546 endif 547 endif 548 enddo 549! 550! wall convective contributions 551! 552 do i=1, nload 553 if((sideload(i)(3:4)).eq.'FC')then 554 node=nelemload(2,i) 555 if (nactdog(0,node).ne.0) then 556 nacteq(0,node)=1 557 endif 558 endif 559 enddo 560! 561! restoring the dependent nodes of the MPC's 562! only for mass flow and pressure 563! 564 if(networkmpcs.eq.1) then 565 do i=1,nmpc 566 if(labmpc(i)(1:7).ne.'NETWORK') cycle 567 index=ipompc(i) 568 idir=nodempc(2,index) 569 if((idir.eq.1).or.(idir.eq.2)) then 570 node=nodempc(1,index) 571 call nident(itg,node,ntg,id) 572 if(id.gt.0) then 573 if(itg(id).eq.node) then 574 nactdog(idir,node)=1 575 endif 576 endif 577 endif 578 enddo 579 endif 580! 581! removing the energy equation from those end nodes for which 582! the temperature constitutes the first term in a network MPC 583! 584 if(networkmpcs.eq.1) then 585 do i=1,nmpc 586 if(labmpc(i)(1:7).ne.'NETWORK') cycle 587 index=ipompc(i) 588 idir=nodempc(2,index) 589 if(idir.eq.0) then 590 node=nodempc(1,index) 591 call nident(itg,node,ntg,id) 592 if(id.gt.0) then 593 if(itg(id).eq.node) then 594 nacteq(0,node)=0 595 endif 596 endif 597 endif 598 enddo 599 endif 600! 601! check whether all mass flow is known 602! 603 do i=1,ntg 604 node=itg(i) 605 if((nactdog(1,node).ne.0).or.(nactdog(3,node).ne.0)) then 606 massflowbcall=.false. 607 exit 608 endif 609 enddo 610! 611! check whether all pressures are known 612! 613 do i=1,ntg 614 node=itg(i) 615 if(nactdog(2,node).ne.0) then 616 pressurebcall=.false. 617 exit 618 endif 619 enddo 620! 621! check for special cases 622! 623 if(massflowbcall.and.((.not.pressurebc).or.(pressurebcall))) then 624! 625! purely thermal (only set to 1 if D-type elements are present) 626! 627 if(network.gt.1) network=2 628 do i=1,nflow 629 nelem=ieg(i) 630 index=ipkon(nelem) 631 node1=kon(index+1) 632 nodem=kon(index+2) 633 node2=kon(index+3) 634 nacteq(2,nodem)=0 635 if(node1.ne.0) nactdog(2,node1)=0 636 if(node2.ne.0) nactdog(2,node2)=0 637 enddo 638 elseif((.not.temperaturebc).and.(.not.walltemp)) then 639! 640! pure liquid dynamics 641! 642 write(*,*) '*INFO in envtemp: no thermal boundary conditions' 643 write(*,*) ' detected; the network is considered to be' 644 write(*,*) ' athermal and no gas temperatures will be' 645 write(*,*) ' calculated' 646 network=4 647 do i=1,ntg 648 node=itg(i) 649 nactdog(0,node)=0 650 nacteq(0,node)=0 651 enddo 652 elseif((.not.temperaturebc).and.walltemp) then 653 write(*,*) '*ERROR in envtemp: at least one temperature' 654 write(*,*) ' boundary condition must be given' 655 call exit(201) 656 elseif(.not.pressurebc) then 657 write(*,*) '*ERROR in envtemp: at least one pressure' 658 write(*,*) ' boundary condition must be given' 659 call exit(201) 660 endif 661! 662! check whether a specific gas constant was defined for all fluid 663! elements (except for purely thermal calculations) 664! 665 if(network.gt.2) then 666 do i=1,nflow 667 nelem=ieg(i) 668 if((lakon(nelem)(2:3).eq.'LI').or. 669 & (lakon(nelem)(2:3).eq.'LP').or. 670 & (lakon(nelem)(2:3).eq.' ')) cycle 671 imat=ielmat(1,nelem) 672 r=shcon(3,1,imat) 673 if(r.lt.1.d-10) then 674 write(*,*)'*ERROR in envtemp: specific gas', 675 & 'constant is close to zero' 676 call exit(201) 677 endif 678 enddo 679 endif 680! 681! check whether the temperature at each inlet or outlet node 682! is given 683! 684 if(network.lt.4) then 685 do i=1,nflow 686 nelem=ieg(i) 687 index=ipkon(nelem) 688 node1=kon(index+1) 689 node2=kon(index+3) 690 if(node1.eq.0) then 691 if(nactdog(0,node2).ne.0) then 692 write(*,*) '*WARNING in envtemp: it is advised to' 693 write(*,*) ' define the temperature at all' 694 write(*,*) ' inlets and outlets by a boundary' 695 write(*,*) ' condition. This is lacking for' 696 write(*,*) ' node ',node2,' of element ',nelem 697 endif 698 endif 699 if(node2.eq.0) then 700 if(nactdog(0,node1).ne.0) then 701 write(*,*) '*WARNING in envtemp: it is advised to' 702 write(*,*) ' define the temperature at all' 703 write(*,*) ' inlets and outlets by a boundary' 704 write(*,*) ' condition. This is lacking for' 705 write(*,*) ' node ',node1,' of element ',nelem 706 endif 707 endif 708 enddo 709 endif 710! 711! numbering the active equations 712! 713 nteq=0 714 do i=1,ntg 715 node=itg(i) 716 do j=0,2 717 if (nacteq(j,node).ne.0) then 718 nteq=nteq+1 719 nacteq(j,node)=nteq 720 endif 721 enddo 722! write(30,*) 'unknowns ',node,(nactdog(j,node),j=0,3) 723 enddo 724! do i=1,ntg 725! node=itg(i) 726! write(30,*) 'equations',node,(nacteq(j,node),j=0,2) 727! enddo 728! 729! taking network MPC's into account 730! 731 if(networkmpcs.eq.1) then 732 do i=1,nmpc 733 if(labmpc(i)(1:7).eq.'NETWORK') nteq=nteq+1 734 enddo 735 endif 736! 737! numbering the active degrees of freedom 738! 739 ntq=0 740 do i=1,ntg 741 node=itg(i) 742 do j=0,3 743 if (nactdog(j,node).ne.0) then 744 ntq=ntq+1 745 nactdog(j,node)=ntq 746 endif 747 enddo 748 enddo 749c 750c open(30,file='dummy',status='unknown') 751c write(30,*) 'nactdog' 752c do i=1,ntg 753c write(30,*) itg(i),(nactdog(j,itg(i)),j=0,3) 754c enddo 755c 756c write(30,*) '' 757c write(30,*) 'nacteq' 758c do i=1,ntg 759c write(30,*) itg(i),(nacteq(j,itg(i)),j=0,3) 760c enddo 761c close(30) 762! 763 if(ntq.ne.nteq) then 764 write(*,*) '*ERROR in envtemp:' 765 write(*,*) '*****number of network equations is not equal to' 766 write(*,*) ' number of active degrees of freedom*****' 767 write(*,*) ' # of network equations = ',nteq 768 write(*,*) ' # of active degrees of freedom= ',ntq 769 call exit(201) 770 endif 771! 772! for isothermal gas pipes the energy equation in the 773! topologically downstream node is replaced by an equation 774! expressing the equality of the static temperature at both 775! ends of the pipe. To this end these downstream nodes are 776! referring in nacteq(3,*) to the topologically upstream node 777! 778! if the temperature in the downstream node is a boundary 779! condition (i.e. the energy equation is not built), the 780! energy equation in the upstream node is replaced. In that 781! case nacteq(3,upstreamnode) refers to the downstream node. 782! 783! if both nodes are boundary conditions, nothing is done 784! 785 do i=1,nflow 786 nelem=ieg(i) 787 if((lakon(nelem)(1:4).eq."DGAP") 788 & .and.(lakon(nelem)(6:6).eq."I")) then 789! 790 index=ipkon(nelem) 791 node1=kon(index+1) 792 node2=kon(index+3) 793 if((node1.eq.0).or.(node2.eq.0)) cycle 794! 795 if(nacteq(0,node2).ne.0) then 796 nacteq(3,node2)=node1 797 elseif(nacteq(0,node1).ne.0) then 798 nacteq(3,node1)=node2 799 endif 800 endif 801 enddo 802! 803 return 804 end 805 806 807