genmap.f 119 KB
Newer Older
Paul Fischer's avatar
Paul Fischer committed
1
c-----------------------------------------------------------------------
Stefan's avatar
Stefan committed
2
3

   VERSION WITH SIMPLE BISECTION AND GROWTH OF CONN SET, COMBO (klh 07/09)
4
c
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
c   genmap.f performs the following operations:
c   
c   1)  Read a nek5k .rea or (.rea/.re2) file describing problem geometry - 2D or 3D.
c   
c   2)  Construct a standard finite element mesh from this file, i.e.,
c   
c         a. a list of coordinates, X1, X2,...,X_npts,  where X_i = (xi,yi,zi) in 3D
c            or X_i = (xi,yi) in 2D
c   
c   
c         b. a list of E cells, each cell comprising 2^d integers, d=2 or 3, where
c            E is the number of spectral elements specified in the .rea file.
c            Each integer in a cell points to one of the coordinate X_i above.
c   
c   3)  Generate the graph Laplacian, G, for the _dual_ graph corresponding to the 
c       FE mesh.   Each node in the dual graph is an element.  Each edge in the
c       graph is an element G_ij, for any element i connected to j.  
c   
c          . G_ij = 1  if element i shares a vertex with element j
c          . G_ij = 2  if element i shares an edge with element j
c          . G_ij = 4  if element i shares an face with element j
c   
c          . G_ii = - \sum j^=i G_ij
c      
c   
c   4)  Partition this graph using the Fiedler vector generated by a restarted
c       Lanczos iteration.
c   
c   5)  Recur on (3)-(4) in a depth-first traversal of the subgraphs until
c       P branches have been identified, where P = floor(log_2 E).
c   
c   
c   6)  Output the results
c   
c-----------------------------------------------------------------------
Paul Fischer's avatar
Paul Fischer committed
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
c
c  genmap() uses the symmetric vertex ordering
c
c      3 ----- 4
c      |       |
c      |       |
c      |       |
c      1 ----- 2
c
c  which can be extended to an arbitrary number of space dimensions
c  (like 3).
c
c-----------------------------------------------------------------------
      program genmap

c     read nekton .rea file and make a .map file

Stefan K's avatar
Stefan K committed
57
#     include "SIZE"
Paul Fischer's avatar
Paul Fischer committed
58
59
60
61
62
63

      parameter(lpts=8*lelm)
      common /carrayi/ cell (lpts) , pmap (lpts)
     $               , order(lpts) , elist(lpts)
      common /carrayw/ w1   (lpts) , w2   (lpts)
     $               , w3   (lpts) , w4   (lpts)
64
     $               , w5   (lpts)
Paul Fischer's avatar
Paul Fischer committed
65

66
67
68
      real        w14(4*lpts)
      equivalence (w1,w14)

Paul Fischer's avatar
Paul Fischer committed
69
70
71
      common /arrayr/  dx(4*lpts)

      common /carrayr/ bc(5*6*lelm)
72
      real*8 bc
Paul Fischer's avatar
Paul Fischer committed
73
74
75
      common /carrayc/ cbc(6,lelm)
      character*3      cbc

76
      logical ifconn,is_connected,face_conn
Paul Fischer's avatar
Paul Fischer committed
77
78

      integer     cell,pmap,order,elist,w1,w2,w3,w4,depth
79
      integer     e1,c1,f1
Paul Fischer's avatar
Paul Fischer committed
80

81
82
83
84
85
86
87
      wdsize=4
      eps=1.0e-12
      oneeps = 1.0+eps
      if (oneeps.ne.1.0) then
         wdsize=8
      endif
 
88
      call makemesh  (cell,nelv,nelt,irnk,dx,cbc,bc,ndim,w14)
Paul Fischer's avatar
Paul Fischer committed
89
c                                    irnk is # unique points
Paul Fischer's avatar
Paul Fischer committed
90
91
92

      nfc = 2*ndim
      nv  = 2**ndim
93
      mo  = 0   !max order initally zero.
Paul Fischer's avatar
Paul Fischer committed
94

95
96
      nic = lpts
      njc = lpts
97
98
99
100
101
102
103
104
105
106
      if (ndim.eq.3) then
         face_conn = .false.
         do i =1,5
            if (.not.face_conn) then
               call face_chk(face_conn,cell,nv,nelt,nic,njc,w1,w2,w3)
            else
               goto 15
            endif
         enddo 
         write(6,*) "WARNING:Missing Face Connection Not Resolved"
107
         call exitt(1) 
108
109
      endif
  15  continue
Paul Fischer's avatar
Paul Fischer committed
110
111

      call izero (order,irnk)
112
113
114
115
116
117

c     3/3/2010:  NO MORE SPECIAL OUTFLOW TREATMENT (fails w/ conj. ht.  transfer)

