CalculiX  2.13
A Free Software Three-Dimensional Structural Finite Element Program
Functions/Subroutines
resultnet.f File Reference

Go to the source code of this file.

Functions/Subroutines

subroutine resultnet (itg, ieg, ntg, bc, nload, sideload, nelemload, xloadact, lakon, ntmat_, v, shcon, nshcon, ipkon, kon, co, nflow, iinc, istep, dtime, ttime, time, ikforc, ilforc, xforcact, nforc, ielmat, nteq, prop, ielprop, nactdog, nacteq, iin, physcon, camt, camf, camp, rhcon, nrhcon, ipobody, ibody, xbodyact, nbody, dtheta, vold, xloadold, reltime, nmethod, set, mi, ineighe, cama, vamt, vamf, vamp, vama, nmpc, nodempc, ipompc, coefmpc, labmpc, iaxial, qat, qaf, ramt, ramf, ramp, cocon, ncocon, iponoel, inoel, iplausi)
 

Function/Subroutine Documentation

◆ resultnet()

subroutine resultnet ( integer, dimension(*)  itg,
integer, dimension(*)  ieg,
integer  ntg,
real*8, dimension(nteq)  bc,
integer  nload,
character*20, dimension(*)  sideload,
integer, dimension(2,*)  nelemload,
real*8, dimension(2,*)  xloadact,
character*8, dimension(*)  lakon,
integer  ntmat_,
real*8, dimension(0:mi(2),*)  v,
real*8, dimension(0:3,ntmat_,*)  shcon,
integer, dimension(*)  nshcon,
integer, dimension(*)  ipkon,
integer, dimension(*)  kon,
real*8, dimension(3,*)  co,
integer  nflow,
integer  iinc,
integer  istep,
real*8  dtime,
real*8  ttime,
real*8  time,
integer, dimension(*)  ikforc,
integer, dimension(*)  ilforc,
real*8, dimension(*)  xforcact,
integer  nforc,
integer, dimension(mi(3),*)  ielmat,
integer  nteq,
real*8, dimension(*)  prop,
integer, dimension(*)  ielprop,
integer, dimension(0:3,*)  nactdog,
integer, dimension(0:3,*)  nacteq,
integer  iin,
real*8, dimension(*)  physcon,
real*8, dimension(*)  camt,
real*8, dimension(*)  camf,
real*8, dimension(*)  camp,
real*8, dimension(0:1,ntmat_,*)  rhcon,
integer, dimension(*)  nrhcon,
integer, dimension(2,*)  ipobody,
integer, dimension(3,*)  ibody,
real*8, dimension(7,*)  xbodyact,
integer  nbody,
real*8  dtheta,
real*8, dimension(0:mi(2),*)  vold,
real*8, dimension(2,*)  xloadold,
real*8  reltime,
integer  nmethod,
character*81, dimension(*)  set,
integer, dimension(*)  mi,
integer, dimension(*)  ineighe,
real*8, dimension(*)  cama,
real*8  vamt,
real*8  vamf,
real*8  vamp,
real*8  vama,
integer  nmpc,
integer, dimension(3,*)  nodempc,
integer, dimension(*)  ipompc,
real*8, dimension(*)  coefmpc,
character*20, dimension(*)  labmpc,
integer  iaxial,
real*8  qat,
real*8  qaf,
real*8  ramt,
real*8  ramf,
real*8  ramp,
real*8, dimension(0:6,ntmat_,*)  cocon,
integer, dimension(2,*)  ncocon,
integer, dimension(*)  iponoel,
integer, dimension(2,*)  inoel,
integer  iplausi 
)
31 !
32  implicit none
33 !
34  logical identity
35  character*8 lakonl,lakon(*)
36  character*20 sideload(*),labmpc(*)
37  character*81 set(*)
38 !
39  integer mi(*),itg(*),ieg(*),ntg,nteq,nflow,nload,
40  & ielmat(mi(3),*),iflag,ider,iaxial,nalt,nalf,
41  & nelemload(2,*),nope,nopes,mint2d,i,j,k,l,nrhcon(*),
42  & node,imat,ntmat_,id,ifaceq(8,6),ifacet(6,4),numf,
43  & ifacew(8,5),node1,node2,nshcon(*),nelem,ig,index,konl(20),
44  & ipkon(*),kon(*),idof,ineighe(*),idir,ncocon(2,*),
45  & iinc,istep,jltyp,nfield,ikforc(*),ipobody(2,*),
46  & ilforc(*),nforc,nodem,idirf(8),ieq,nactdog(0:3,*),nbody,
47  & nacteq(0:3,*),ielprop(*),nodef(8),iin,kflag,ibody(3,*),icase,
48  & inv, index2,nmethod,nelem0,nodem0,nelem1,nodem1,nelem2,
49  & nodem2,nelemswirl,nmpc,nodempc(3,*),ipompc(*),
50  & iponoel(*),inoel(2,*),iplausi,indexe
51 !
52  real*8 bc(nteq),xloadact(2,*),cp,h(2),physcon(*),r,dvi,rho,
53  & xl2(3,8),coords(3),dxsj2,temp,xi,et,weight,xsj2(3),
54  & gastemp,v(0:mi(2),*),shcon(0:3,ntmat_,*),co(3,*),shp2(7,8),
55  & field,prop(*),tg1,tg2,dtime,ttime,time,g(3),eta,
56  & xforcact(*),areaj,xflow,tvar(2),f,df(8),camt(*),camf(*),
57  & camp(*),tl2(8),cama(*),vamt,vamf,vamp,vama,term,
58  & rhcon(0:1,ntmat_,*),xbodyact(7,*),sinktemp,kappa,a,t,tt,pt,
59  & dtheta,ts1,ts2,xs2(3,7),xk1,xk2,xdenom1,xdenom2,expon,pt1,
60  & pt2,dt1,dt2,xcst,xnum1,xnum2,qred_crit,xflow_crit,
61  & xflow0,xflow1,reltime,coefmpc(*),qat,qaf,ramt,ramf,ramp,
62  & xflow2,r1,r2,rout,rin,uout,uin,heat,cocon(0:6,ntmat_,*),
63  & cp_cor,u,ct,vold(0:mi(2),*),xloadold(2,*),omega,bnac,
64  & xflow360,tt1,tt2,t1,t2,heatnod,heatfac
65 !
66  include "gauss.f"
67 !
68  data ifaceq /4,3,2,1,11,10,9,12,
69  & 5,6,7,8,13,14,15,16,
70  & 1,2,6,5,9,18,13,17,
71  & 2,3,7,6,10,19,14,18,
72  & 3,4,8,7,11,20,15,19,
73  & 4,1,5,8,12,17,16,20/
74  data ifacet /1,3,2,7,6,5,
75  & 1,2,4,5,9,8,
76  & 2,3,4,6,10,9,
77  & 1,4,3,8,10,7/
78  data ifacew /1,3,2,9,8,7,0,0,
79  & 4,5,6,10,11,12,0,0,
80  & 1,2,5,4,7,14,10,13,
81  & 2,3,6,5,8,15,11,14,
82  & 4,6,3,1,12,15,9,13/
83  data iflag /2/
84 !
85  kflag=2
86  ider=0
87 !
88  tvar(1)=time
89  tvar(2)=ttime+time
90 !
91  heatnod=0.d0
92 !
93 ! calculating the maximum correction to the solution in the
94 ! present network iteration
95 !
96  camt(1)=0.d0
97  camf(1)=0.d0
98  camp(1)=0.d0
99  cama(1)=0.d0
100 !
101  camt(2)=0.5d0
102  camf(2)=0.5d0
103  camp(2)=0.5d0
104  cama(2)=0.5d0
105 !
106 ! maximum change in the solution since the start of the
107 ! network iterations (= entering radflowload)
108 !
109  vamt=0.d0
110  vamf=0.d0
111  vamp=0.d0
112  vama=0.d0
113 !
114  do i=1,ntg
115  node=itg(i)
116  do j=0,3
117  if(nactdog(j,node).eq.0) cycle
118  idof=nactdog(j,node)
119 !
120  if(dabs(v(j,node)).gt.1.d-30) then
121  if(dabs(bc(idof)/v(j,node)).lt.1.d-13) bc(idof)=0.d0
122  endif
123  bnac=bc(idof)
124 !
125  if(j.eq.0) then
126  if(dabs(bnac).gt.camt(1)) then
127  camt(1)=dabs(bc(idof))
128  camt(2)=node+0.5d0
129  endif
130  elseif(j.eq.1) then
131  if(dabs(bnac).gt.camf(1)) then
132  camf(1)=dabs(bc(idof))
133  camf(2)=node+0.5d0
134  endif
135  elseif(j.eq.2) then
136  if(dabs(bnac).gt.camp(1)) then
137  camp(1)=dabs(bc(idof))
138  camp(2)=node+0.5d0
139  endif
140  else
141  if(dabs(bnac).gt.cama(1)) then
142  cama(1)=dabs(bc(idof))
143  cama(2)=node+0.5d0
144  endif
145  endif
146  enddo
147  enddo
148 !
149 ! updating v
150 ! vold is the solution at the entry of radflowload (= at
151 ! the start of the network iterations)
152 !
153  do i=1,ntg
154  node=itg(i)
155  do j=0,2
156  if(nactdog(j,node).eq.0) cycle
157  v(j,node)=v(j,node)+bc(nactdog(j,node))*dtheta
158 !
159  if((j.eq.0).and.(dabs(v(j,node)-vold(j,node)).gt.vamt)) then
160  vamt=dabs(v(j,node)-vold(j,node))
161  elseif((j.eq.1).and.
162  & (dabs(v(j,node)-vold(j,node)).gt.vamf)) then
163  vamf=dabs(v(j,node)-vold(j,node))
164  elseif((j.eq.2).and.
165  & (dabs(v(j,node)-vold(j,node)).gt.vamp)) then
166  vamp=dabs(v(j,node)-vold(j,node))
167  endif
168 !
169  enddo
170  enddo
171 !
172 ! update geometry changes
173 !
174  do i=1,nflow
175  if(lakon(ieg(i))(6:7).eq.'GV') then
176  index=ipkon(ieg(i))
177  node=kon(index+2)
178  if(nactdog(3,node).eq.0) cycle
179  index=ielprop(ieg(i))
180  v(3,node)=v(3,node)+bc(nactdog(3,node))*dtheta
181  if(v(3,node).gt.0.99999d0) then
182  v(3,node)=0.99999d0
183  elseif(v(3,node).lt.0.12501) then
184  v(3,node)=0.12501d0
185  endif
186  if(dabs(v(3,node)).gt.vama) vama=dabs(v(3,node))
187 !
188 ! Geometry factor of an ACC tube
189 ! for all versions of the ACC tube
190 !
191  elseif(lakon(ieg(i))(2:7).eq.'ACCTUB') then
192  index=ipkon(ieg(i))
193  node=kon(index+2)
194  if(nactdog(3,node).eq.0) cycle
195  index=ielprop(ieg(i))
196 ! Interval factor unknown
197  if(nint(prop(index+1)).eq.2) then
198 ! Using a smaller step gets better convergence(factor 0.5)
199  v(3,node)=v(3,node)+bc(nactdog(3,node))*
200  & dtheta
201  if(v(3,node).lt.0.1)then
202  v(3,node) = 0.1
203  endif
204 ! Hole diameter factor unknown
205  elseif(nint(prop(index+1)).eq.3) then
206  v(3,node)=v(3,node)+bc(nactdog(3,node))*dtheta
207  if(v(3,node).lt.0.1)then
208  v(3,node) = 0.1
209  endif
210  endif
211  if(dabs(v(3,node)).gt.vama) vama=dabs(v(3,node))
212 !
213 ! update location of hydraulic jump
214 !
215  elseif(lakon(ieg(i))(2:5).eq.'LICH') then
216  if((lakon(ieg(i))(6:7).eq.'SG').or.
217  & (lakon(ieg(i))(6:7).eq.'WE').or.
218  & (lakon(ieg(i))(6:7).eq.'DS')) then
219  index=ipkon(ieg(i))
220  node=kon(index+2)
221  if(nactdog(3,node).eq.0) cycle
222  index=ielprop(ieg(i))
223  if(lakon(ieg(i))(6:7).eq.'SG') then
224  eta=prop(index+4)+bc(nactdog(3,node))*dtheta
225  prop(index+4)=eta
226  nelem=nint(prop(index+7))
227  elseif(lakon(ieg(i))(6:7).eq.'WE') then
228  eta=prop(index+4)+bc(nactdog(3,node))*dtheta
229  prop(index+4)=eta
230  nelem=nint(prop(index+7))
231  elseif(lakon(ieg(i))(6:7).eq.'DS') then
232  eta=prop(index+7)+bc(nactdog(3,node))*dtheta
233  prop(index+7)=eta
234  nelem=nint(prop(index+9))
235  endif
236  v(3,node)=eta
237  vama=eta
238 !
239 ! check whether 0<=eta<=1. If not, the hydraulic jump
240 ! does not take place in the element itself and has to
241 ! be forced out of the element by adjusting the
242 ! water depth of one of the end nodes
243 !
244  if((eta.lt.0.d0).or.(eta.gt.1.d0)) then
245  index=ipkon(nelem)
246  node1=kon(index+1)
247  nodem=kon(index+2)
248  node2=kon(index+3)
249  xflow=v(1,nodem)
250 !
251 ! determining the temperature for the
252 ! material properties
253 !
254  if(xflow.gt.0) then
255  if(node1.eq.0) then
256  gastemp=v(0,node2)
257  else
258  gastemp=v(0,node1)
259  endif
260  else
261  if(node2.eq.0) then
262  gastemp=v(0,node1)
263  else
264  gastemp=v(0,node2)
265  endif
266  endif
267 !
268  imat=ielmat(1,nelem)
269 !
270  call materialdata_tg(imat,ntmat_,gastemp,shcon,nshcon,
271  & cp,r,dvi,rhcon,nrhcon,rho)
272 !
273  do j=1,3
274  g(j)=0.d0
275  enddo
276  if(nbody.gt.0) then
277  index=nelem
278  do
279  j=ipobody(1,index)
280  if(j.eq.0) exit
281  if(ibody(1,j).eq.2) then
282  g(1)=g(1)+xbodyact(1,j)*xbodyact(2,j)
283  g(2)=g(2)+xbodyact(1,j)*xbodyact(3,j)
284  g(3)=g(3)+xbodyact(1,j)*xbodyact(4,j)
285  endif
286  index=ipobody(2,index)
287  if(index.eq.0) exit
288  enddo
289  endif
290 !
291  kflag=3
292  call flux(node1,node2,nodem,nelem,lakon,kon,ipkon,
293  & nactdog,identity,
294  & ielprop,prop,kflag,v,xflow,f,nodef,idirf,df,
295  & cp,r,rho,physcon,g,co,dvi,numf,vold,set,shcon,
296  & nshcon,rhcon,nrhcon,ntmat_,mi,ider,ttime,time,
297  & iaxial)
298  kflag=2
299 !
300  endif
301  endif
302  endif
303  enddo
304 !
305 c write(*,*) 'resultnet'
306 c do i=1,ntg
307 c node=itg(i)
308 c write(*,'(i10,3(1x,e11.4))') node,(v(j,node),j=0,2)
309 c enddo
310 !
311 ! testing the validity of the pressures
312 !
313  do i=1,ntg
314  node=itg(i)
315  if(v(2,node).lt.0) then
316  write(*,*) 'wrong pressure; node ',node
317  iin=0
318  return
319  endif
320  enddo
321 !
322 ! testing validity of temperatures
323 !
324  do i=1,ntg
325  node=itg(i)
326  if(v(0,node).lt.0) then
327  iin=0
328  return
329  endif
330  enddo
331 !
332 ! testing the validity of the solution for branches elements
333 ! and restrictor. Since the element properties are dependent on
334 ! a predefined flow direction a change of this will lead to
335 ! wrong head losses
336 !
337  do i=1,nflow
338  nelem=ieg(i)
339  if ((lakon(nelem)(4:5).eq.'ATR').or.
340  & (lakon(nelem)(4:5).eq.'RTA')) then
341  xflow=v(1,kon(ipkon(nelem)+2))
342  if(xflow.lt.0d0)then
343  Write(*,*)'*WARNING in resultnet.f'
344  write(*,*)'Element',nelem,'of TYPE ABSOLUTE TO RELATIVE'
345  write(*,*)'The flow direction is no more conform '
346  write(*,*)'to element definition'
347  write(*,*)'Check the pertinence of the results'
348  endif
349  elseif(lakon(nelem)(2:3).eq.'RE') then
350 !
351  if(lakon(nelem)(4:5).ne.'BR') then
352  nodem=kon(ipkon(nelem)+2)
353  xflow=v(1,nodem)
354  if (xflow.lt.0) then
355  Write(*,*)'*WARNING in resultnet.f'
356  write(*,*)'Element',nelem,'of TYPE RESTRICTOR'
357  write(*,*)'The flow direction is no more conform '
358  write(*,*)'to element definition'
359  write(*,*)'Check the pertinence of the results'
360  endif
361 !
362  elseif(lakon(nelem)(4:5).eq.'BR') then
363  index=ielprop(nelem)
364 !
365  nelem0=nint(prop(index+1))
366  nodem0=kon(ipkon(nelem0)+2)
367  xflow0=v(1,nodem0)
368 !
369  nelem1=nint(prop(index+2))
370  nodem1=kon(ipkon(nelem1)+2)
371  xflow1=v(1,nodem1)
372 !
373  nelem2=nint(prop(index+3))
374  nodem2=kon(ipkon(nelem2)+2)
375  xflow2=v(1,nodem2)
376 !
377  if((xflow0.lt.0).or.(xflow1.lt.0).or.(xflow2.lt.0)) then
378  Write(*,*)'*WARNING in resultnet.f'
379  write(*,*)'Element',nelem,'of TYPE BRANCH'
380  write(*,*)'The flow direction is no more conform '
381  write(*,*)'to element definition'
382  write(*,*)'Check the pertinence of the results'
383  endif
384  endif
385  endif
386  enddo
387 !
388 ! determining the static temperature and checking for
389 ! critical conditions (only for non-chambers, i.e. for
390 ! nodes purely belonging to gas pipes or restrictors except
391 ! restrictor wall orifice)
392 !
393  do i=1,ntg
394  node=itg(i)
395  nelem=ineighe(i)
396  if(nelem.eq.-1) then
397  v(3,node)=v(0,node)
398  elseif(nelem.gt.0) then
399 !
400  nodem=kon(ipkon(nelem)+2)
401  t=v(3,node)
402  tt=v(0,node)
403  pt=v(2,node)
404  xflow=v(1,nodem)
405 !
406  icase=0
407  inv=1
408  imat=ielmat(1,nelem)
409  call materialdata_tg(imat,ntmat_,v(3,node),
410  & shcon,nshcon,cp,r,dvi,rhcon,nrhcon,rho)
411 !
412  index=ielprop(nelem)
413  indexe=ipkon(nelem)
414 !
415  kappa=(cp/(cp-r))
416 !
417  if(lakon(nelem)(2:5).eq.'GAPF') then
418  a=prop(index+1)
419  if(lakon(nelem)(2:6).eq.'GAPFA') then
420  icase=0
421  elseif(lakon(nelem)(2:6).eq.'GAPFI') then
422  icase=1
423  endif
424 !
425  elseif(lakon(nelem)(2:3).eq.'RE') then
426  node1=kon(indexe+1)
427  node2=kon(indexe+3)
428 !
429  if(lakon(nelem)(4:5).eq.'EX') then
430  if(lakon(nint(prop(index+4)))(2:6).eq.'GAPFA') then
431  icase=0
432  elseif(lakon(nint(prop(index+4)))(2:6).eq.'GAPFI')then
433  icase=1
434  endif
435  else
436  icase=0
437  endif
438 !
439 ! defining the sections
440 !
441  if(lakon(nelem)(4:5).eq.'BE') then
442  a=prop(index+1)
443 !
444  elseif(lakon(nelem)(4:5).eq.'BR') then
445  if(lakon(nelem)(4:6).eq.'BRJ') then
446  if(nelem.eq.nint(prop(index+2)))then
447  a=prop(index+5)
448  elseif(nelem.eq.nint(prop(index+3))) then
449  a=prop(index+6)
450  endif
451  elseif(lakon(nelem)(4:6).eq.'BRS') then
452  if(nelem.eq.nint(prop(index+2)))then
453  a=prop(index+5)
454  elseif(nelem.eq.nint(prop(index+3))) then
455  a=prop(index+6)
456  endif
457  endif
458 !
459  else
460 !
461  if(node.eq.node1) then
462  a=prop(index+1)
463  elseif(node.eq.node2) then
464  a=prop(index+2)
465  endif
466  endif
467  elseif(lakon(nelem)(2:3).eq.'UP') then
468 !
469 ! user elements whose names start with UP are assumed
470 ! to be Pipe-like (static and total temperatures differ)
471 ! The cross area is assumed to be the first property
472 !
473  a=prop(index+1)
474  icase=0
475  endif
476 !
477  xflow360=dabs(xflow*iaxial)
478  call ts_calc(xflow360,tt,pt,kappa,r,a,t,icase)
479 !
480  v(3,node)=t
481 !
482  if(lakon(nelem)(2:3).eq.'RE') then
483 !
484 ! check whether the mass flow does not exceed
485 ! critical conditions (only for restrictor elements)
486 !
487  if(xflow.lt.0d0) then
488  inv=-1
489  else
490  inv=1
491  endif
492 !
493  if(icase.eq.0) then
494  qred_crit=dsqrt(kappa/r)*(1.+0.5*(kappa-1.))
495  & **(-0.5d0*(kappa+1.)/(kappa-1.))
496  else
497  qred_crit=dsqrt(1/r)*(1.+0.5*(kappa-1.)/kappa)
498  & **(-0.5d0*(kappa+1.)/(kappa-1.))
499  endif
500  xflow_crit=qred_crit*pt*a/dsqrt(tt)
501 !
502  if(xflow360.ge.xflow_crit) then
503  v(1,nodem)=inv*xflow_crit/iaxial
504  if(icase.eq.1) then
505 !
506  if(nactdog(0,node2).ne.0) then
507  node1=kon(indexe+1)
508  node2=kon(indexe+3)
509  v(3,node2)=v(3,node1)
510  v(0,node2)=v(3,node2)
511  & *(1+0.5d0*(kappa-1)/kappa)
512 
513  endif
514  endif
515  endif
516  endif
517 !
518  endif
519 !
520  enddo
521 !
522 ! testing if the flow direction is conform to the pressure
523 ! gradient
524 !
525  iplausi=1
526  do i=1,nflow
527  nelem=ieg(i)
528  indexe=ipkon(nelem)
529  node1=kon(indexe+1)
530  nodem=kon(indexe+2)
531  node2=kon(indexe+3)
532 !
533  if((node1.ne.0).and.(node2.ne.0)) then
534 !
535 ! the mass flow rates must always be oriented in the direction
536 ! of the decreasing pressure gradient except for ACCTUBE,
537 ! CROSS SPLIT, MOEHRING, ROTATING CAVITY, RADIAL OUTFLOW and
538 ! INLET, RIMSEAL, S-PUMP and VORTEX
539 !
540  if((lakon(nelem)(2:8).ne.'ACCTUBE').and.
541  & (lakon(nelem)(2:5).ne.'CROS').and.
542  & (lakon(nelem)(2:4).ne.'MRG').and.
543  & (lakon(nelem)(2:4).ne.'RCV').and.
544  & (lakon(nelem)(2:3).ne.'RO').and.
545  & (lakon(nelem)(2:5).ne.'RIMS').and.
546  & (lakon(nelem)(2:6).ne.'SPUMP').and.
547  & (lakon(nelem)(2:3).ne.'VO')) then
548  if(nactdog(1,nodem).ne.0) then
549  if(((v(1,nodem).gt.1.d-12).and.
550  & (v(2,node1).lt.0.999d0*v(2,node2))).or.
551  & ((v(1,nodem).lt.-1.d-12).and.
552  & (0.999d0*v(2,node1).gt.v(2,node2)))) then
553  iplausi=0
554  write(*,*) '*WARNING in resultnet: the flow'
555  write(*,*) ' direction is not conform'
556  write(*,*) ' to the pressure gradient'
557  write(*,*) ' element',nelem
558  endif
559  endif
560  endif
561  endif
562  enddo
563 !
564 ! reinitialisation of the Bc matrix
565 !
566  do i=1,nteq
567  bc(i)=0.d0
568  enddo
569 !
570  qat=0.d0
571  qaf=0.d0
572  nalt=0
573  nalf=0
574 !
575 ! determining the residual
576 !
577  do i=1,nflow
578  nelem=ieg(i)
579  index=ipkon(nelem)
580  node1=kon(index+1)
581  nodem=kon(index+2)
582  node2=kon(index+3)
583  xflow=v(1,nodem)
584 !
585 ! gas: the property temperature is the static temperature
586 !
587  if((lakon(nelem)(2:3).ne.'LP').and.
588  & (lakon(nelem)(2:3).ne.'LI')) then
589  if(node1.eq.0) then
590  tg1=v(0,node2)
591  tg2=tg1
592  ts1=v(3,node2)
593  ts2=ts1
594  elseif(node2.eq.0) then
595  tg1=v(0,node1)
596  tg2=tg1
597  ts1=v(3,node1)
598  ts2=ts1
599  else
600  tg1=v(0,node1)
601  tg2=v(0,node2)
602  ts1=v(3,node1)
603  ts2=v(3,node2)
604  endif
605  gastemp=(ts1+ts2)/2.d0
606  else
607 !
608 ! liquid: only one temperature
609 !
610  if(xflow.gt.0) then
611  if(node1.eq.0) then
612  gastemp=v(0,node2)
613  else
614  gastemp=v(0,node1)
615  endif
616  else
617  if(node2.eq.0) then
618  gastemp=v(0,node1)
619  else
620  gastemp=v(0,node2)
621  endif
622  endif
623 !
624  if(node1.eq.0) then
625  tg2=v(0,node2)
626  tg1=tg2
627  elseif(node2.eq.0) then
628  tg1=v(0,node1)
629  tg2=tg1
630  else
631  tg1=v(0,node1)
632  tg2=v(0,node2)
633  endif
634  endif
635 !
636  imat=ielmat(1,nelem)
637 !
638  call materialdata_tg(imat,ntmat_,gastemp,shcon,nshcon,cp,r,dvi,
639  & rhcon,nrhcon,rho)
640  kappa=cp/(cp-r)
641 !
642 ! Definitions of the constant for isothermal flow elements
643 !
644  if(lakon(nelem)(2:6).eq.'GAPFI') then
645  if((node1.ne.0).and.(node2.ne.0)) then
646  a=prop(ielprop(nelem)+1)
647  pt1=v(2,node1)
648  pt2=v(2,node2)
649 !
650  xflow360=xflow*iaxial
651  icase=1
652 !
653  if(pt1.ge.pt2)then
654  if(dabs(tg2/ts2-(1+0.5*(kappa-1)/kappa)).lt.1e-5) then
655  pt2=dabs(xflow360)*dsqrt(tg2*r)/a
656  & *(1+0.5*(kappa-1)/kappa)
657  & **(0.5*(kappa+1)/(kappa-1))
658 !
659  endif
660  tt1=v(0,node1)-physcon(1)
661  tt2=v(0,node2)-physcon(1)
662  t1=v(3,node1)
663  call ts_calc(xflow360,tt1,pt1,kappa,r,a,t1,icase)
664  call ts_calc(xflow360,tt2,pt2,kappa,r,a,t2,icase)
665  v(3,node1)=t1
666  v(3,node2)=t2
667  else
668  pt1=v(2,node2)
669  pt2=v(2,node1)
670  if(v(2,nodem).ge.(pt2/pt1))then
671  pt2=v(2,nodem)*pt1
672  endif
673 !
674  tt1=v(0,node2)-physcon(1)
675  call ts_calc(xflow360,tt1,pt1,kappa,r,a,t1,icase)
676  tt2=v(0,node1)-physcon(1)
677  call ts_calc(xflow360,tt2,pt2,kappa,r,a,t2,icase)
678  endif
679 !
680  endif
681  endif
682 !
683  if(node1.ne.0) then
684 !
685 ! energy equation contribution node1
686 !
687  if (nacteq(0,node1).ne.0) then
688  ieq=nacteq(0,node1)
689 !
690  if(nacteq(3,node1).eq.0) then
691  if (xflow.lt.0d0)then
692  term=cp*(tg1-tg2)*xflow
693  bc(ieq)=bc(ieq)+term
694  qat=qat+dabs(term)
695  nalt=nalt+1
696  endif
697 !
698  elseif(lakon(nelem)(2:6).eq.'GAPFI') then
699  if((nacteq(3,node1).eq.node2)) then
700 !
701  bc(ieq)=(t2-t1)
702 !
703  endif
704  endif
705  endif
706 !
707 ! mass equation contribution node1
708 !
709  if (nacteq(1,node1).ne.0) then
710  ieq=nacteq(1,node1)
711  bc(ieq)=bc(ieq)-xflow
712  qaf=qaf+dabs(xflow)
713  nalf=nalf+1
714  endif
715  endif
716 !
717  if(node2.ne.0) then
718 !
719 ! energy equation contribution node2
720 !
721  if (nacteq(0,node2).ne.0) then
722  ieq=nacteq(0,node2)
723 !
724  if(nacteq(3,node2).eq.0) then
725  if (xflow.gt.0d0)then
726  term=cp*(tg2-tg1)*xflow
727  bc(ieq)=bc(ieq)-term
728  qat=qat+dabs(term)
729  nalt=nalt+1
730  endif
731 !
732  elseif(lakon(nelem)(2:6).eq.'GAPFI') then
733  if(nacteq(3,node2).eq.node1) then
734 !
735  bc(ieq)=(t2-t1)
736 !
737  endif
738  endif
739  endif
740 !
741 ! mass equation contribution node2
742 ! only in the case of an ACCTUBE but not in the case of an ACCTUBO
743 !
744  if(lakon(nelem)(2:8).ne.'ACCTUBE') then
745  if (nacteq(1,node2).ne.0) then
746  ieq=nacteq(1,node2)
747  bc(ieq)=bc(ieq)+xflow
748  qaf=qaf+dabs(xflow)
749  nalf=nalf+1
750  endif
751  else
752  if (nacteq(1,node2).ne.0) then
753  if(nelem.ne.nint(prop(ielprop(nelem)+14))) then
754  ieq=nacteq(1,node2)
755  bc(ieq)=bc(ieq)+xflow-v(0,nodem)
756  qaf=qaf+dabs(xflow)
757  nalf=nalf+1
758  else
759  ieq=nacteq(1,node2)
760  bc(ieq)=bc(ieq)+xflow
761  qaf=qaf+dabs(xflow)
762  nalf=nalf+1
763  endif
764  endif
765  endif
766  endif
767 !
768 ! element equation
769 !
770  if (nacteq(2,nodem).ne.0) then
771  ieq=nacteq(2,nodem)
772 !
773 ! for liquids: determine the gravity vector
774 !
775  if(lakon(nelem)(2:3).eq.'LI') then
776  do j=1,3
777  g(j)=0.d0
778  enddo
779  if(nbody.gt.0) then
780  index=nelem
781  do
782  j=ipobody(1,index)
783  if(j.eq.0) exit
784  if(ibody(1,j).eq.2) then
785  g(1)=g(1)+xbodyact(1,j)*xbodyact(2,j)
786  g(2)=g(2)+xbodyact(1,j)*xbodyact(3,j)
787  g(3)=g(3)+xbodyact(1,j)*xbodyact(4,j)
788  endif
789  index=ipobody(2,index)
790  if(index.eq.0) exit
791  enddo
792  endif
793  endif
794 !
795  call flux(node1,node2,nodem,nelem,lakon,kon,ipkon,
796  & nactdog,identity,
797  & ielprop,prop,kflag,v,xflow,f,nodef,idirf,df,
798  & cp,r,rho,physcon,g,co,dvi,numf,vold,set,shcon,
799  & nshcon,rhcon,nrhcon,ntmat_,mi,ider,ttime,time,iaxial)
800  bc(ieq)=-f
801  endif
802  enddo
803 !
804 ! convection with the walls: contribution to the energy equations
805 !
806  do i=1,nload
807  if(sideload(i)(3:4).eq.'FC') then
808  nelem=nelemload(1,i)
809  index=ipkon(nelem)
810  if(index.lt.0) cycle
811  lakonl=lakon(nelem)
812  node=nelemload(2,i)
813  ieq=nacteq(0,node)
814  if(ieq.eq.0) then
815  cycle
816  endif
817 !
818  call nident(itg,node,ntg,id)
819 !
820 ! calculate the area
821 !
822  read(sideload(i)(2:2),'(i1)') ig
823 !
824 ! number of nodes and integration points in the face
825 !
826  if(lakonl(4:4).eq.'2') then
827  nope=20
828  nopes=8
829  elseif(lakonl(4:4).eq.'8') then
830  nope=8
831  nopes=4
832  elseif(lakonl(4:5).eq.'10') then
833  nope=10
834  nopes=6
835  elseif(lakonl(4:4).eq.'4') then
836  nope=4
837  nopes=3
838  elseif(lakonl(4:5).eq.'15') then
839  nope=15
840  else
841  nope=6
842  endif
843 !
844  if(lakonl(4:5).eq.'8R') then
845  mint2d=1
846  elseif((lakonl(4:4).eq.'8').or.(lakonl(4:6).eq.'20R'))
847  & then
848  if(lakonl(7:7).eq.'A') then
849  mint2d=2
850  else
851  mint2d=4
852  endif
853  elseif(lakonl(4:4).eq.'2') then
854  mint2d=9
855  elseif(lakonl(4:5).eq.'10') then
856  mint2d=3
857  elseif(lakonl(4:4).eq.'4') then
858  mint2d=1
859  endif
860 !
861  if(lakonl(4:4).eq.'6') then
862  mint2d=1
863  if(ig.le.2) then
864  nopes=3
865  else
866  nopes=4
867  endif
868  endif
869  if(lakonl(4:5).eq.'15') then
870  if(ig.le.2) then
871  mint2d=3
872  nopes=6
873  else
874  mint2d=4
875  nopes=8
876  endif
877  endif
878 !
879 ! connectivity of the element
880 !
881  do k=1,nope
882  konl(k)=kon(index+k)
883  enddo
884 !
885 ! coordinates of the nodes belonging to the face
886 !
887  if((nope.eq.20).or.(nope.eq.8)) then
888  do k=1,nopes
889  tl2(k)=v(0,konl(ifaceq(k,ig)))
890  do j=1,3
891  xl2(j,k)=co(j,konl(ifaceq(k,ig)))+
892  & v(j,konl(ifaceq(k,ig)))
893  enddo
894  enddo
895  elseif((nope.eq.10).or.(nope.eq.4)) then
896  do k=1,nopes
897  tl2(k)=v(0,konl(ifacet(k,ig)))
898  do j=1,3
899  xl2(j,k)=co(j,konl(ifacet(k,ig)))+
900  & v(j,konl(ifacet(k,ig)))
901  enddo
902  enddo
903  else
904  do k=1,nopes
905  tl2(k)=v(0,konl(ifacew(k,ig)))
906  do j=1,3
907  xl2(j,k)=co(j,konl(ifacew(k,ig)))+
908  & v(j,konl(ifacew(k,ig)))
909  enddo
910  enddo
911  endif
912 !
913 ! integration to obtain the area and the mean
914 ! temperature
915 !
916  do l=1,mint2d
917  if((lakonl(4:5).eq.'8R').or.
918  & ((lakonl(4:4).eq.'6').and.(nopes.eq.4))) then
919  xi=gauss2d1(1,l)
920  et=gauss2d1(2,l)
921  weight=weight2d1(l)
922  elseif((lakonl(4:4).eq.'8').or.
923  & (lakonl(4:6).eq.'20R').or.
924  & ((lakonl(4:5).eq.'15').and.(nopes.eq.8))) then
925  xi=gauss2d2(1,l)
926  et=gauss2d2(2,l)
927  weight=weight2d2(l)
928  elseif(lakonl(4:4).eq.'2') then
929  xi=gauss2d3(1,l)
930  et=gauss2d3(2,l)
931  weight=weight2d3(l)
932  elseif((lakonl(4:5).eq.'10').or.
933  & ((lakonl(4:5).eq.'15').and.(nopes.eq.6))) then
934  xi=gauss2d5(1,l)
935  et=gauss2d5(2,l)
936  weight=weight2d5(l)
937  elseif((lakonl(4:4).eq.'4').or.
938  & ((lakonl(4:4).eq.'6').and.(nopes.eq.3))) then
939  xi=gauss2d4(1,l)
940  et=gauss2d4(2,l)
941  weight=weight2d4(l)
942  endif
943 !
944  if(nopes.eq.8) then
945  call shape8q(xi,et,xl2,xsj2,xs2,shp2,iflag)
946  elseif(nopes.eq.4) then
947  call shape4q(xi,et,xl2,xsj2,xs2,shp2,iflag)
948  elseif(nopes.eq.6) then
949  call shape6tri(xi,et,xl2,xsj2,xs2,shp2,iflag)
950  else
951  call shape3tri(xi,et,xl2,xsj2,xs2,shp2,iflag)
952  endif
953 !
954  dxsj2=dsqrt(xsj2(1)*xsj2(1)+xsj2(2)*xsj2(2)+
955  & xsj2(3)*xsj2(3))
956  areaj=dxsj2*weight
957 !
958  temp=0.d0
959  do k=1,3
960  coords(k)=0.d0
961  enddo
962  do j=1,nopes
963  temp=temp+tl2(j)*shp2(4,j)
964  do k=1,3
965  coords(k)=coords(k)+xl2(k,j)*shp2(4,j)
966  enddo
967  enddo
968 !
969  sinktemp=v(0,node)
970  if(sideload(i)(5:6).ne.'NU') then
971  h(1)=xloadact(1,i)
972  else
973  read(sideload(i)(2:2),'(i1)') jltyp
974  jltyp=jltyp+10
975  call film(h,sinktemp,temp,istep,
976  & iinc,tvar,nelem,l,coords,jltyp,field,nfield,
977  & sideload(i),node,areaj,v,mi,ipkon,kon,lakon,
978  & iponoel,inoel,ielprop,prop,ielmat,shcon,nshcon,
979  & rhcon,nrhcon,ntmat_,cocon,ncocon,
980  & ipobody,xbodyact,ibody,heatnod,heatfac)
981 c if(nmethod.eq.1) h(1)=xloadold(1,i)+
982 c & (h(1)-xloadold(1,i))*reltime
983  endif
984  if(lakonl(5:7).eq.'0RA') then
985  term=2.d0*((temp-sinktemp)*h(1)+heatnod)*dxsj2*weight
986  bc(ieq)=bc(ieq)+term
987  qat=qat+dabs(term)
988  nalt=nalt+1
989  else
990  term=((temp-sinktemp)*h(1)+heatnod)*dxsj2*weight
991  bc(ieq)=bc(ieq)+term
992  qat=qat+dabs(term)
993  nalt=nalt+1
994  endif
995  enddo
996  endif
997  enddo
998 !
999 ! prescribed heat generation: contribution the energy equations
1000 !
1001  do i=1,ntg
1002  node=itg(i)
1003  idof=8*(node-1)
1004  call nident(ikforc,idof,nforc,id)
1005  if(id.gt.0) then
1006  if(ikforc(id).eq.idof) then
1007  ieq=nacteq(0,node)
1008  if(ieq.ne.0) then
1009  term=xforcact(ilforc(id))
1010  bc(ieq)=bc(ieq)+term
1011  qat=qat+dabs(term)
1012  nalt=nalt+1
1013  endif
1014  cycle
1015  endif
1016  endif
1017  enddo
1018 !
1019 ! in the case of forced vortices, when temperature change
1020 ! is required, additional heat input is added in the energy
1021 ! equation for the downstream node
1022 !
1023  do i=1,nflow
1024  nelem=ieg(i)
1025  if(lakon(nelem)(2:3).ne.'VO') cycle
1026 !
1027 ! free vortex and no temperature change
1028 !
1029  if((lakon(nelem)(2:5).eq.'VOFR').and.
1030  & (nint(prop(ielprop(nelem)+8)).eq.0)) cycle
1031 !
1032 ! free vortex and temperature change in the absolute system
1033 !
1034  if((lakon(nelem)(2:5).eq.'VOFR').and.
1035  & (nint(prop(ielprop(nelem)+8)).eq.1)) cycle
1036 !
1037 ! forced vortex and no temperature change
1038 !
1039  if((lakon(nelem)(2:5).eq.'VOFO').and.
1040  & (nint(prop(ielprop(nelem)+6)).eq.0)) cycle
1041 !
1042  nodem=kon(ipkon(nelem)+2)
1043  xflow=v(1,nodem)
1044  if(xflow.gt.0d0) then
1045  node1=kon(ipkon(nelem)+1)
1046  node2=kon(ipkon(nelem)+3)
1047  else
1048  node1=kon(ipkon(nelem)+1)
1049  node2=kon(ipkon(nelem)+3)
1050  endif
1051 !
1052  if(xflow.gt.0d0) then
1053  r1=prop(ielprop(nelem)+2)
1054  r2=prop(ielprop(nelem)+1)
1055  if(r1.gt.r2) then
1056  rout=r2
1057  rin=r1
1058  else
1059  rout=r2
1060  rin=r1
1061  endif
1062  else
1063  r1=prop(ielprop(nelem)+2)
1064  r2=prop(ielprop(nelem)+1)
1065  if(r1.gt.r2) then
1066  rout=r1
1067  rin=r2
1068  else
1069  rout=r1
1070  rin=r2
1071  endif
1072  endif
1073 !
1074 ! computing temperature corrected Cp=Cp(T) coefficient
1075 !
1076  tg1=v(0,node1)
1077  tg2=v(0,node2)
1078  if((lakon(nelem)(2:3).ne.'LP').and.
1079  & (lakon(nelem)(2:3).ne.'LI')) then
1080  gastemp=(tg1+tg2)/2.d0
1081  else
1082  if(xflow.gt.0) then
1083  gastemp=tg1
1084  else
1085  gastemp=tg2
1086  endif
1087  endif
1088 !
1089  imat=ielmat(1,nelem)
1090  call materialdata_tg(imat,ntmat_,gastemp,
1091  & shcon,nshcon,cp,r,dvi,rhcon,nrhcon,rho)
1092 !
1093  call cp_corrected(cp,tg1,tg2,cp_cor)
1094 !
1095  uout=prop(ielprop(nelem)+5)*rout
1096  uin=prop(ielprop(nelem)+5)*rin
1097 !
1098 ! free and forced vortices with temperature
1099 ! change in the relative system of coordinates
1100 !
1101  if((lakon(nelem)(2:5).eq.'VOFR') .and.
1102  & (nint(prop(ielprop(nelem)+8)).eq.(-1))) then
1103 !
1104  uout=prop(ielprop(nelem)+7)*rout
1105  uin=prop(ielprop(nelem)+7)*rin
1106 !
1107  heat=0.5d0*cp/cp_cor*(uout**2-uin**2)*xflow
1108 !
1109  elseif (((lakon(nelem)(2:5).eq.'VOFO')
1110  & .and.(nint(prop(ielprop(nelem)+6)).eq.(-1)))) then
1111 !
1112  uout=prop(ielprop(nelem)+5)*rout
1113  uin=prop(ielprop(nelem)+5)*rin
1114 !
1115  heat=0.5d0*cp/cp_cor*(uout**2-uin**2)*xflow
1116 !
1117 ! forced vortices with temperature change in the absolute system
1118 !
1119  elseif((lakon(nelem)(2:5).eq.'VOFO')
1120  & .and.((nint(prop(ielprop(nelem)+6)).eq.1))) then
1121 !
1122  heat=cp/cp_cor*(uout**2-uin**2)*xflow
1123 !
1124  endif
1125 !
1126 ! including the resulting additional heat flux in the energy equation
1127 !
1128  if(xflow.gt.0d0)then
1129  ieq=nacteq(0,node2)
1130  if(ieq.ne.0) then
1131  bc(ieq)=bc(ieq)+heat
1132  qat=qat+dabs(heat)
1133  nalt=nalt+1
1134  endif
1135  else
1136  ieq=nacteq(0,node1)
1137  if(ieq.ne.0) then
1138  bc(ieq)=bc(ieq)+heat
1139  qat=qat+dabs(heat)
1140  nalt=nalt+1
1141  endif
1142  endif
1143  enddo
1144 !
1145 ! transfer element ABSOLUTE TO RELATIVE / RELATIVE TO ABSOLUTE
1146 !
1147  do i= 1, nflow
1148  nelem=ieg(i)
1149 !
1150  if((lakon(nelem)(2:4).eq.'ATR').or.
1151  & (lakon(nelem)(2:4).eq.'RTA')) then
1152 !
1153  nodem=kon(ipkon(nelem)+2)
1154  xflow=v(1,nodem)
1155  if(xflow.gt.0d0) then
1156  node1=kon(ipkon(nelem)+1)
1157  node2=kon(ipkon(nelem)+3)
1158  else
1159  node1=kon(ipkon(nelem)+1)
1160  node2=kon(ipkon(nelem)+3)
1161  endif
1162 !
1163 ! computing temperature corrected Cp=Cp(T) coefficient
1164 !
1165  tg1=v(0,node1)
1166  tg2=v(0,node2)
1167  if((lakon(nelem)(2:3).ne.'LP').and.
1168  & (lakon(nelem)(2:3).ne.'LI')) then
1169  gastemp=(tg1+tg2)/2.d0
1170  else
1171  if(xflow.gt.0) then
1172  gastemp=tg1
1173  else
1174  gastemp=tg2
1175  endif
1176  endif
1177 !
1178  imat=ielmat(1,nelem)
1179  call materialdata_tg(imat,ntmat_,gastemp,
1180  & shcon,nshcon,cp,r,dvi,rhcon,nrhcon,rho)
1181 !
1182  call cp_corrected(cp,tg1,tg2,cp_cor)
1183 !
1184  index=ielprop(nelem)
1185  u=prop(index+1)
1186  ct=prop(index+2)
1187 !
1188 ! if a swirl element was given by the user, this takes
1189 ! precendence to a value of ct
1190 !
1191 c if(ct.eq.0) then
1192  if(nint(prop(index+3)).ne.0) then
1193  nelemswirl=nint(prop(index+3))
1194  index2=ielprop(nelemswirl)
1195 !
1196 ! previous element is a preswirl nozzle
1197 !
1198  if(lakon(nelemswirl)(2:5).eq.'ORPN') then
1199  ct=prop(index2+4)
1200 !
1201 ! previous element is a forced vortex
1202 !
1203  elseif(lakon(nelemswirl)(2:5).eq.'VOFO') then
1204  ct=prop(index2+7)
1205 !
1206 ! previous element is a free vortex
1207 !
1208  elseif(lakon(nelemswirl)(2:5).eq.'VOFR') then
1209  ct=prop(index2+9)
1210  endif
1211  endif
1212 !
1213  if(lakon(nelem)(2:4).eq.'ATR') then
1214  heat=cp/cp_cor*(0.5d0*(u**2-2d0*u*ct)*xflow)
1215 !
1216  elseif(lakon(nelem)(2:4).eq.'RTA') then
1217  heat=cp/cp_cor*(-0.5d0*(u**2-2d0*u*ct)*xflow)
1218  endif
1219 !
1220 ! including the resulting additional heat flux in the energy equation
1221 !
1222  if(xflow.gt.0d0)then
1223  ieq=nacteq(0,node2)
1224  if(ieq.ne.0) then
1225  bc(ieq)=bc(ieq)+heat
1226  qat=qat+dabs(heat)
1227  nalt=nalt+1
1228  endif
1229  else
1230  ieq=nacteq(0,node1)
1231  if(ieq.ne.0) then
1232  bc(ieq)=bc(ieq)+heat
1233  qat=qat+dabs(heat)
1234  nalt=nalt+1
1235  endif
1236  endif
1237  endif
1238  enddo
1239 !
1240 ! additional multiple point constraints
1241 !
1242  j=nteq+1
1243  do i=nmpc,1,-1
1244  if(labmpc(i)(1:7).ne.'NETWORK') cycle
1245  j=j-1
1246  if(labmpc(i)(8:8).eq.' ') then
1247 !
1248 ! linear equation
1249 !
1250  index=ipompc(i)
1251 !
1252  do
1253  node=nodempc(1,index)
1254  idir=nodempc(2,index)
1255  bc(j)=bc(j)-v(idir,node)*coefmpc(index)
1256  index=nodempc(3,index)
1257  if(index.eq.0) exit
1258  enddo
1259  else
1260 !
1261 ! user-defined network equation
1262 !
1263  call networkmpc_rhs(i,ipompc,nodempc,coefmpc,labmpc,
1264  & v,bc,j,mi)
1265  endif
1266  enddo
1267 !
1268 ! determining the typical energy flow
1269 !
1270  if(nalt.gt.0) qat=qat/nalt
1271 !
1272 ! determining the typical mass flow
1273 !
1274  if(nalf.gt.0) qaf=qaf/nalf
1275 !
1276 ! max. residual energy flow, residual mass flow
1277 ! and residuals from the element equations
1278 !
1279  ramt=0.d0
1280  ramf=0.d0
1281  ramp=0.d0
1282  do i=1,ntg
1283  node=itg(i)
1284  if(nacteq(0,node).ne.0) then
1285  ramt=max(ramt,dabs(bc(nacteq(0,node))))
1286  endif
1287  if(nacteq(1,node).ne.0) then
1288  ramf=max(ramf,dabs(bc(nacteq(1,node))))
1289  endif
1290  if(nacteq(2,node).ne.0) then
1291  ramp=max(ramp,dabs(bc(nacteq(2,node))))
1292  endif
1293  enddo
1294 !
1295 c write(30,*) 'bc in resultnet'
1296 c do i=1,9
1297 c write(30,'(1x,e11.4)') bc(i)
1298 c enddo
1299 !
1300  return
shape8q
subroutine shape8q(xi, et, xl, xsj, xs, shp, iflag)
Definition: shape8q.f:20
max
#define max(a, b)
Definition: cascade.c:32
shape3tri
subroutine shape3tri(xi, et, xl, xsj, xs, shp, iflag)
Definition: shape3tri.f:20
df
subroutine df(x, u, uprime, rpar, nev)
Definition: subspace.f:133
flux
subroutine flux(node1, node2, nodem, nelem, lakon, kon, ipkon, nactdog, identity, ielprop, prop, kflag, v, xflow, f, nodef, idirf, df, cp, R, rho, physcon, g, co, dvi, numf, vold, set, shcon, nshcon, rhcon, nrhcon, ntmat_, mi, ider, ttime, time, iaxial)
Definition: flux.f:24
ts_calc
subroutine ts_calc(xflow, Tt, Pt, kappa, r, A, Ts, icase)
Definition: ts_calc.f:20
r1
static double * r1
Definition: filtermain.c:48
shape4q
subroutine shape4q(xi, et, xl, xsj, xs, shp, iflag)
Definition: shape4q.f:20
cp_corrected
subroutine cp_corrected(cp, Tg1, Tg2, cp_cor)
Definition: cp_corrected.f:24
nident
subroutine nident(x, px, n, id)
Definition: nident.f:26
networkmpc_rhs
subroutine networkmpc_rhs(i, ipompc, nodempc, coefmpc, labmpc, v, bc, j, mi)
Definition: networkmpc_rhs.f:23
shape6tri
subroutine shape6tri(xi, et, xl, xsj, xs, shp, iflag)
Definition: shape6tri.f:20
film
subroutine film(h, sink, temp, kstep, kinc, time, noel, npt, coords, jltyp, field, nfield, loadtype, node, area, vold, mi, ipkon, kon, lakon, iponoel, inoel, ielprop, prop, ielmat, shcon, nshcon, rhcon, nrhcon, ntmat_, cocon, ncocon, ipobody, xbody, ibody, heatnod, heatfac)
Definition: film.f:24
materialdata_tg
subroutine materialdata_tg(imat, ntmat_, t1l, shcon, nshcon, sph, r, dvi, rhcon, nrhcon, rho)
Definition: materialdata_tg.f:20
Hosted by OpenAircraft.com, (Michigan UAV, LLC)