← Back to branch summary

~fluidity-core/fluidity/embedded_models

« back to all changes in this revision

Viewing changes to libadaptivity/adapt3d/src/mkfixd.F

  • Committer: Timothy Bond
  • Date: 2011-04-14 15:40:24 UTC
  • Revision ID: timothy.bond@imperial.ac.uk-20110414154024-116ci9gq6mwigmaw
Following the move from svn to bzr we change the nature of inclusion of these
four software libraries. Previously, they were included as svn externals and
pulled at latest version for trunk, pinned to specific versions for release
and stable trunk. Since bzr is less elegant at dealing with externals we have
made the decision to include the packages directly into the trunk instead.

At this import the versions are:

libadaptivity: r163
libvtkfortran: r67
libspud: r545
libmba2d: r28

Show diffs side-by-side

added added

removed removed

Lines of Context:
libadaptivity/adapt3d/src/mkfixd.F
 
1
C Copyright (C) 2006 Imperial College London and others.
 
2
C 
 
3
C Please see the AUTHORS file in the main source directory for a full list
 
4
C of copyright holders.
 
5
C 
 
6
C Adrian Umpleby
 
7
C Applied Modelling and Computation Group
 
8
C Department of Earth Science and Engineering
 
9
C Imperial College London
 
10
C 
 
11
C adrian@Imperial.ac.uk
 
12
C 
 
13
C This library is free software; you can redistribute it and/or
 
14
C modify it under the terms of the GNU Lesser General Public
 
15
C License as published by the Free Software Foundation; either
 
16
C version 2.1 of the License.
 
17
C 
 
18
C This library is distributed in the hope that it will be useful,
 
19
C but WITHOUT ANY WARRANTY; without even the implied warranty of
 
20
C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 
21
C Lesser General Public License for more details.
 
22
C 
 
23
C You should have received a copy of the GNU Lesser General Public
 
24
C License along with this library; if not, write to the Free Software
 
25
C Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307
 
26
C USA
 
27
#include "ewrite.h"
 
28
      SUBROUTINE MKFIXD( BIGLST, NODLST, NDNMBR,
 
29
     :                   NODX, NODY, NODZ, ORGMTX, SIZMTX,
 
30
     :                   PRDNDS, NPRDND,
 
31
     :                   NODLOC, SIZLOC, SUROWN, SIZOWN,
 
32
     :                   ELMREG, ENLBAS, ENLIST,
 
33
     :                   SNLBAS, SNLIST, SURFID,
 
34
     :                   NELM, NNOD, NSELM, SZENLS, SZSNLS, GEOM3D,
 
35
     :                   GATHER, SCATER, NGATH, NHALO,
 
36
     :                   ATOSEN, ATOREC, NPROC )
 
37
C-----------------------------------------------------------------------
 
38
C
 
39
C - This subroutine forms the element node list (ENLIST) and its base
 
40
C - pointer (ENLBAS) for the new adapted mesh, as well as returning the
 
41
C - element material/region id in ELMREG for each element. The nodal
 
42
C - positions for the adapted mesh is returned in NODX, NODY and NODZ,
 
43
C - and the interpolated error metrics are returned in ORGMTX if SIZMTX
 
44
C - is sent as 9*NNOD (otherwise it should be sent as 1).
 
45
C
 
46
C-----------------------------------------------------------------------
 
47
      use write_log
 
48
      IMPLICIT NONE
 
49
C
 
50
      INTEGER NELM, SZENLS, NNOD, SZSNLS, NSELM, SIZOWN, SIZLOC, SIZMTX,
 
51
     :        NHALO, NPROC, NGATH, NPRDND
 
52
C
 
53
      INTEGER ELMREG(NELM), ENLIST(SZENLS), ENLBAS(NELM+1), 
 
54
     :        SNLBAS(NSELM+1), SNLIST(SZSNLS), SUROWN(SIZOWN), 
 
55
     :        NODLOC(SIZLOC), GATHER(NGATH), SCATER(NHALO),
 
56
     :        ATOSEN(NPROC), ATOREC(NPROC), PRDNDS(NPRDND),
 
57
     :        SURFID(NSELM)
 
58
C
 
59
      REAL NODX(NNOD), NODY(NNOD), NODZ(NNOD), ORGMTX(SIZMTX)
 
60
C
 
61
      LOGICAL GEOM3D, NDNMBR
 
62
C
 
63
      INCLUDE 'blknew.i'
 
64
C
 
65
      INCLUDE 'blkbig.i'
 
66
C
 
67
      INCLUDE 'blkerr.i'
 
68
C
 
69
      REAL XSHIFT, YSHIFT, ZSHIFT, XSCALE, YSCALE, ZSCALE
 
70
      COMMON / BLKSCL / XSHIFT, YSHIFT, ZSHIFT, XSCALE, YSCALE, ZSCALE
 
71
C
 
72
      INTEGER IPTR, CONELM(4), EDGLST(6), INOD, IELM, IFLAG, IREG,
 
73
     :        INEW, NXT, LST, NODS(4), IPOS, I, J, K, SPOS, SELM,
 