c     Determine number of outflow points and order them last
c     call set_outflow
c    $     (no,order,mo,cell,nv,nelv,irnk,cbc,nfc,w1,i0,i1,w2)
Paul Fischer's avatar
Paul Fischer committed
118

Paul Fischer's avatar
Paul Fischer committed
119
c     Find all periodic connections, based on cbc info.
Paul Fischer's avatar
Paul Fischer committed
120
      call periodic_vtx
121
     $               (cell,nv,nelt,irnk,dx,ndim,cbc,bc,nfc,w14,w5)
122
      
123

Paul Fischer's avatar
Paul Fischer committed
124
c     Recursive bisection of element graph; reverse-order interface points
Paul Fischer's avatar
Paul Fischer committed
125
      call rec_bisect (elist,pmap,order,mo,cell,nv,nelv,ndim
Paul Fischer's avatar
Paul Fischer committed
126
     $                                           ,w1,w2,w3,w4,w5)
Paul Fischer's avatar
Paul Fischer committed
127

Paul Fischer's avatar
Paul Fischer committed
128
c     Clean up 
Paul Fischer's avatar
Paul Fischer committed
129
130
      call isort     (elist,w1,nelv)
      call iswap_ip  (pmap ,w1,nelv)
Paul Fischer's avatar
Paul Fischer committed
131

Paul Fischer's avatar
Paul Fischer committed
132
133

      if (nelt.gt.nelv) then
134
135
136
137
138
         nels = nelt-nelv ! number of elements in solid
         e1 = 1 + nelv
         c1 = 1 + nelv*nv

         write(6,*) 'Starting rec_bisect2: nels = ',nels
Stefan Kerkemeier's avatar
Stefan Kerkemeier committed
139
         write(6,*)  
140
         write(6,*) 
Stefan Kerkemeier's avatar
Stefan Kerkemeier committed
141
     $     'NOTE: rec_bisect works only when solid is contiguous',
142
143
144
145
146
147
148
149
     $     'We could fix this with a connected graph test.'
         write(6,*) 

         call rec_bisect 
     $        (elist,pmap(e1),order,mo,cell(c1),nv,nels,ndim
     $                                           ,w1,w2,w3,w4,w5)
         call isort     (elist   ,w1,nels)
         call iswap_ip  (pmap(e1),w1,nels)
150
151
152

c         call maptest   (pmap,nelv,nelt,'map test A',w1,w2)

Paul Fischer's avatar
Paul Fischer committed
153
154
155
      endif

      npts = nv*nelt
Paul Fischer's avatar
Paul Fischer committed
156
      call iranku       (cell,nrnk,npts,w1)
Paul Fischer's avatar
Paul Fischer committed
157
      call self_chk     (cell,nv,nelt,1)    ! check for not self-ptg.
Paul Fischer's avatar
Paul Fischer committed
158
159

      call fill_order   (order,mo,cell,nv,nelt)
Paul Fischer's avatar
Paul Fischer committed
160
      call assign_order (cell,nv,nelt,order)
Paul Fischer's avatar
Paul Fischer committed
161
162

      call iranku       (cell,nrnk,npts,w1) ! make cell numbering contiguous
Paul Fischer's avatar
Paul Fischer committed
163
      call self_chk     (cell,nv,nelt,2)    ! check for not self-ptg.
164
      if (nelv.eq.nelt) call reverse_p (pmap,nelt) ! lightly load node 0
165
c      call maptest      (pmap,nelv,nelt,'map test B',w1,w2)
Paul Fischer's avatar
Paul Fischer committed
166
167

c     Output to .map file:
168
169
c     noutflow    = no    ! for now - no outflow bcs
      noutflow    = 0     ! for now - no outflow bcs  (handled in nek?)
170
171
      call out_mapfile (pmap,nelt,cell,nv,nrnk,noutflow)
c      call maptest     (pmap,nelv,nelt,'map test C',w1,w2)
Paul Fischer's avatar
Paul Fischer committed
172

173
174
      call mult_chk(dx,ndim,nv,nelt,cell,nrnk)

Paul Fischer's avatar
Paul Fischer committed
175
176
c     call out_geofile (dx,ndim,nv,nelt,pmap,39)
c     call out_geofile2(dx,ndim,nv,nelt,cell,nrnk)
Paul Fischer's avatar
Paul Fischer committed
177

Paul Fischer's avatar
Paul Fischer committed
178
179
180
181
182
c     call outmati(pmap,13,9,'pmap  ',nelt,1)
c     open(unit=22,file='p.dat')
c     write(22,1) (pmap(k),k=1,nelt)
c   1 format(i9)
c     close(unit=22)
Paul Fischer's avatar
Paul Fischer committed
183
184
185

      end
c-----------------------------------------------------------------------
186
      subroutine makemesh(cell,nelv,nelt,irnk,dx,cbc,bc,ndim,wk)
