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

Go to the source code of this file.

Functions/Subroutines

subroutine preinitialnet (ieg, lakon, v, ipkon, kon, nflow, prop, ielprop, ielmat, ntmat_, shcon, nshcon, rhcon, nrhcon, mi, iponoel, inoel, itg, ntg)
 

Function/Subroutine Documentation

◆ preinitialnet()

subroutine preinitialnet ( integer, dimension(*)  ieg,
character*8, dimension(*)  lakon,
real*8, dimension(0:mi(2),*)  v,
integer, dimension(*)  ipkon,
integer, dimension(*)  kon,
integer  nflow,
real*8, dimension(*)  prop,
integer, dimension(*)  ielprop,
integer, dimension(mi(3),*)  ielmat,
integer  ntmat_,
real*8, dimension(0:3,ntmat_,*)  shcon,
integer, dimension(*)  nshcon,
real*8, dimension(0:1,ntmat_,*)  rhcon,
integer, dimension(*)  nrhcon,
integer, dimension(*)  mi,
integer, dimension(*)  iponoel,
integer, dimension(2,*)  inoel,
integer, dimension(*)  itg,
integer  ntg 
)
22 !
23 ! this routine only applies to compressible networks
24 !
25 ! determination of initial values based on the boundary conditions
26 ! and the initial values given by the user by propagating these
27 ! through the network (using information on the mass flow direction
28 ! derived from unidirectional network elements or mass flow given
29 ! by the user (boundary conditions or initial conditions))
30 !
31 ! a mass flow with size 1.d-30 is used to propagate the sign
32 ! of the mass flow (in case only the sign and not the size is known)
33 !
34 ! it is assumed that pressures, temperatures and mass flows cannot
35 ! be identically zero (a zero pressure does not make sense,
36 ! temperature has to be the absolute temperature,
37 ! a zero mass flow leads to convergence problems).
38 !
39  implicit none
40 !
41  character*8 lakon(*)
42 !
43  integer mi(*),ieg(*),nflow,i,ielmat(mi(3),*),ntmat_,node1,node2,
44  & nelem,index,nshcon(*),ipkon(*),kon(*),nodem,imat,ielprop(*),
45  & nrhcon(*),neighbor,ichange,iponoel(*),inoel(2,*),indexe,
46  & itg(*),ntg,node,imin,imax,iel,nodemnei,ierror,nelemnei,
47  & nodenei,ibranch,numel,noderef,nelemref,ierr
48 !
49  real*8 prop(*),shcon(0:3,ntmat_,*),xflow,v(0:mi(2),*),cp,r,
50  & dvi,rho,rhcon(0:1,ntmat_,*),kappa,cti,ti,ri,ro,p1zp2,omega,
51  & p2zp1,xmin,xmax,fluxtot,ratio,pref,prefnew
52 !
53 ! the user should assign an initial pressure to any
54 ! - node which is connected to an inlet or an outlet
55 ! - node belonging to more than 2 network elements
56 ! this is checked in the next lines
57 !
58  ierror=0
59 !
60  do i=1,ntg
61  node=itg(i)
62  index=iponoel(node)
63 !
64 ! node not connected to any element
65 !
66  if(index.eq.0) cycle
67  nelem=inoel(1,index)
68 !
69 ! midside node
70 !
71  if(kon(ipkon(nelem)+2).eq.node) cycle
72 !
73 ! initial pressure assigned
74 !
75  if(v(2,node).gt.0.d0) cycle
76 !
77 ! check whether node belongs to more than 1 element
78 !
79  index=inoel(2,index)
80  if(index.eq.0) then
81  write(*,*) '*ERROR in preinitialnet:'
82  write(*,*) ' node',node,
83  & ' is an inlet or outlet, yet'
84  write(*,*) ' no initial pressure was assigned'
85  ierror=1
86  cycle
87  endif
88 !
89  call networkneighbor(nelem,node,nelemnei,nodenei,ibranch,
90  & iponoel,inoel,ipkon,kon)
91  if(nodenei.eq.0) then
92  write(*,*) '*ERROR in preinitialnet:'
93  write(*,*) ' node',node,
94  & ' is an inlet or outlet, yet'
95  write(*,*) ' no initial pressure was assigned'
96  ierror=1
97  cycle
98  endif
99 !
100  index=inoel(2,index)
101  if(index.ne.0) then
102  write(*,*) '*ERROR in preinitialnet:'
103  write(*,*) ' node',node,
104  & ' belongs to more than 2 network elements, yet'
105  write(*,*) ' no initial pressure was assigned'
106  ierror=1
107  endif
108  enddo
109  if(ierror.eq.1) call exit(201)
110 !
111 ! for directional elements: small mass flow if none specified
112 !
113  do i=1,nflow
114  nelem=ieg(i)
115  indexe=ipkon(nelem)
116 !
117  nodem=kon(indexe+2)
118 !
119  if(lakon(nelem)(2:3).eq.'OR') then
120 !
121 ! orifice
122 !
123  if(v(1,nodem).eq.0.d0) v(1,nodem)=1.d-30
124  elseif(lakon(nelem)(2:4).eq.'LAB') then
125 !
126 ! stepped labyrinth
127 !
128  if((lakon(nelem)(5:6).eq.'SP').or.
129  & (lakon(nelem)(6:7).eq.'SP')) then
130  if(v(1,nodem).eq.0.d0) v(1,nodem)=1.d-30
131  endif
132  elseif(lakon(nelem)(2:3).eq.'RE') then
133 !
134 ! enlargement, contraction, wall orifice, entrance or exit
135 !
136  if((lakon(nelem)(4:5).eq.'EL').or.
137  & (lakon(nelem)(4:5).eq.'CO').or.
138  & (lakon(nelem)(4:7).eq.'WAOR').or.
139  & (lakon(nelem)(4:5).eq.'EN').or.
140  & (lakon(nelem)(4:5).eq.'EX')) then
141  if(v(1,nodem).eq.0.d0) v(1,nodem)=1.d-30
142  endif
143  elseif(lakon(nelem)(4:5).eq.'BR') then
144 !
145 ! joint or split
146 !
147  if((lakon(nelem)(6:6).eq.'J').or.
148  & (lakon(nelem)(6:6).eq.'S')) then
149  if(v(1,nodem).eq.0.d0) v(1,nodem)=1.d-30
150  endif
151  elseif(lakon(nelem)(2:7).eq.'CROSPL') then
152 !
153 ! cross split
154 !
155  if(v(1,nodem).eq.0.d0) v(1,nodem)=1.d-30
156  elseif(lakon(nelem)(2:3).eq.'MR') then
157 !
158 ! Moehring
159 !
160  if(v(1,nodem).eq.0.d0) v(1,nodem)=1.d-30
161  endif
162  enddo
163 !
164  do
165  ichange=0
166 !
167 ! propagation of the pressure through the network
168 !
169  do i=1,nflow
170  nelem=ieg(i)
171  indexe=ipkon(nelem)
172 !
173  node1=kon(indexe+1)
174  if(node1.eq.0) cycle
175  nodem=kon(indexe+2)
176  node2=kon(indexe+3)
177  if(node2.eq.0) cycle
178 !
179 ! at least one total pressure value unknown in the element
180 !
181  if(((v(2,node1).ne.0.d0).and.(v(2,node2).eq.0.d0)).or.
182  & ((v(2,node1).eq.0.d0).and.(v(2,node2).ne.0.d0))) then
183 !
184  if(lakon(nelem)(2:3).eq.'VO') then
185 !
186 ! vortex: pressure ratio can be determined
187 ! from geometry
188 !
189  index=ielprop(nelem)
190 !
191  if(prop(index+1).lt.prop(index+2)) then
192 !
193 ! r2 < r1
194 !
195  if(v(0,node2).ne.0.d0) then
196  ri=prop(index+1)
197  ro=prop(index+2)
198  ti=v(0,node2)
199  if(lakon(nelem)(4:5).eq.'FO') then
200  omega=prop(index+5)
201  else
202  omega=prop(index+7)
203  endif
204 !
205  imat=ielmat(1,nelem)
206  call materialdata_tg(imat,ntmat_,ti,shcon,
207  & nshcon,cp,
208  & r,dvi,rhcon,nrhcon,rho)
209 !
210  kappa=cp/(cp-r)
211  cti=omega*ri
212 !
213  if(lakon(nelem)(4:5).eq.'FO') then
214 !
215 ! forced vortex
216 !
217  p1zp2=(1.d0+cti**2*((ro/ri)**2-1.d0)/
218  & (2.d0*cp*ti))**(kappa/(kappa-1.d0))
219  else
220 !
221 ! free vortex
222 !
223  p1zp2=(1.d0+cti**2*(1.d0-(ri/ro)**2)/
224  & (2.d0*cp*ti))**(kappa/(kappa-1.d0))
225  endif
226 !
227  if(v(2,node1).eq.0.d0) then
228  v(2,node1)=v(2,node2)*p1zp2
229  else
230  v(2,node2)=v(2,node1)/p1zp2
231  endif
232  ichange=1
233  endif
234  else
235 !
236 ! r1 <= r2
237 !
238  if(v(0,node1).ne.0.d0) then
239  ri=prop(index+2)
240  ro=prop(index+1)
241  ti=v(0,node1)
242  if(lakon(nelem)(4:5).eq.'FO') then
243  omega=prop(index+5)
244  else
245  omega=prop(index+7)
246  endif
247 !
248  imat=ielmat(1,nelem)
249  call materialdata_tg(imat,ntmat_,ti,shcon,
250  & nshcon,cp,
251  & r,dvi,rhcon,nrhcon,rho)
252 !
253  kappa=cp/(cp-r)
254  cti=omega*ri
255 !
256  if(lakon(nelem)(4:5).eq.'FO') then
257 !
258 ! forced vortex
259 !
260  p2zp1=(1.d0+cti**2*((ro/ri)**2-1.d0)/
261  & (2.d0*cp*ti))**(kappa/(kappa-1.d0))
262  else
263 !
264 ! free vortex
265 !
266  p2zp1=(1.d0+cti**2*(1.d0-(ri/ro)**2)/
267  & (2.d0*cp*ti))**(kappa/(kappa-1.d0))
268  endif
269 !
270  if(v(2,node1).eq.0.d0) then
271  v(2,node1)=v(2,node2)/p2zp1
272  else
273  v(2,node2)=v(2,node1)*p2zp1
274  endif
275  ichange=1
276  endif
277  endif
278  elseif(v(1,nodem).ne.0.d0) then
279 !
280 ! mass flow is given (either size or just the
281 ! direction): small slope
282 !
283  ierror=0
284  if(v(1,nodem).gt.0.d0) then
285  if(v(2,node1).eq.0.d0) then
286  v(2,node1)=v(2,node2)*1.01d0
287  call networkneighbor(nelem,node1,nelemnei,
288  & nodenei,ibranch,iponoel,inoel,ipkon,kon)
289  if(v(2,nodenei).le.v(2,node1)) ierror=1
290  else
291  v(2,node2)=v(2,node1)*0.99d0
292  call networkneighbor(nelem,node2,nelemnei,
293  & nodenei,ibranch,iponoel,inoel,ipkon,kon)
294  if(v(2,nodenei).ge.v(2,node2)) ierror=2
295  endif
296  else
297  if(v(2,node1).eq.0.d0) then
298  v(2,node1)=v(2,node2)*0.99d0
299  call networkneighbor(nelem,node1,nelemnei,
300  & nodenei,ibranch,iponoel,inoel,ipkon,kon)
301  if(v(2,nodenei).ge.v(2,node1)) ierror=1
302  else
303  v(2,node2)=v(2,node1)*1.01d0
304  call networkneighbor(nelem,node2,nelemnei,
305  & nodenei,ibranch,iponoel,inoel,ipkon,kon)
306  if(v(2,nodenei).le.v(2,node2)) ierror=2
307  endif
308  endif
309 !
310 ! if a discontinuity is detected in the pressure
311 ! the complete branch (i.e. the elements between
312 ! branch points and/or inlets and/or outlets) is
313 ! reanalyzed
314 !
315  if(ierror.ne.0) then
316 !
317 ! looking in the direction of node 1 until a
318 ! branch point/inlet/outlet
319 !
320  node=node1
321  do
322  call networkneighbor(nelem,node,nelemnei,
323  & nodenei,ibranch,iponoel,inoel,ipkon,kon)
324  if((ibranch.eq.1).or.(nodenei.eq.0)) exit
325  node=nodenei
326  nelem=nelemnei
327  enddo
328 !
329 ! noderef is branch point/inlet/outlet
330 ! the adjacent element into the branch is nelemref
331 !
332  noderef=node
333  nelemref=nelem
334 !
335  indexe=ipkon(nelemref)
336  if(kon(indexe+1).eq.noderef) then
337  node=kon(indexe+3)
338  else
339  node=kon(indexe+1)
340  endif
341 !
342 ! looking in the other direction until
343 ! branch point/inlet/outlet
344 ! counting the non-vortex elements
345 ! storing the pressure ratio over the vortices
346 !
347  ratio=1.d0
348  numel=0
349 !
350  if(lakon(nelem)(2:3).eq.'VO') then
351  ratio=ratio*v(2,node)/v(2,noderef)
352  else
353  numel=numel+1
354  endif
355 !
356  ierr=0
357  do
358  call networkneighbor(nelem,node,nelemnei,
359  & nodenei,ibranch,iponoel,inoel,ipkon,kon)
360  if((ibranch.eq.1).or.(nodenei.eq.0)) exit
361  if(lakon(nelemnei)(2:3).eq.'VO') then
362  if((v(2,node).eq.0.d0).or.
363  & (v(2,nodenei).eq.0.d0)) then
364  ierr=1
365  exit
366  endif
367  ratio=ratio*v(2,nodenei)/v(2,node)
368  else
369  numel=numel+1
370  endif
371  node=nodenei
372  nelem=nelemnei
373  enddo
374 !
375  if(ierr.eq.0) then
376 !
377 ! determining the required pressure ratio over the
378 ! non-vortex elements from the pressure at the ends of
379 ! the branch and the pressure ratio over the vortices
380 !
381  ratio=(v(2,node)/(v(2,noderef)*ratio))
382  & **(1.d0/numel)
383 !
384 ! going through the branch again; determining the
385 ! initial values
386 !
387  indexe=ipkon(nelemref)
388  if(kon(indexe+1).eq.noderef) then
389  node=kon(indexe+3)
390  else
391  node=kon(indexe+1)
392  endif
393 !
394  nelem=nelemref
395  pref=v(2,node)
396  if(lakon(nelem)(2:3).ne.'VO') then
397  v(2,node)=v(2,noderef)*ratio
398  endif
399 !
400  do
401  call networkneighbor(nelem,node,nelemnei,
402  & nodenei,ibranch,iponoel,inoel,ipkon,kon)
403  if((ibranch.eq.1).or.(nodenei.eq.0)) exit
404  if(lakon(nelemnei)(2:3).eq.'VO') then
405  prefnew=v(2,nodenei)
406  v(2,nodenei)=v(2,node)*v(2,nodenei)/pref
407  pref=prefnew
408  else
409  pref=v(2,nodenei)
410  v(2,nodenei)=v(2,node)*ratio
411 c write(*,*) nodenei,v(2,nodenei)
412  endif
413  node=nodenei
414  nelem=nelemnei
415  enddo
416  else
417 !
418 ! the pressure in the nodes adjacent to at least one
419 ! vortex element were not available
420 !
421 ! reset values; no change;
422 !
423  if(ierror.eq.1) then
424  v(2,node1)=0.d0
425  else
426  v(2,node2)=0.d0
427  endif
428  cycle
429  endif
430  endif
431 !
432  ichange=1
433  endif
434  endif
435  enddo
436 !
437 ! propagation of the mass flow through the network
438 !
439  loop1: do i=1,nflow
440  nelem=ieg(i)
441  indexe=ipkon(nelem)
442  nodem=kon(indexe+2)
443 !
444  if(dabs(v(1,nodem)).le.1.d-30) then
445 !
446 ! no initial mass flow given yet
447 ! check neighbors for mass flow (only if not
448 ! branch nor joint)
449 !
450 ! first end node
451 !
452  node1=kon(indexe+1)
453 !
454  if(node1.ne.0) then
455  index=iponoel(node1)
456 !
457  if(inoel(2,inoel(2,index)).eq.0) then
458 !
459 ! no branch nor joint; determine neighboring element
460 !
461  if(inoel(1,index).eq.nelem) then
462  neighbor=inoel(1,inoel(2,index))
463  else
464  neighbor=inoel(1,index)
465  endif
466 !
467 ! initial mass flow in neighboring element
468 !
469  xflow=v(1,kon(ipkon(neighbor)+2))
470 !
471  if(dabs(v(1,nodem)).gt.0.d0) then
472 !
473 ! propagate initial mass flow
474 !
475  if(dabs(xflow).gt.1.d-30) then
476  if(kon(ipkon(neighbor)+1).eq.node1) then
477  v(1,nodem)=-xflow
478  else
479  v(1,nodem)=xflow
480  endif
481  ichange=1
482  cycle
483  endif
484  else
485 !
486 ! propagate only the sign of the mass flow
487 !
488  if(dabs(xflow).gt.0.d0) then
489  if(kon(ipkon(neighbor)+1).eq.node1) then
490  v(1,nodem)=-xflow
491  else
492  v(1,nodem)=xflow
493  endif
494  ichange=1
495  cycle
496  endif
497  endif
498 !
499 ! if more than 2 elements meet: check whether
500 ! the flux in all but the element at stake is
501 ! known. If so, apply the mass balance
502 !
503  fluxtot=0.d0
504  do
505  if(inoel(1,index).ne.nelem) then
506  iel=inoel(1,index)
507  nodemnei=kon(ipkon(iel)+2)
508  if(dabs(v(1,nodemnei)).le.1.d-30) exit
509 !
510 ! convention: inflow = positive
511 !
512  if(kon(ipkon(iel)+1).eq.node1) then
513  fluxtot=fluxtot-v(1,nodemnei)
514  else
515  fluxtot=fluxtot+v(1,nodemnei)
516  endif
517  endif
518  if(inoel(2,index).eq.0) then
519  v(1,nodem)=fluxtot
520  ichange=1
521  cycle loop1
522  else
523  index=inoel(2,index)
524  endif
525  enddo
526 !
527  endif
528  endif
529 !
530 ! second end node
531 !
532  node2=kon(indexe+3)
533 !
534  if(node2.ne.0) then
535  index=iponoel(node2)
536 !
537  if(inoel(2,inoel(2,index)).eq.0) then
538 !
539 ! no branch nor joint; determine neighboring element
540 !
541  if(inoel(1,index).eq.nelem) then
542  neighbor=inoel(1,inoel(2,index))
543  else
544  neighbor=inoel(1,index)
545  endif
546 !
547 ! initial mass flow in neighboring element
548 !
549  xflow=v(1,kon(ipkon(neighbor)+2))
550 !
551  if(dabs(v(1,nodem)).gt.0.d0) then
552 !
553 ! propagate initial mass flow
554 !
555  if(dabs(xflow).gt.1.d-30) then
556  if(kon(ipkon(neighbor)+3).eq.node2) then
557  v(1,nodem)=-xflow
558  else
559  v(1,nodem)=xflow
560  endif
561  ichange=1
562  cycle
563  endif
564  else
565 !
566 ! propagate only the sign of the mass flow
567 !
568  if(dabs(xflow).gt.0.d0) then
569  if(kon(ipkon(neighbor)+3).eq.node2) then
570  v(1,nodem)=-xflow
571  else
572  v(1,nodem)=xflow
573  endif
574  ichange=1
575  cycle
576  endif
577  endif
578 !
579 ! if more than 2 elements meet: check whether
580 ! the flux in all but the element at stake is
581 ! known. If so, apply the mass balance
582 !
583  fluxtot=0.d0
584  do
585  if(inoel(1,index).ne.nelem) then
586  iel=inoel(1,index)
587  nodemnei=kon(ipkon(iel)+2)
588  if(dabs(v(1,nodemnei)).le.1.d-30) exit
589 !
590 ! convention: outflow = positive
591 !
592  if(kon(ipkon(iel)+3).eq.node2) then
593  fluxtot=fluxtot-v(1,nodemnei)
594  else
595  fluxtot=fluxtot+v(1,nodemnei)
596  endif
597  endif
598  if(inoel(2,index).eq.0) then
599  v(1,nodem)=fluxtot
600  ichange=1
601  cycle loop1
602  else
603  index=inoel(2,index)
604  endif
605  enddo
606 !
607  endif
608  endif
609  endif
610  enddo loop1
611 !
612 ! propagation of the temperature
613 !
614  do i=1,nflow
615  nelem=ieg(i)
616  indexe=ipkon(nelem)
617  node1=kon(indexe+1)
618  if(node1.eq.0) cycle
619  node2=kon(indexe+3)
620  if(node2.eq.0) cycle
621 !
622  if(((v(0,node1).ne.0.d0).and.(v(0,node2).ne.0.d0)).or.
623  & ((v(0,node1).eq.0.d0).and.(v(0,node2).eq.0.d0))) cycle
624 !
625 ! If the element is a adiabatic gas pipe the
626 ! total temperature at both ends is equal
627 !
628  if(lakon(nelem)(2:6).eq.'GAPFA') then
629  if(v(0,node1).eq.0.d0) then
630  v(0,node1)=v(0,node2)
631  else
632  v(0,node2)=v(0,node1)
633  endif
634  ichange=1
635  cycle
636  endif
637 !
638  nodem=kon(indexe+2)
639 !
640  if(v(1,nodem).eq.0.d0) then
641 !
642 ! direction of mass flow unknown in the element
643 !
644  cycle
645  elseif(v(1,nodem).gt.0.d0) then
646 !
647 ! positive mass flow (i.e. going from node1 to node2)
648 !
649  if(v(0,node1).eq.0.d0) cycle
650 !
651 ! propagating the temperature to node2
652 !
653  v(0,node2)=v(0,node1)
654  ichange=1
655  cycle
656  else
657 !
658 ! negative mass flow (i.e. going from node2 to node1)
659 !
660  if(v(0,node2).eq.0.d0) cycle
661 !
662 ! propagating the temperature to node1
663 !
664  v(0,node1)=v(0,node2)
665  ichange=1
666  cycle
667  endif
668  enddo
669 c write(*,*) 'preinitialnet '
670 c do i=1,ntg
671 c write(*,'(i10,3(1x,e11.4))') itg(i),(v(j,itg(i)),j=0,2)
672 c enddo
673  if(ichange.eq.0) exit
674  enddo
675 !
676 ! set of mass flow of +-1.d-30 to zero
677 !
678  do i=1,nflow
679  nelem=ieg(i)
680  indexe=ipkon(nelem)
681  nodem=kon(indexe+2)
682  if(dabs(v(1,nodem)).eq.1.d-30) v(1,nodem)=0.d0
683  enddo
684 !
685 ! check the pressures: set pressures to zero (i.e. no initial condtion)
686 ! which lie in between the neighboring pressures => Laplace method
687 ! is applied in initialnet
688 !
689  loop: do i=1,ntg
690  node=itg(i)
691 !
692 ! neighboring elements (excluding nodes which do not belong to
693 ! any network element or just to one network element (middle nodes)
694 !
695  index=iponoel(node)
696  if((index.eq.0).or.(inoel(2,index).eq.0)) cycle
697 !
698  imin=node
699  imax=node
700  xmin=v(2,node)
701  xmax=v(2,node)
702 !
703  do
704  nelem=inoel(1,index)
705  if(lakon(nelem)(2:3).eq.'VO') cycle loop
706  indexe=ipkon(nelem)
707 !
708 ! neighboring vertex node
709 !
710  if(kon(indexe+1).ne.node) then
711  neighbor=kon(indexe+1)
712  else
713  neighbor=kon(indexe+3)
714  endif
715  if(neighbor.eq.0) cycle loop
716 !
717 ! check its value
718 !
719  if(dabs(v(2,neighbor)).lt.xmin) then
720  xmin=dabs(v(2,neighbor))
721  imin=neighbor
722  elseif(dabs(v(2,neighbor)).gt.xmax) then
723  xmax=dabs(v(2,neighbor))
724  imax=neighbor
725  endif
726 !
727  index=inoel(2,index)
728  if(index.eq.0) exit
729  enddo
730 !
731 ! if value lies in between the neighboring values => assign a
732 ! negative sign as marker
733 !
734  if((imin.ne.node).and.(imax.ne.node)) then
735  v(2,node)=-v(2,node)
736  endif
737  enddo loop
738 !
739 ! set marked values to zero => Laplace equation will be used
740 ! in initialnet.f
741 !
742  do i=1,ntg
743  node=itg(i)
744  index=iponoel(node)
745  if((index.eq.0).or.(inoel(2,index).eq.0)) cycle
746 c if(v(2,node).lt.0.d0) v(2,node)=0.d0
747  if(v(2,node).lt.0.d0) v(2,node)=-v(2,node)
748  enddo
749 c write(*,*) 'preinitialnet end '
750 c do i=1,ntg
751 c write(*,'(i10,3(1x,e11.4))') itg(i),(v(j,itg(i)),j=0,2)
752 c enddo
753 !
754 ! same for temperatures?
755 !
756  return
networkneighbor
subroutine networkneighbor(nelem, node, nelemnei, nodenei, ibranch, iponoel, inoel, ipkon, kon)
Definition: networkneighbor.f:21
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)