74
     :        NODCON(4), NXTC, SUM, CNTGMY, CNTINT, CNTSPL, JNOD, KNOD,
 
75
     :        minelv, maxelv, minelr, maxela, i1, i2
 
76
C
 
77
      REAL X(4), Y(4), Z(4), TETVOL, avol, arad, asp, m(9), areas(4),
 
78
     :     minvol, minrad, minrvl, minvrd, maxasp, maxavl, maxard,
 
79
     :     minras, minvas, L(6), maxvol, maxvrd, maxvas
 
80
      integer mxsfid, nodcnt, nusfid
 
81
C      
 
82
      integer, dimension(:), allocatable:: srfcnt
 
83
      real, dimension(:), allocatable:: srfavx, srfavy, srfavz
 
84
      integer, parameter:: MAX_SFID_STATS=1000
 
85
C
 
86
      REAL*8 VOLSUM, ASPAVE
 
87
C
 
88
      LOGICAL GTNDFL, FACS(4), getflg
 
89
C
 
90
      IF( IERR .NE. 0 ) RETURN
 
91
C
 
92
c      print*,'NNOD,NELM,NSELM,GEOM3D: ',NNOD,NELM,NSELM,IS3DMS
 
93
      IF( GEOM3D .NEQV. IS3DMS ) THEN
 
94
        WRITE(0,*) '*** MKFIXD: Inconsistent GEOM3D!'
 
95
        IF( GEOM3D ) THEN
 
96
          WRITE(0,*) "   Didn't do 3D adapt before, but GEOM3D now true"
 
97
        ELSE
 
98
          WRITE(0,*) '   Did 3D adapt before, but GEOM3D now false'
 
99
        END IF
 
100
        IERR = -7835
 
101
        RETURN
 
102
      END IF
 
103
c      print*,'SIZMTX: ',SIZMTX
 
104
C
 
105
      CNTGMY = 0
 
106
      CNTSPL = 0
 
107
      CNTINT = 0
 
108
      IPTR = STTNOD
 
109
      INOD = 0
 
110
      JNOD = NNOD + 1
 
111
      NODCNT = 0
 
112
C
 
113
  10  IF( GTNDFL(BIGLST,NODLST,IPTR,1) ) THEN
 
114
C
 
115
         if( gtndfl(biglst,nodlst,iptr,2) ) then
 
116
            cntgmy = cntgmy + 1
 
117
         else if( gtndfl(biglst,nodlst,iptr,4) ) then
 
118
            cntspl = cntspl + 1
 
119
         else if( gtndfl(biglst,nodlst,iptr,3) ) then
 
120
            cntint = cntint + 1
 
121
         end if
 
122
c
 
123
         nodcnt = nodcnt + 1
 
124
         if(nodcnt.gt.nnod) then
 
125
            WRITE(0,*) '*** MKFIXD: Too many nodes!'
 
126
            WRITE(0,*) nodcnt,nnod,iptr,topnod
 
127
            ierr = -9797
 
128
            return
 
129
         end if
 
130
C
 
131
         IF( NDNMBR ) THEN
 
132
            IF( GTNDFL(BIGLST,NODLST,IPTR,8) ) THEN
 
133
               JNOD = JNOD - 1
 
134
               KNOD = JNOD
 
135
            ELSE
 
136
               INOD = INOD + 1
 
137
               KNOD = INOD
 
138
            END IF
 
139
         ELSE
 
140
            KNOD = INT(NODLST(17,IPTR))
 
141
            if( knod .le. 0 ) then
 
142
               WRITE(0,*) '*** MKFIXD: Got node again!'
 
143
               WRITE(0,*) knod,iptr,inod,jnod,nnod
 
144
               ierr = -1987
 
145
               return
 
146
            else if( knod .gt. nnod ) then
 
147
               WRITE(0,*)'*** MKFIXD: Got out-of-range KNOD: ',knod,nnod
 
148
               ierr = -1987
 
149
               return
 
150
            end if
 
151
         END IF
 
152
C
 
153
C - rescale and shift the domain back to what it was
 
154
C
 
155
         NODX(KNOD) = NODLST(1,IPTR)*XSCALE + XSHIFT
 
156
         NODY(KNOD) = NODLST(2,IPTR)*YSCALE + YSHIFT
 
157
         NODZ(KNOD) = NODLST(3,IPTR)*ZSCALE + ZSHIFT
 
158
C
 
159
         IF( SIZMTX .EQ. 9*NNOD ) THEN
 
160
            I = KNOD*9-9
 
161
            ORGMTX(I+1) = NODLST( 7,IPTR)/XSCALE/XSCALE
 
162
            ORGMTX(I+2) = NODLST( 8,IPTR)/XSCALE/YSCALE
 
163
            ORGMTX(I+3) = NODLST( 9,IPTR)/XSCALE/ZSCALE
 
164
            ORGMTX(I+4) = NODLST(10,IPTR)/YSCALE/XSCALE
 
165
            ORGMTX(I+5) = NODLST(11,IPTR)/YSCALE/YSCALE
 
166
            ORGMTX(I+6) = NODLST(12,IPTR)/YSCALE/ZSCALE
 