Paul Fischer's avatar
Paul Fischer committed
187

188
c     read nekton .rea file and make a mesh
Paul Fischer's avatar
Paul Fischer committed
189

Stefan K's avatar
Stefan K committed
190
#     include "SIZE" 
Stefan Kerkemeier's avatar
Stefan Kerkemeier committed
191

192
      integer      cell(1)
Paul Fischer's avatar
Paul Fischer committed
193
      character*3  cbc (6,1)
194
      real*8       bc  (5,6,1)
Paul Fischer's avatar
Paul Fischer committed
195
      real         dx  (1)
196
      real         wk  (1)
Paul Fischer's avatar
Paul Fischer committed
197
198
199
      integer e,f

      character*3 cbt(6)
200
      real*8      bt(5,6)
Paul Fischer's avatar
Paul Fischer committed
201
202
203
204

      parameter(lpts=8*lelm)

      common /arrayi/ i_n(lpts) , j_n(4*lpts)
205
     $                          , j_o(4*lpts)
Paul Fischer's avatar
Paul Fischer committed
206

207
208
      logical ifbinary,ifbswap
      integer buf(30)
Paul Fischer's avatar
Paul Fischer committed
209
210
211
212

      integer eface(6)  ! return Nekton preprocessor face ordering
      save    eface
      data    eface / 4 , 2 , 1 , 3 , 5 , 6 /
213
         
214
      io = 10
215
216
      ifma2 = .false.

217
218
219
      call getreafile('Input .rea / .re2 name:$',ifbinary,io,ierr)
      if (ierr.gt.0) then 
         write(6,'(A)') 'Error no .rea / .re2 file found!'
220
221
222
c         call linearmsh(cell,nelv,nelt,ndim)
c         return
         call exitt(1) 
223
      endif
224
225

      write(6,'(A)') 'Input mesh tolerance (default 0.2):'
226
227
      write(6,'(A,A)') 'NOTE: smaller is better, but generous is more ',
     &                 'forgiving for bad meshes.'
228

229
230
231
232
233
234
235
236
      read(5,'(f7.2)') qin
      if(qin.gt.0) then
        q = qin
      else
        write(6,'(A,2f7.2)') ' using default value'
        q = 0.2
      endif
  
Paul Fischer's avatar
Paul Fischer committed
237
      call cscan_dxyz (dx,nelt,nelv,ndim,ifbinary,ifbswap)
Paul Fischer's avatar
Paul Fischer committed
238

239
      ierr = 0
240
      if (ifbinary) then
241
         ifma2 = .true.
242
         ! skip curved side data
243
244
         if(wdsizi.eq.8) then 
           call byte_read(rcurve,2,ierr)
245
           if (ifbswap) call byte_reverse8(rcurve,2,ierr)
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
           ncurve = rcurve
           if(ierr.ne.0) call exitti
     $         ('Error reading ncurve in makemesh ',ierr)
           do k = 1,ncurve
              call byte_read(buf,16,ierr)
              if(ierr.ne.0) call exitti
     $         ('Error reading curve data in makemesh ',ierr)
           enddo
         else
           call byte_read(ncurve,1,ierr)
           if (ifbswap) call byte_reverse(ncurve,1,ierr)
           if(ierr.ne.0) call exitti
     $         ('Error reading ncurve in makemesh ',ierr)
           do k = 1,ncurve
              call byte_read(buf,8,ierr)
              if(ierr.ne.0) call exitti
     $         ('Error reading curve data in makemesh ',ierr)
           enddo
         endif
Paul Fischer's avatar
Paul Fischer committed
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
c        For current version of genmap, only need the fluid bcs.
c        Later, for more complex MHD configs, we'll need fluid + induct.

c        Also, if we start to have fluid/thermal domains with differing
c        periodic bc topologies, then we'll need something completely
c        different (two distinct numberings).
c
c        In fact, we need the thermal bcs because of periodicity... 
c        This is not pretty --- there are many cases to be considered.
c
c        For now, we default to solid topology if nelt > nelv
c


         call rd_bc_bin(cbc,bc,nelv,nelt,ifbswap)

281
282
         call byte_close(ierr)
         if(ierr.ne.0) call exitti
283
     $       ('Error closing file in makemesh ',ierr)
284
      else
Paul Fischer's avatar
Paul Fischer committed
285
         call cscan_bcs   (cbc,bc,nelv,nelt,ndim)
286
      endif
Paul Fischer's avatar
Paul Fischer committed
287
c     call outbc(cbc,bc,nelt,ndim,' CBC 1')
288
289
      close (unit=10)  ! close .rea file

Paul Fischer's avatar
Paul Fischer committed
290
291
292
      nface = 2*ndim
      do e=1,nelt !  SWAP TO PREPROCESSOR NOTATION
         call chcopy(cbt,cbc(1,e)  ,nface*3)
293
         call copy8 ( bt, bc(1,1,e),nface*5)
Paul Fischer's avatar
Paul Fischer committed
294
         do f=1,nface
295
            call copy8 ( bc(1,f,e), bt(1,eface(f)),5)
Paul Fischer's avatar
Paul Fischer committed
296
297
298
299
300
            call chcopy(cbc(  f,e),cbt(  eface(f)),3)
         enddo
      enddo
c     call outbc(cbc,bc,nelt,ndim,' CBC 2')

Paul Fischer's avatar
Paul Fischer committed
301
302

c     Compress vertices based on coordinates
303
      call unique_vertex2(cell,dx,ndim,nelt,q,i_n,j_n,j_o,wk)
Paul Fischer's avatar
Paul Fischer committed
304
305

      nv   = 2**ndim
Paul Fischer's avatar
Paul Fischer committed
306
      npts = nelt*nv
307

Paul Fischer's avatar
Paul Fischer committed
308
      call iranku    (cell,irnk,npts,i_n)
309
      call self_chk  (cell,nv,nelt,32)       ! check for not self-ptg.
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
      return
      end
c-----------------------------------------------------------------------
      subroutine linearmsh(cell,nelv,nelt,ndim)
      
      integer      cell(2,1)
      integer e

      ndim = 1

      write(6,*) 'Input number of linear finite elements:'
      read (5,*) nelt

      do e = 1,nelt
         cell(1,e) = e
         cell(2,e) = e+1
      enddo

      nelv=nelt


      
Paul Fischer's avatar
Paul Fischer committed
332
333
334
      return
      end
c-----------------------------------------------------------------------
335
336
      subroutine exitti(name,ie)
      character*40 name
Stefan's avatar
Stefan committed
337
338
      write(6,*) name, ie
      stop
339
      end
340
c-----------------------------------------------------------------------
Paul Fischer's avatar
Paul Fischer committed
341
      subroutine exitt(ie)
Stefan's avatar
Stefan committed
342
343
      write(6,*) 'exit status', ie
      stop
Paul Fischer's avatar
Paul Fischer committed
344
      end
345
c-----------------------------------------------------------------------
Paul Fischer's avatar
Paul Fischer committed
346
      subroutine cscan_dxyz (dx,nelt,nelv,ndim,ifbinary,ifbswap)
Paul Fischer's avatar
Paul Fischer committed
347
348
349
c
c     Scan for xyz data, read it, and set characteristic length, d2
c
350

Stefan K's avatar
Stefan K committed
351
#     include "SIZE" 
Stefan Kerkemeier's avatar
Stefan Kerkemeier committed
352
     
Paul Fischer's avatar
Paul Fischer committed
353
354
355
356
      character*80 string
c
      real dx(1)
      real x(8),y(8),z(8)
357
      integer e,buf(50)
Paul Fischer's avatar
Paul Fischer committed
358
359
360
361
362

      integer h2s(8) ! hypercube to strange ordering
      save    h2s
      data    h2s / 1,2,4,3,5,6,8,7 /

363
364
365
366
      logical ifbinary,ifbswap

      ifbswap  = .false.

367
368
369
370
371
372
      write(6,*) 'reading mesh data ...'

      if (.not. ifbinary) then
         call cscan(string,'MESH DATA',9)
         read (10,*) nelt,ndim,nelv
      endif
373
       
374
      if (nelt.lt.0 .or. ifbinary) then
375
         ifbinary = .true.
376
         call open_bin_file(ifbswap,nelgtr,ndimr,nelgvr,wdsizi)
377
378
379
380
381
         if(wdsize.eq.4.and.wdsizi.eq.8) then
             write(6,*) "Double Precision .rea not supported ",
     $                  "in Single Precision mode, compile with -r8"
             call exitt(wdsize)
         endif
382
         nelt = nelgtr
383
         ndim = ndimr
384
         nelv = nelgvr
385
         nwds = (1 + ndim*(2**ndim))*(wdsizi/4) ! group + 2x4 for 2d, 3x8 for 3d
386
387
      endif

388
c      write(6,*) nelt,ndim,nelv,ifbinary, ' nelt,ndim,nelv,ifre2 '
389

Stefan K's avatar
Stefan K committed
390
391
392
      if (nelt.gt.lelm) then
        write(6,*) 'Abort: number of elements too large', nelt
        write(6,*) 'change MAXNEL and recompile' 
393
        call exitt(1)
394
      endif
Paul Fischer's avatar
Paul Fischer committed
395
396
397

      b = 1.e22
      l = 1
398
399

      ierr = 0
Paul Fischer's avatar
Paul Fischer committed
400
      if (ndim.eq.3) then