167
            ORGMTX(I+7) = NODLST(13,IPTR)/ZSCALE/XSCALE
 
168
            ORGMTX(I+8) = NODLST(14,IPTR)/ZSCALE/YSCALE
 
169
            ORGMTX(I+9) = NODLST(15,IPTR)/ZSCALE/ZSCALE
 
170
C            DO I = 1, 9
 
171
C               ORGMTX(I+KNOD*9-9) = NODLST(6+I,IPTR)
 
172
C            END DO
 
173
         END IF
 
174
C
 
175
         NODLST(17,IPTR) = -FLOAT(KNOD)
 
176
C
 
177
         IF( SIZLOC .EQ. NNOD ) NODLOC(KNOD) = INT(NODLST(16,IPTR))
 
178
c         print*,'Node ',nodx(Knod),nody(Knod),nodz(Knod),iptr
 
179
         if(nodlst(nxtnod,iptr).eq.iptr) then
 
180
            WRITE(0,*) '*** MKFIXD: Node points to itself!'
 
181
            WRITE(0,*) iptr,inod,jnod,knod,nnod
 
182
            ierr = -1876
 
183
            return
 
184
         end if
 
185
         IPTR = NODLST(NXTNOD,IPTR)
 
186
         if(iptr.le.topnod .and. iptr.gt.0) goto 10
 
187
C
 
188
         if(nodcnt.lt.nnod) then
 
189
            WRITE(0,*) '*** MKFIXD: Reached end of node list too early'
 
190
            WRITE(0,*) nodcnt,nnod,iptr,topnod
 
191
            ierr = -9696
 
192
            return
 
193
         end if
 
194
C
 
195
      END IF
 
196
C
 
197
      IF( IERR .NE. 0 ) THEN
 
198
         WRITE(0,*) '*** MKFIXD: GOT ERROR FROM GTNDFL'
 
199
         RETURN
 
200
      END IF
 
201
c
 
202
      iptr = sttnod
 
203
      inod = 0
 
204
  11  if( gtndfl(biglst,nodlst,iptr,1) ) then
 
205
         inod = inod + 1
 
206
         knod = int(nodlst(17,iptr))
 
207
         if(knod.ge.0) then
 
208
            WRITE(0,*) '*** MKFIXD: Found +ve entry!'
 
209
            WRITE(0,*) nodlst(17,iptr),iptr,inod,nnod
 
210
            ierr = -17856
 
211
            return
 
212
         end if
 
213
         if(inod.gt.nnod) then
 
214
            WRITE(0,*) '*** MKFIXD: Beyond end of nodes!'
 
215
            WRITE(0,*) nodlst(17,iptr),iptr,inod,nnod
 
216
            ierr = -16847
 
217
            return
 
218
         end if
 
219
         knod = -knod
 
220
         if(knod.gt.nnod) then
 
221
            WRITE(0,*) '*** MKFIXD: Node number too big!'
 
222
            WRITE(0,*) nodlst(17,iptr),iptr,inod,nnod
 
223
            ierr = -13968
 
224
            return
 
225
         end if
 
226
         nodlst(17,iptr) = float(knod)
 
227
         iptr = nodlst(nxtnod,iptr)
 
228
         if(iptr.le.topnod .and. iptr.gt.0) goto 11
 
229
      end if
 
230
C
 
231
      IF( NPRDND .GT. 1 ) THEN
 
232
C
 
233
         DO I = 1, NPRDND
 
234
            j = prdnds(i)
 
235
            PRDNDS(I) = INT( NODLST(17,PRDNDS(I)) )
 
236
            ewrite(2,*) 
 
237
     :           '  periodic renumber: ',j,prdnds(i)
 
238
         END DO
 
239
C
 
240
      END IF
 
241
c
 
242
      ewrite(2,*) 
 
243
     :     'Finished new node data'
 
244
      ewrite(2,*) 
 
245
     :     'Geom,split,int: ',cntgmy,cntspl,cntint
 
246
C
 
247
      CNTINT = 0
 
248
      CNTGMY = 0
 
249
      IPTR = STTBIG
 
250
      IELM = 0
 
251
      IPOS = 0
 
252
      SELM = 0
 
253
      SPOS = 0
 
254
      ENLBAS(IELM+1) = IPOS
 
255
      SNLBAS(SELM+1) = SPOS
 
256
      mxsfid = 0
 
257
      allocate( srfcnt(0:MAX_SFID_STATS) )
 
258
      allocate( srfavx(0:MAX_SFID_STATS) )
 
259
      allocate( srfavy(0:MAX_SFID_STATS) )
 
260
      allocate( srfavz(0:MAX_SFID_STATS) )
 
261
      srfavx = 0.0
 
262
      srfavy = 0.0
 
263
      srfavz = 0.0
 
264
      srfcnt = 0
 
265
c      ewrite(3,*) 'cleared srfavs'
 
266
C
 
267
  20  IF( BIGLST(4,IPTR) .GT. 0 ) THEN
 
268
C
 
269
         IELM = IELM + 1
 
270
C
 
271
         IF( IELM .GT. NELM ) THEN
 
272
            WRITE(0,*) '*** MKFIXD: ERROR COUNTING ELEMENTS'
 
273
            WRITE(0,*) NELM,IELM,IPTR
 
274
            IERR = -8989
 
275
            deallocate(srfcnt, srfavx, srfavy, srfavz)
 
276
            RETURN
 
277
         END IF
 
278
C
 
279
         CALL ELMINF( BIGLST, NODLST, IPTR, CONELM, EDGLST,
 
280
     :                INEW, IFLAG, IREG, NXT, LST )
 
281
C
 
282
         ELMREG(IELM) = IREG
 
283
C
 
284
         CALL ELNODS( BIGLST, NODLST, IPTR, NODS, NXT, LST, .TRUE. )
 
285
C
 
286
         i1 = 1
 
287
         i2 = 2
 
288
         IF( IS3DMS ) THEN
 
289
            X(1) = NODLST(1,NODS(1))
 
290
            X(2) = NODLST(1,NODS(2))
 
291
            X(3) = NODLST(1,NODS(3))
 
292
            X(4) = NODLST(1,NODS(4))
 
293
            Y(1) = NODLST(2,NODS(1))
 
294
            Y(2) = NODLST(2,NODS(2))
 
295
            Y(3) = NODLST(2,NODS(3))
 
296
            Y(4) = NODLST(2,NODS(4))
 
297
            Z(1) = NODLST(3,NODS(1))
 
298
            Z(2) = NODLST(3,NODS(2))
 
299
            Z(3) = NODLST(3,NODS(3))
 
300
            Z(4) = NODLST(3,NODS(4))
 
301
            IF( TETVOL( X, Y, Z ) .LT. 0.0 ) THEN
 
302
c               print*,'*** MKFIXD: Got inside-out tet...'
 
303
c               STOP
 
304
c               print*,ielm,iptr
 
305
c               print*,nods
 
306
c               print*,x
 
307
c               print*,y
 
308
c               print*,z
 
309
               i1 = 2
 
310
               i2 = 1
 
311
            END IF
 
312
         END IF
 
313
C
 
314
         ENLIST(IPOS+1) = NODLST(17,NODS(i1))
 
315
         ENLIST(IPOS+2) = NODLST(17,NODS(i2))
 
316
         ENLIST(IPOS+3) = NODLST(17,NODS(3))
 
317
C
 
318
         IPOS = IPOS + 3
 
319
         IF( IS3DMS ) THEN
 
320
            ENLIST(IPOS+1) = NODLST(17,NODS(4))
 
321
            IPOS = IPOS + 1
 
322
         END IF
 
323
         if( ielm .lt. 11 .and. debug ) then
 
324
            print*,'old elm ',iptr,(nods(i),i=1,4)
 
325
            print*,'elm ',ielm,(enlist(ipos+i-4),i=1,4)
 
326
         end if
 
327
C
 
328
         IF( IPOS .GT. SZENLS+1 ) THEN
 
329
            WRITE(0,*) '*** MKFIXD: BEYOND END OF ENLIST'
 
330
            WRITE(0,*) SZENLS,IPOS,IELM,NELM,IPTR
 
331
            IERR = -9898
 
332
            deallocate(srfcnt, srfavx, srfavy, srfavz)
 
333
            RETURN
 
334
         END IF
 
335
C
 
336
         ENLBAS(IELM+1) = IPOS
 
337
C
 
338
         IF( IS3DMS ) THEN
 
339
C
 
340
c            IF( CONELM(1) .LE. 0 .OR. CONELM(2) .LE. 0 .OR.
 
341
c     :          CONELM(3) .LE. 0 .OR. CONELM(4) .EQ. 0 ) THEN
 
342
C
 
343
               IF( CONELM(1) .LE. 0 ) THEN
 
344
                  SELM = SELM + 1
 
345
                  SNLIST(SPOS+i1) = NODLST(17,NODS(2))
 
346
                  SNLIST(SPOS+i2) = NODLST(17,NODS(1))
 
347
                  SNLIST(SPOS+3)  = NODLST(17,NODS(3))
 
348
c                  if( selm .lt. 11 ) then
 
349
c                     print*,'srf ',selm,(snlist(spos+i),i=1,3)
 
350
c                  end if
 
351
                  SPOS = SPOS + 3
 
352
                  SNLBAS(SELM+1) = SPOS
 
353
                  SURFID(SELM) = -CONELM(1)-1
 
354
                  IF( SIZOWN .EQ. NSELM ) SUROWN(SELM) = IELM
 
355
                  if( ielm.gt.nelm .or. ielm.le.0 )
 
356
     :                WRITE(0,*) 'MKFIXD: IELM too big!',ielm,nelm,selm
 
357
                  mxsfid = max(mxsfid,surfid(selm))
 
358
                  if( surfid(selm) .lt. 0 ) then
 
359
                     WRITE(0,*) 'BAD SURFID RANGE: ',surfid(selm),selm
 
360
                  else if( surfid(selm) .le. MAX_SFID_STATS ) then
 
361
c                     if(surfid(selm).gt.0) elmreg(ielm) = 3
 
362
                     srfcnt(surfid(selm))=srfcnt(surfid(selm))+1
 