Paul Fischer's avatar
Paul Fischer committed
401
         do e=1,nelt
402
            if(ifbinary) then
403
404
              call byte_read(buf,nwds,ierr)
              if(ierr.ne.0) goto 100
405
              call buf_to_xyz(buf,x,y,z,e,ifbswap,ndim,wdsizi)
406
407
408
409
410
411
412
413
414
            else 
              read (10,80) string
              read (10,*)   (x(k),k=1,4)
              read (10,*)   (y(k),k=1,4)
              read (10,*)   (z(k),k=1,4)
              read (10,*)   (x(k),k=5,8)
              read (10,*)   (y(k),k=5,8)
              read (10,*)   (z(k),k=5,8)
            endif
415
416
417
418
c           write(6,*) e
c           do ii=1,8
c          write(6,*) x(ii),y(ii),z(ii)
c           enddo
Paul Fischer's avatar
Paul Fischer committed
419
420
421
422
423
424
425
426
427
            do k=1,8
               dx(l+0) = b
               dx(l+1) = x(h2s(k))
               dx(l+2) = y(h2s(k))
               dx(l+3) = z(h2s(k))
               l = l + (ndim+1)
            enddo
         enddo
      else
Paul Fischer's avatar
Paul Fischer committed
428
         do e=1,nelt
429
            if(ifbinary) then
430
431
              call byte_read(buf,nwds,ierr)
              if(ierr.ne.0) goto 100
432
              call buf_to_xyz(buf,x,y,z,e,ifbswap,ndim,wdsizi)
433
434
435
436
437
            else
              read (10,80) string
              read (10,*)   (x(k),k=1,4)
              read (10,*)   (y(k),k=1,4)
            endif
Paul Fischer's avatar
Paul Fischer committed
438
439
440
441
442
443
444
445
446
            do k=1,4
               dx(l+0) = b
               dx(l+1) = x(h2s(k))
               dx(l+2) = y(h2s(k))
               l = l + (ndim+1)
            enddo
         enddo
      endif
   80 format(a80)
447
 
Paul Fischer's avatar
Paul Fischer committed
448
      nvrt = 2**ndim
Paul Fischer's avatar
Paul Fischer committed
449
      call set_d2(dx,nvrt,nelt,ndim)
450
451
452
453
 
      return
 100  write(6,*) "Error reading xyz byte data in scan_dxyz. Abort"
      call exitt(ierr)
Paul Fischer's avatar
Paul Fischer committed
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
      return
      end
c-----------------------------------------------------------------------
      subroutine set_d2(dx,nvrt,nel,ndim)
      real dx(0:ndim,nvrt,nel)

      integer neigh(3,8),e
      save    neigh
      data    neigh / 2,3,5 , 1,4,6 , 1,4,7 , 2,3,8   ! symm. ordering
     $              , 1,6,7 , 2,5,8 , 3,5,8 , 4,6,7 /

      b = 1.e22

      do e = 1,nel
      do i = 1,nvrt

         dx(0,i,e) = b

         do k = 1,ndim

            n   = neigh(k,i)
            d2l = 0.

            do j=1,ndim
               d2l = d2l + (dx(j,n,e)-dx(j,i,e))**2
            enddo

            dx(0,i,e) = min (dx(0,i,e),d2l)
c           write(6,9) i,e,' dx',(dx(j,i,e),j=0,ndim)
c  9        format(2i6,a3,1p4e11.3)

         enddo

      enddo
      enddo
c     call exitt(5)

      return
      end
c-----------------------------------------------------------------------
Paul Fischer's avatar
Paul Fischer committed
494
      subroutine cscan_bcs  (cbc,bc,nelv,nelt,ndim)
Paul Fischer's avatar
Paul Fischer committed
495
496
497
c
c     Scan for cbc data and read it
c
Paul Fischer's avatar
Paul Fischer committed
498
      character*3  cbc(6,nelt)
499
      real*8       bc (5,6,nelt)
Paul Fischer's avatar
Paul Fischer committed
500
501
502
503
504
505
      character*80 string

      integer e

      call cscan(string,'BOUNDARY',8) ! for now, fluid only

Paul Fischer's avatar
Paul Fischer committed
506
507
      npass = 1
      if (nelt.gt.nelv) npass = 2
Paul Fischer's avatar
Paul Fischer committed
508

Paul Fischer's avatar
Paul Fischer committed
509
510
511
512
      do kpass = 1,npass

         read (10,80) string
   80    format(a80)
Paul Fischer's avatar
Paul Fischer committed
513

Paul Fischer's avatar
Paul Fischer committed
514
515
516
517
518
519
520
         ifield = kpass
         if (indx1(string,'NO ',3).ne.0) then
            ifield = ifield+1
            call cscan(string,'BOUNDARY',8) ! then, temp only
         endif

         nel=nelv