363
                     srfavx(surfid(selm))=srfavx(surfid(selm))
 
364
     :                                   +x(1)+x(2)+x(3)
 
365
                     srfavy(surfid(selm))=srfavy(surfid(selm))
 
366
     :                                   +y(1)+y(2)+y(3)
 
367
                     srfavz(surfid(selm))=srfavz(surfid(selm))
 
368
     :                                   +z(1)+z(2)+z(3)
 
369
                  end if
 
370
               END IF
 
371
C
 
372
               IF( CONELM(2) .LE. 0 ) THEN
 
373
                  SELM = SELM + 1
 
374
                  SNLIST(SPOS+i1) = NODLST(17,NODS(1))
 
375
                  SNLIST(SPOS+i2) = NODLST(17,NODS(2))
 
376
                  SNLIST(SPOS+3)  = NODLST(17,NODS(4))
 
377
c                  if( selm .lt. 11 ) then
 
378
c                     print*,'srf ',selm,(snlist(spos+i),i=1,3)
 
379
c                  end if
 
380
                  SPOS = SPOS + 3
 
381
                  SNLBAS(SELM+1) = SPOS
 
382
                  SURFID(SELM) = -CONELM(2)-1
 
383
                  IF( SIZOWN .EQ. NSELM ) SUROWN(SELM) = IELM
 
384
                  if( ielm.gt.nelm .or. ielm.le.0 )
 
385
     :                WRITE(0,*) 'MKFIXD: IELM too big!',ielm,nelm,selm
 
386
                  mxsfid = max(mxsfid,surfid(selm))
 
387
                  if( surfid(selm) .lt. 0 ) then
 
388
                     WRITE(0,*) 'BAD SURFID RANGE: ',surfid(selm),selm
 
389
                  else if( surfid(selm) .le. MAX_SFID_STATS ) then
 
390
c                     if(surfid(selm).gt.0) elmreg(ielm) = 3
 
391
                     srfcnt(surfid(selm))=srfcnt(surfid(selm))+1
 
392
                     srfavx(surfid(selm))=srfavx(surfid(selm))
 
393
     :                                   +x(1)+x(2)+x(3)
 
394
                     srfavy(surfid(selm))=srfavy(surfid(selm))
 
395
     :                                   +y(1)+y(2)+y(3)
 
396
                     srfavz(surfid(selm))=srfavz(surfid(selm))
 
397
     :                                   +z(1)+z(2)+z(3)
 
398
                  end if
 
399
               END IF
 
400
C
 
401
               IF( CONELM(3) .LE. 0 ) THEN
 
402
                  SELM = SELM + 1
 
403
                  SNLIST(SPOS+i1) = NODLST(17,NODS(3))
 
404
                  SNLIST(SPOS+i2) = NODLST(17,NODS(1))
 
405
                  SNLIST(SPOS+3)  = NODLST(17,NODS(4))
 
406
c                  if( selm .lt. 11 ) then
 
407
c                     print*,'srf ',selm,(snlist(spos+i),i=1,3)
 
408
c                  end if
 
409
                  SPOS = SPOS + 3
 
410
                  SNLBAS(SELM+1) = SPOS
 
411
                  SURFID(SELM) = -CONELM(3)-1
 
412
                  IF( SIZOWN .EQ. NSELM ) SUROWN(SELM) = IELM
 
413
                  if( ielm.gt.nelm .or. ielm.le.0 )
 
414
     :                WRITE(0,*) 'MKFIXD: IELM too big!',ielm,nelm,selm
 
415
                  mxsfid = max(mxsfid,surfid(selm))
 
416
                  if( surfid(selm) .lt. 0 ) then
 
417
                     WRITE(0,*) 'BAD SURFID RANGE: ',surfid(selm),selm
 
418
                  else if( surfid(selm) .le. MAX_SFID_STATS ) then
 
419
c                     if(surfid(selm).gt.0) elmreg(ielm) = 3
 
420
                     srfcnt(surfid(selm))=srfcnt(surfid(selm))+1
 
421
                     srfavx(surfid(selm))=srfavx(surfid(selm))
 
422
     :                                   +x(1)+x(2)+x(3)
 
423
                     srfavy(surfid(selm))=srfavy(surfid(selm))
 
424
     :                                   +y(1)+y(2)+y(3)
 
425
                     srfavz(surfid(selm))=srfavz(surfid(selm))
 
426
     :                                   +z(1)+z(2)+z(3)
 
427
                  end if
 
428
               END IF
 
429
C
 
430
               IF( CONELM(4) .LE. 0 ) THEN
 
431
                  SELM = SELM + 1
 
432
c                  if(facs(4).eq.4) then
 
433
                     WRITE(0,*) '+++ MKFIXD: Eh??? Surface on face 4!'
 
434
                     WRITE(0,*) ielm,selm
 
435
                     WRITE(0,*) conelm
 
436
                     WRITE(0,*) nods
 
437
c                  end if
 
438
                  SNLIST(SPOS+i1) = NODLST(17,NODS(2))
 
439
                  SNLIST(SPOS+i2) = NODLST(17,NODS(3))
 
440
                  SNLIST(SPOS+3)  = NODLST(17,NODS(4))
 
441
c                  if( selm .lt. 11 ) then
 
442
c                     print*,'srf ',selm,(snlist(spos+i),i=1,3)
 
443
c                  end if
 
444
                  SPOS = SPOS + 3
 
445
                  SNLBAS(SELM+1) = SPOS
 
446
                  SURFID(SELM) = -CONELM(4)-1
 
447
                  IF( SIZOWN .EQ. NSELM ) SUROWN(SELM) = IELM
 
448
                  if( ielm.gt.nelm .or. ielm.le.0 )
 
449
     :                WRITE(0,*) 'MKFIXD: IELM too big!',ielm,nelm,selm
 
450
                  mxsfid = max(mxsfid,surfid(selm))
 
451
                  if( surfid(selm) .lt. 0 ) then
 
452
                     WRITE(0,*) 'BAD SURFID RANGE: ',surfid(selm),selm
 
453
                  else if( surfid(selm) .le. MAX_SFID_STATS ) then
 
454
c                     if(surfid(selm).gt.0) elmreg(ielm) = 3
 
455
                     srfcnt(surfid(selm))=srfcnt(surfid(selm))+1
 
456
                     srfavx(surfid(selm))=srfavx(surfid(selm))
 
457
     :                                   +x(1)+x(2)+x(3)
 
458
                     srfavy(surfid(selm))=srfavy(surfid(selm))
 
459
     :                                   +y(1)+y(2)+y(3)
 
460
                     srfavz(surfid(selm))=srfavz(surfid(selm))
 
461
     :                                   +z(1)+z(2)+z(3)
 
462
                  end if
 
463
               END IF
 
464
C
 
465
               IF( SELM .GT. NSELM ) THEN
 
466
                 WRITE(0,*)'*** MKFIXD: ERROR COUNTING SURFACE ELEMENTS'
 
467
                 WRITE(0,*) NSELM,SELM
 
468
                 STOP
 
469
               END IF
 
470
C
 
471
               IF( SPOS .GT. SZSNLS+1 ) THEN
 
472
                  WRITE(0,*) '*** MKFIXD: BEYOND END OF SNLIST'
 
473
                  WRITE(0,*) SZSNLS,SPOS
 
474
                  STOP
 
475
               END IF
 
476
C
 
477
c            END IF
 
478
C
 
479
C         ELSE
 
480
C
 
481
C
 
482
         END IF
 
483
C
 
484
         IPTR = NXT
 
485
         IF( IPTR .LE. TOPBIG .AND. IPTR .GT. 0 ) GOTO 20
 
486
C
 
487
         IF( IELM .NE. NELM ) THEN
 
488
           WRITE(0,*)'*** MKFIXD: Reached end of element list too early'
 
489
           WRITE(0,*) ielm,nelm,iptr,topbig
 
490
           ierr = -9595
 
491
           deallocate(srfcnt, srfavx, srfavy, srfavz)
 
492
           return
 
493
         end if
 
494
C
 
495
      ELSE IF( BIGLST(4,IPTR) .LT. 0 ) THEN
 
496
C
 
497
C         CALL EDGINF( BIGLST, NODLST, IPTR, NODS,
 
498
C     :                INEW, IFLAG, NXT, LST )
 
499
C
 
500
         if( getflg(biglst,nodlst,iptr,2) ) then
 
501
            cntgmy = cntgmy + 1
 
502
         else if( getflg(biglst,nodlst,iptr,3) ) then
 
503
            cntint = cntint + 1
 
504
         end if
 
505
C
 
506
         IPTR = BIGLST( NXTBIG, IPTR )
 
507
         IF( IPTR .LE. TOPBIG .AND. IPTR .GT. 0 ) GOTO 20
 
508
C
 
509
         IF( IELM .NE. NELM ) THEN
 
510
           WRITE(0,*)'*** MKFIXD: Reached end of element list too early'
 
511
           WRITE(0,*) ielm,nelm,iptr,topbig
 
512
           ierr = -9595
 
513
           deallocate(srfcnt, srfavx, srfavy, srfavz)
 
514
           return
 
515
         end if
 
516
C
 
517
      END IF
 
518
c
 
519
      ewrite(2,*) 
 
520
     :     'Finished new element data'
 
521
      ewrite(2,*) 
 
522
     :     'Geom ed, int ed: ',cntgmy,cntint
 
523
      nusfid = count(srfcnt>0)
 
524
      if(nusfid.lt.200) then
 
525
         ewrite(2,*)
 
526
     :       'Surface averages  ( max id :',mxsfid,' )'
 
527
         do i = 0, MAX_SFID_STATS
 
528
            if( srfcnt(i) .gt. 0 ) then
 
529
               ewrite(2,*) 
 
530
     :              i,srfcnt(i),srfavx(i)/srfcnt(i)/3,
 
531
     :              srfavy(i)/srfcnt(i)/3,srfavz(i)/srfcnt(i)/3
 
532
            end if
 
533
         end do
 
534
      else
 
535
         ewrite(2,*) 
 
536
     :        'Check in MKFIXD:  NSELM, MXSFID, NUSFID =',
 