521
         if (kpass.eq.2) nel=nelt
Paul Fischer's avatar
Paul Fischer committed
522
523
524
         call rd_bc(cbc,bc,nel,ndim,ifield,10)

      enddo
Paul Fischer's avatar
Paul Fischer committed
525
526
527
528

      return
      end
c-----------------------------------------------------------------------
Paul Fischer's avatar
Paul Fischer committed
529
      subroutine rd_bc(cbc,bc,nel,ndim,ifield,io)
Paul Fischer's avatar
Paul Fischer committed
530
531
532

c     .Read Boundary Conditions (and connectivity data)

Paul Fischer's avatar
Paul Fischer committed
533
      character*3 cbc(6,nel)
534
      real*8      bc(5,6,nel)
Paul Fischer's avatar
Paul Fischer committed
535
      integer e,f
536
537
538

c     write(6,*) 'inside rd_bc: ',nel,ndim,ifield,io

Paul Fischer's avatar
Paul Fischer committed
539
      nbcrea = 5
Paul Fischer's avatar
Paul Fischer committed
540
      nface  = 2*ndim
Paul Fischer's avatar
Paul Fischer committed
541
      do e=1,nel
Paul Fischer's avatar
Paul Fischer committed
542
      do f=1,nface
Paul Fischer's avatar
Paul Fischer committed
543
         if (nel.lt.1000) then
544
            read(io,50,err=510,end=600)    
Paul Fischer's avatar
Paul Fischer committed
545
546
547
     $      chtemp,
     $      cbc(f,e),id1,id2,
     $      (bc(ii,f,e),ii=1,nbcrea)
548
   50       format(a1,a3,2i3,5g14.6)
549
         elseif (nel.lt.100 000) then
550
            read(io,51,err=520,end=600)    
Paul Fischer's avatar
Paul Fischer committed
551
552
553
     $      chtemp,
     $      cbc(f,e),id1,id2,
     $      (bc(ii,f,e),ii=1,nbcrea)
554
   51       format(a1,a3,i5,i1,5g14.6)
555
556
557
         elseif (nel.lt.1 000 000) then
            read(io,52,err=530,end=600)    
     $      cbc(f,e),id1,(bc(ii,f,e),ii=1,nbcrea)
558
   52       format(1x,a3,i6,5g14.6)
Paul Fischer's avatar
Paul Fischer committed
559
         else
560
            read(io,53,err=540,end=600)    
Paul Fischer's avatar
Paul Fischer committed
561
     $      cbc(f,e),id1,(bc(ii,f,e),ii=1,nbcrea)
562
   53       format(1x,a3,i12,5g18.11)
Paul Fischer's avatar
Paul Fischer committed
563
         endif
564
c        write(6,*) e,f,' ',cbc(f,e),' BC IN?',nel,ifield
Paul Fischer's avatar
Paul Fischer committed
565
566
567
568
569
570
571
      enddo
      enddo

      return
C
C     Error handling:
C
572
  500 format(2x,'ERROR: error reading ',i4,2i12,/,
Stefan Kerkemeier's avatar
Stefan Kerkemeier committed
573
     $       2x,'aborting ',a3,' in routine rdbdry.')
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
  510 write(6,500) ifield,e,nel,'510'
      call exitt(ifield)
      return

  520 write(6,500) ifield,e,nel,'520'
      call exitt(ifield)
      return

  530 write(6,500) ifield,e,nel,'530'
      call exitt(ifield)
      return

  540 write(6,500) ifield,e,nel,'540'
      call exitt(ifield)
      return

Paul Fischer's avatar
Paul Fischer committed
590
      return
591
592
593

  600 continue
      write(6,601) ifield,e,nel
594
  601 FORMAT(2X,'ERROR: end of file',i4,2i12,/,
Stefan Kerkemeier's avatar
Stefan Kerkemeier committed
595
     $       2X,'ABORTING 600 IN ROUTINE RDBDRY.')
596
597
598
      call exitt(ifield)
      return

Paul Fischer's avatar
Paul Fischer committed
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
      end
c-----------------------------------------------------------------------
      subroutine blank(s,n)
      character*1 s(1)
      do i=1,n
        s(i)=' '
      enddo
      return
      end
c-----------------------------------------------------------------------
      function ltrunc(s,n)
      character*1 s(1)
      ltrunc = 0
      do j=n,1,-1
         if (s(j).ne.' ') then
            ltrunc = j 
            return
         endif
      enddo
      return
      end
c-----------------------------------------------------------------------
      integer function indx1(s1,s2,l2)
      character*80 s1,s2
C
      n1=80-l2+1
      indx1=0
      if (n1.lt.1) return
C
      do 300 i=1,n1
         i2=i+l2-1
         if (s1(i:i2).eq.s2(1:l2)) then
            indx1=i
            return
         endif