537
     :                                       nselm,mxsfid,nusfid
 
538
      end if
 
539
C
 
540
      IF( NPROC .GT. 1 ) THEN
 
541
C
 
542
         ewrite(2,*) 
 
543
     :        'Working out new gather array...'
 
544
         DO I = 1, NGATH
 
545
            INOD = GATHER(I)
 
546
            INOD = INT(NODLST(17,INOD))
 
547
            if(inod.gt.nnod) then
 
548
               WRITE(0,*) '  old:',gather(i),'  new:',inod,
 
549
     :                                   ' ***TOO BIG!!!',nnod
 
550
            else if(inod.lt.1) then
 
551
               WRITE(0,*) '  old:',gather(i),'  new:',inod,
 
552
     :                                   ' ***TOO SMALL!!!'
 
553
            else if( debug ) then
 
554
              WRITE(0,*) '  old:',gather(i),'  new:',
 
555
     :                                     inod,nodlst(6,inod)
 
556
            end if
 
557
            GATHER(I) = INOD
 
558
         END DO
 
559
C
 
560
         ewrite(2,*) 
 
561
     :        'Working out new scatter array...'
 
562
         DO I = 1, NHALO
 
563
            INOD = SCATER(I)
 
564
            INOD = INT(NODLST(17,INOD))
 
565
            if(inod.gt.nnod) then
 
566
               WRITE(0,*) '  old:',scater(i),'  new:',inod,
 
567
     :                                   ' ***TOO BIG!!!',nnod
 
568
            else if(inod.lt.1) then
 
569
               WRITE(0,*) '  old:',scater(i),'  new:',inod,
 
570
     :                                   ' ***TOO SMALL!!!'
 
571
            else if( debug ) then
 
572
               ewrite(2,*) 
 
573
     :              '  old:',scater(i),'  new:',
 
574
     :              inod,nodlst(6,inod)
 
575
            end if
 
576
            SCATER(I) = INOD
 
577
         END DO
 
578
C
 
579
      END IF
 
580
      ewrite(2,*) 
 
581
     :     'Checking local node ordering...'
 
582
C
 
583
      do i = 2, 8
 
584
         m(i) = 0.0
 
585
      end do
 
586
      m(1) = 1.0
 
587
      m(5) = 1.0
 
588
      m(9) = 1.0
 
589
c
 
590
      MAXASP = 0.0
 
591
      ASPAVE = 0.0
 
592
      VOLSUM = 0.0
 
593
      MAXVOL = 0.0
 
594
      MINVOL = 1E+20
 
595
      MINRAD = 1E+20
 
596
C - to avoid compiler warnings...
 
597
      MINELR = 0
 
598
      MAXELA = 0
 
599
      MAXELV = 0
 
600
      MINELV = 0
 
601
C
 
602
      DO I = 1, NELM
 
603
         IPOS = 0
 
604
c         print*,'element ',i,enlbas(i)+1
 
605
         DO J = ENLBAS(I)+1, ENLBAS(I+1)
 
606
            IPOS = IPOS + 1
 
607
            if(enlist(j).gt.nnod.or.enlist(j).le.0)
 
608
     :         WRITE(0,*) 'MKFIXD: bad node: ',enlist(j),i,j,ipos
 
609
            X(IPOS) = NODX(ENLIST(J))
 
610
            Y(IPOS) = NODY(ENLIST(J))
 
611
            Z(IPOS) = NODZ(ENLIST(J))
 
612
c            print*,'   node ',ipos,enlist(j),x(ipos),y(ipos),z(ipos)
 
613
         END DO
 
614
         CALL MTETIN( X, Y, Z, M, AVOL, AREAS, L, ARAD, ASP )
 
615
c         print*,'   vol,rad,asp ',avol,arad,asp
 
616
         ASPAVE = ASPAVE + ASP
 
617
         VOLSUM = VOLSUM + AVOL
 
618
         IF( AVOL .LE. 0.0 ) THEN
 
619
            WRITE(0,*) '+++ FOUND NEGATIVE VOLUME!',avol
 
620
            WRITE(0,*) I,(ENLIST(J),J=ENLBAS(I)+1,ENLBAS(I+1))
 
621
            J = ENLBAS(I) + 1
 
622
            IPOS = ENLIST(J)
 
623
            ENLIST(J)   = ENLIST(J+1)
 
624
            ENLIST(J+1) = IPOS
 
625
         END IF
 
626
         IF( ASP .GT. MAXASP ) THEN
 
627
            MAXASP = ASP
 
628
            MAXELA = I
 
629
            MAXAVL = AVOL
 
630
            MAXARD = ARAD
 
631
         END IF
 
632
         IF( AVOL .LT. MINVOL ) THEN
 
633
            MINVOL = AVOL
 
634
            MINELV = I
 
635
            MINVRD = ARAD
 
636
            MINVAS = ASP
 
637
         END IF
 
638
         IF( AVOL .GT. MAXVOL ) THEN
 
639
            MAXVOL = AVOL
 
640
            MAXELV = I
 
641
            MAXVRD = ARAD
 
642
            MAXVAS = ASP
 