300   continue
c
      return
      end
c-----------------------------------------------------------------------
      subroutine cscan(sout,key,nk)
      character*80 sout,key
      character*80 string
      character*1  string1(80)
      equivalence (string1,string)
c
      do i=1,99999999      
         call blank(string,80)
         read (10,80,end=100,err=100) string
         call chcopy(sout,string,80)
c        write (6,*) string
         if (indx1(string,key,nk).ne.0) return
      enddo
  100 continue
c
   80 format(a80)
      return
      end
c
c-----------------------------------------------------------------------
      subroutine readwrite(sout,key,nk)
      character*80 sout,key
      character*80 string
      character*1  string1(80)
      equivalence (string1,string)
c
      do i=1,90000      
         call blank(string,80)
         read (10,80,end=100,err=100) string
         len = ltrunc(string,80)
         write(11,81) (string1(k),k=1,len)
         if (indx1(string,key,nk).ne.0) return
      enddo
  100 continue
c
   80 format(a80)
   81 format(80a1)
      return
      end
c
c-----------------------------------------------------------------------
      subroutine readwrite2(sout,key1,nk1,key2,nk2)
      character*80 sout,key1,key2
      character*80 string
      character*1  string1(80)
      equivalence (string1,string)
c
      do i=1,90000      
         call blank(string,80)
         read (10,80,end=100,err=100) string
         len = ltrunc(string,80)
         write(11,81) (string1(k),k=1,len)
c        write(6 ,81) (string1(k),k=1,len)
         if (indx1(string,key1,nk1).ne.0) return
         if (indx1(string,key2,nk2).ne.0) return
      enddo
  100 continue
   80 format(a80)
   81 format(80a1)
      return
      end
c
c-----------------------------------------------------------------------
      integer function log2(k)
      rk=(k)
      rlog=log10(rk)
      rlog2=log10(2.0)
      rlog=rlog/rlog2  + 1.e-6  !  + 0.5  ! don't round up!
      log2=int(rlog)
      return
      end
710
711
712
713
714
715
716
717
718
719
c-----------------------------------------------------------------------
      function iglmin(a,n)
      integer a(1),tmin
      tmin=999999999
      do 100 i=1,n
         tmin=min(tmin,a(i))
  100 continue
      iglmin=tmin
      return
      end
720
721
722
723
724
725
726
727
728
729
c-----------------------------------------------------------------------
      function ivlmax(a,n)
      integer a(1),tmax
      tmax=-999999999
      do 100 i=1,n
         tmax=max(tmax,a(i))
  100 continue
      ivlmax=tmax
      return
      end
Paul Fischer's avatar
Paul Fischer committed
730
731
732
c-----------------------------------------------------------------------
      function iglmax(a,n)
      integer a(1),tmax
733
      tmax=-999999999
Paul Fischer's avatar
Paul Fischer committed
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
      do 100 i=1,n
         tmax=max(tmax,a(i))
  100 continue
      iglmax=tmax
      return
      end
c-----------------------------------------------------------------------
      function glmax(a,n)
      real a(1)
      tmax=-99.0E20
      do 100 i=1,n
         tmax=max(tmax,a(i))
  100 continue
      glmax=tmax
      return
      end
c-----------------------------------------------------------------------
      function glmin(a,n)
      real a(1)
      tmin=99.0E20
      do 100 i=1,n
         tmin=min(tmin,a(i))
  100 continue
      glmin = tmin
      return
      end
c-----------------------------------------------------------------------
      subroutine icopy(x,y,n)
      integer x(1),y(1)
      do i=1,n
         x(i) = y(i)
      enddo
      return
      end
c-----------------------------------------------------------------------
      subroutine cmult2(x,y,c,n)
      real x(1),y(1)
      do i=1,n
         x(i) = c*y(i)
      enddo
      return
      end
c-----------------------------------------------------------------------
      subroutine copy(x,y,n)
      real x(1),y(1)
      do i=1,n
         x(i) = y(i)
      enddo
      return
      end
784
785
786
787
788
789
790
791
c-----------------------------------------------------------------------
      subroutine copy8(x,y,n)
      real*8 x(1),y(1)
      do i=1,n
         x(i) = y(i)
      enddo
      return
      end
Paul Fischer's avatar
Paul Fischer committed
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
c-----------------------------------------------------------------------
      subroutine izero(x,n)
      integer x(1)
      do i=1,n
         x(i) = 0
      enddo
      return
      end
c-----------------------------------------------------------------------
      subroutine cfill(x,c,n)
      real x(1)
      do i=1,n
         x(i) = c
      enddo
      return
      end
c-----------------------------------------------------------------------
      subroutine rint(x,n)
      real x(1)
      do i=1,n
         x(i) = i
      enddo
      return
      end