643
         END IF
 
644
         IF( ARAD .LT. MINRAD ) THEN
 
645
            MINRAD = ARAD
 
646
            MINELR = I
 
647
            MINRVL = AVOL
 
648
            MINRAS = ASP
 
649
         END IF
 
650
      END DO
 
651
      ewrite(2,*) 
 
652
     :     'VOLSUM,ASPAVE: ',VOLSUM,ASPAVE/NELM
 
653
      ewrite(2,*) 
 
654
     :     'Maximum asp,vol,rad: ',maxasp,maxavl,maxard
 
655
      if( debug ) then
 
656
        ipos = enlbas(maxela)
 
657
        inod = enlist(ipos+1)
 
658
        print*,'  node 1 : ',nodx(inod),nody(inod),nodz(inod)
 
659
        inod = enlist(ipos+2)
 
660
        print*,'  node 2 : ',nodx(inod),nody(inod),nodz(inod)
 
661
        inod = enlist(ipos+3)
 
662
        print*,'  node 3 : ',nodx(inod),nody(inod),nodz(inod)
 
663
        inod = enlist(ipos+4)
 
664
        print*,'  node 4 : ',nodx(inod),nody(inod),nodz(inod)
 
665
      end if
 
666
      ewrite(2,*) 
 
667
     :     'Maximum vol,rad,asp: ',maxvol,maxvrd,maxvas
 
668
      if( debug ) then
 
669
        ipos = enlbas(maxelv)
 
670
        inod = enlist(ipos+1)
 
671
        print*,'  node 1 : ',nodx(inod),nody(inod),nodz(inod)
 
672
        inod = enlist(ipos+2)
 
673
        print*,'  node 2 : ',nodx(inod),nody(inod),nodz(inod)
 
674
        inod = enlist(ipos+3)
 
675
        print*,'  node 3 : ',nodx(inod),nody(inod),nodz(inod)
 
676
        inod = enlist(ipos+4)
 
677
        print*,'  node 4 : ',nodx(inod),nody(inod),nodz(inod)
 
678
      end if
 
679
      ewrite(2,*) 
 
680
     :     'Minimum vol,rad,asp: ',minvol,minvrd,minvas
 
681
      if( debug ) then
 
682
        ipos = enlbas(minelv)
 
683
        inod = enlist(ipos+1)
 
684
        print*,'  node 1 : ',nodx(inod),nody(inod),nodz(inod)
 
685
        inod = enlist(ipos+2)
 
686
        print*,'  node 2 : ',nodx(inod),nody(inod),nodz(inod)
 
687
        inod = enlist(ipos+3)
 
688
        print*,'  node 3 : ',nodx(inod),nody(inod),nodz(inod)
 
689
        inod = enlist(ipos+4)
 
690
        print*,'  node 4 : ',nodx(inod),nody(inod),nodz(inod)
 
691
      end if
 
692
      ewrite(2,*) 
 
693
     :     'Minimum rad,vol,asp: ',minrad,minrvl,minras
 
694
      if( debug ) then
 
695
        ipos = enlbas(minelr)
 
696
        inod = enlist(ipos+1)
 
697
        print*,'  node 1 : ',nodx(inod),nody(inod),nodz(inod)
 
698
        inod = enlist(ipos+2)
 
699
        print*,'  node 2 : ',nodx(inod),nody(inod),nodz(inod)
 
700
        inod = enlist(ipos+3)
 
701
        print*,'  node 3 : ',nodx(inod),nody(inod),nodz(inod)
 
702
        inod = enlist(ipos+4)
 
703
        print*,'  node 4 : ',nodx(inod),nody(inod),nodz(inod)
 
704
      end if
 
705
c
 
706
c                do i = 1, nelm
 
707
c                   iptr = 0
 
708
c                   do j = enlbas(i)+1, enlbas(i+1)
 
709
c                      iptr = iptr + 1
 
710
c                      xxx(iptr) = nodx(enlist(j))
 
711
c                      yyy(iptr) = nody(enlist(j))
 
712
c                      zzz(iptr) = nodz(enlist(j))
 
713
c                   end do
 
714
c                   aa = tetvol( xxx, yyy, zzz )
 
715
c                   if( aa .lt. 0.0 ) then
 
716
c                      j = enlbas(i)+1
 
717
c                      iptr = enlist(j)
 
718
c                      enlist(j) = enlist(j+1)
 
719
c                      enlist(j+1) = iptr
 
720
c                   end if
 
721
c                end do
 
722
c
 
723
c      print*,'ENLBAS: ',(enlbas(i),i=1,nelm+1)
 
724
c      print*,'ENLIST: ',(enlist(i),i=1,szenls)
 
725
c
 
726
      ewrite(2,*) 
 
727
     :     'Finished checking local node ordering'
 
728
C
 
729
  6   format( a, i7, a )
 
730
  7   format( a, 3(1pe13.5) )
 
731
  8   format( 2i7, 3(1pe13.5) )
 
732
  9   format( a, 2i7 )
 
733
C
 
734
      deallocate(srfcnt, srfavx, srfavy, srfavz)
 
735
      RETURN
 
736
      END
 
737
C