c-----------------------------------------------------------------------
      subroutine jjnt(x,n)
      integer x(1)
      do i=1,n
         x(i) = i
      enddo
      return
      end
c-----------------------------------------------------------------------
      subroutine ifill(x,c,n)
      integer x(1),c
      do i=1,n
         x(i) = c
      enddo
      return
      end
c-----------------------------------------------------------------------
      subroutine icadd(x,c,n)
      integer x(1),c
      do i=1,n
         x(i) = x(i)+c
      enddo
      return
      end
c-----------------------------------------------------------------------
      subroutine cadd(x,c,n)
      real x(1)
      do i=1,n
         x(i) = x(i)+c
      enddo
      return
      end
c-----------------------------------------------------------------------
      subroutine cmult(x,c,n)
      real x(1)
      do i=1,n
         x(i) = x(i)*c
      enddo
      return
      end
c-----------------------------------------------------------------------
      subroutine chcopy(x,y,n)
      character*1 x(1),y(1)
      do i=1,n
         x(i) = y(i)
      enddo
      return
      end
c-----------------------------------------------------------------------
865
      subroutine getreafile(prompt,ifbinary,io,ierr)
Paul Fischer's avatar
Paul Fischer committed
866
c
867
868
      character*1 prompt(1)
      logical ifbinary
Paul Fischer's avatar
Paul Fischer committed
869
870
871
872
873
874
875
876
c
      common /sess/ session
      character*80 session

      character*80 file
      character*1  file1(80)
      equivalence (file1,file)

877
      ierr = 0
Stefan's avatar
Stefan committed
878
      ifbinary = .false.
879
880

c     Get file name
Paul Fischer's avatar
Paul Fischer committed
881
882
883
884
885
886
      len = indx1(prompt,'$',1) - 1
      write(6,81) (prompt(k),k=1,len)
   81 format(80a1)
      call blank(session,80)
      read(5,80) session
   80 format(a80)
887

888
889
      if (session.eq.'-1') then
         io = -1
890
         ierr = 1
891
892
893
894
         return
      else
         call chcopy(file,session,80)
         len = ltrunc(file,80)
895
896
897
898
899
900
901
902
903
904
         call chcopy(file1(len+1),'.rea',4)
         open(unit=io, file=file, status='old', iostat=ierr)

         if (ierr.gt.0) then
            call chcopy(file,session,80)
            len = ltrunc(file,80)
            call chcopy(file1(len+1),'.re2',4)
            inquire(file=file, exist=ifbinary)
            if(ifbinary) ierr = 0
         endif
905
      endif
906

Stefan's avatar
Stefan committed
907
908
      write(6,*) 'reading ', file

Paul Fischer's avatar
Paul Fischer committed
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
      return
      end
c-----------------------------------------------------------------------
      subroutine isort(a,ind,n)
C
C     Use Heap Sort (p 231 Num. Rec., 1st Ed.)
C
      integer a(1),ind(1)
      integer aa
C
      dO 10 j=1,n
         ind(j)=j
   10 continue
C
      if (n.le.1) return
      L=n/2+1
      ir=n
  100 continue
         if (l.gt.1) then
            l=l-1
            aa  = a  (l)
            ii  = ind(l)
         else
                 aa =   a(ir)
                 ii = ind(ir)
              a(ir) =   a( 1)
            ind(ir) = ind( 1)
            ir=ir-1
            if (ir.eq.1) then
                 a(1) = aa
               ind(1) = ii
               return
            endif
         endif
         i=l
         j=l+l
  200    continue
         if (j.le.ir) then
            if (j.lt.ir) then
               if ( a(j).lt.a(j+1) ) j=j+1
            endif
            if (aa.lt.a(j)) then
                 a(i) = a(j)
               ind(i) = ind(j)
               i=j
               j=j+j
            else
               j=ir+1
            endif
         GOTO 200
         endif
           a(i) = aa
         ind(i) = ii
      GOTO 100
      end
c-----------------------------------------------------------------------
      subroutine iswap_ip(x,p,n)
      integer x(1),xstart
      integer p(1)
c
c     In-place permutation: x' = x(p)
c
      do k=1,n
         if (p(k).gt.0) then   ! not swapped
            xstart     = x(k)
            loop_start = k
            last       = k
            do j=k,n
               next    = p(last)
               if (next.lt.0) then
                  write(6,*) 'Hey! iswap_ip problem.',j,k,n,next
                  call exitt(0)
               elseif (next.eq.loop_start) then
                  x(last) = xstart
                  p(last) = -p(last)
                  goto 10
               else
                  x(last) = x(next)
                  p(last) = -p(last)
                  last    = next
               endif
            enddo
   10       continue
         endif
      enddo
c
      do k=1,n
         p(k) = -p(k)
      enddo
      return
      end
c-----------------------------------------------------------------------
For faster browsing, not all history is shown. View entire blame