← Back to branch summary

~airpollution/fluidity/fluidity_airpollution

« back to all changes in this revision

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

  • Committer: ziyouzhj
  • Date: 2013-12-09 16:51:29 UTC
  • Revision ID: ziyouzhj@gmail.com-20131209165129-ucoetc3u0atyy05c
airpolution

Show diffs side-by-side

added added

removed removed

Lines of Context:
libadaptivity/adapt3d/src/elmfnc.F
 
1
#include "fdebug_adapt.h"
 
2
 
 
3
!#define LIPNIKOV_FUNCTIONAL
 
4
#ifdef LIPNIKOV_FUNCTIONAL
 
5
      function elmfnc(biglst, nodlst, ielm, nod1, nod2, nod3, nod4,
 
6
     &         undlen) result(functional)
 
7
 
 
8
      use addebug
 
9
      implicit none
 
10
      include 'blknew.i'
 
11
      include 'blkbig.i'
 
12
 
 
13
      integer :: ielm, nod1, nod2, nod3, nod4
 
14
      real :: undlen
 
15
      integer :: cache 
 
16
      integer :: ielm2 
 
17
      real :: functional
 
18
      real, parameter :: m=2.0*huge(0)-3.0, c= -1.0*huge(0)+2.0 ! mapping between reals and ints for the cache
 
19
      real, dimension(9) :: metric
 
20
      real, dimension(3) :: tmppos
 
21
      real, dimension(3, 3) :: positions
 
22
      real :: volume ! volume with respect to the metric
 
23
      real :: metric_determinant ! determinant of the metric
 
24
      real :: tetvol ! volume of the tet in physical space
 
25
      real :: edgelens ! 2-norm of vector of edge lengths in metric space
 
26
      real :: F ! from Agouzal et al, 1999
 
27
      real :: x
 
28
      real, parameter :: factor = 6.0**(2.5) * sqrt(2.0)
 
29
      logical :: negate
 
30
      integer, dimension(4) :: nods
 
31
      integer :: k
 
32
      real, parameter :: bst = 12.0 / sqrt(2.0) ! inverse of ideal tet volume
 
33
      integer :: intval ! undlen cast to an integer
 
34
 
 
35
 
 
36
C - Sort the nodes into ascending order
 
37
C
 
38
      NODS(1) = NOD1
 
39
      NODS(2) = NOD2
 
40
      NODS(3) = NOD3
 
41
      NODS(4) = NOD4
 
42
      NEGATE = .false.
 
43
C
 
44
      K = NODS(1)
 
45
      IF( K .GT. NODS(2) ) THEN
 
46
        NODS(1) = NODS(2)
 
47
        NODS(2) = K
 
48
        NEGATE = .NOT. NEGATE
 
49
      ELSE
 
50
        K = NODS(2)
 
51
      END IF
 
52
      IF( K .GT. NODS(3) ) THEN
 
53
        NODS(2) = NODS(3)
 
54
        NODS(3) = K
 
55
        NEGATE = .NOT. NEGATE
 
56
      ELSE
 
57
        K = NODS(3)
 
58
      END IF
 
59
      IF( K .GT. NODS(4) ) THEN
 
60
        NODS(3) = NODS(4)
 
61
        NODS(4) = K
 
62
        NEGATE = .NOT. NEGATE
 
63
      END IF
 
64
      K = NODS(1)
 
65
      IF( K .GT. NODS(2) ) THEN
 
66
        NODS(1) = NODS(2)
 
67
        NODS(2) = K
 
68
        NEGATE = .NOT. NEGATE
 
69
      ELSE
 
70
        K = NODS(2)
 
71
      END IF
 
72
      IF( K .GT. NODS(3) ) THEN
 
73
        NODS(2) = NODS(3)
 
74
        NODS(3) = K
 
75
        NEGATE = .NOT. NEGATE
 
76
      END IF
 
77
      K = NODS(1)
 
78
      IF( K .GT. NODS(2) ) THEN
 
79
        NODS(1) = NODS(2)
 
80
        NODS(2) = K
 
81
        NEGATE = .NOT. NEGATE
 
82
      END IF
 
83
 
 
84
      if (ielm > 0) then
 
85
        cache = biglst(7, ielm)
 
86
        ielm2 = biglst(nxtbig, ielm)
 
87
      else
 
88
        cache = 0
 
89
        ielm2 = -1
 
90
      end if
 
91
 
 
92
      ! if we have no cached information, then we need to compute it, na ja?
 
93
      if (cache == 0) then
 
94
        ! first things first, let us fetch the metric and positions
 
95
        tmppos = nodlst(1:3, nods(1))
 
96
        metric = nodlst(7:15, nods(1))
 
97
 
 
98
        positions(:, 1) = nodlst(1:3, nods(2)) - tmppos
 
99
        metric = metric + nodlst(7:15, nods(2))
 
100
 
 
101
        positions(:, 2) = nodlst(1:3, nods(3)) - tmppos
 
102
        metric = metric + nodlst(7:15, nods(3))
 
103
 
 
104
        positions(:, 3) = nodlst(1:3, nods(4)) - tmppos
 
105
        metric = metric + nodlst(7:15, nods(4))
 
106
 
 
107
        metric = metric / 4.0
 
108
 
 
109
        metric_determinant = 
 
110
     &  metric(1) * (metric(5) * metric(9) - metric(6) * metric(8))
 
111
     &- metric(2) * (metric(4) * metric(9) - metric(6) * metric(7))
 
112
     &+ metric(3) * (metric(4) * metric(8) - metric(5) * metric(7))
 
113
 
 
114
        tetvol = 
 
115
     &  positions(1,1) * (positions(2,2) * positions(3,3) - positions(2,3) * positions(3,2))
 
116
     &- positions(1,2) * (positions(2,1) * positions(3,3) - positions(2,3) * positions(3,1))
 
117
     &+ positions(1,3) * (positions(2,1) * positions(3,2) - positions(2,2) * positions(3,1))
 
118
        tetvol = tetvol / 6.0
 
119
 
 
120
        undlen = sqrt(metric_determinant) * tetvol
 
121
!        if(negate) undlen = - undlen
 
122
        volume = abs(undlen)
 
123
        undlen = undlen * bst
 
124
 
 
125
        !write(0,*) "computed undlen: ", undlen
 
126
        tmppos = positions(:, 1)
 
127
        edgelens =
 
128
     &  tmppos(1) * (metric(1) * tmppos(1) + metric(2) * tmppos(2) + metric(3) * tmppos(3))
 
129
     &+ tmppos(2) * (metric(4) * tmppos(1) + metric(5) * tmppos(2) + metric(6) * tmppos(3))
 
130
     &+ tmppos(3) * (metric(7) * tmppos(1) + metric(8) * tmppos(2) + metric(9) * tmppos(3))
 
131
 
 
132
        tmppos = positions(:, 2)
 
133
        edgelens = edgelens +
 
134
     &  tmppos(1) * (metric(1) * tmppos(1) + metric(2) * tmppos(2) + metric(3) * tmppos(3))
 
135
     &+ tmppos(2) * (metric(4) * tmppos(1) + metric(5) * tmppos(2) + metric(6) * tmppos(3))
 
136
     &+ tmppos(3) * (metric(7) * tmppos(1) + metric(8) * tmppos(2) + metric(9) * tmppos(3))
 
137
 
 
138
        tmppos = positions(:, 3)
 
139
        edgelens = edgelens +
 
140
     &  tmppos(1) * (metric(1) * tmppos(1) + metric(2) * tmppos(2) + metric(3) * tmppos(3))
 
141
     &+ tmppos(2) * (metric(4) * tmppos(1) + metric(5) * tmppos(2) + metric(6) * tmppos(3))
 
142
     &+ tmppos(3) * (metric(7) * tmppos(1) + metric(8) * tmppos(2) + metric(9) * tmppos(3))
 
143
 
 
144
        tmppos = positions(:, 1) - positions(:, 2)
 
145
        edgelens = edgelens +
 
146
     &  tmppos(1) * (metric(1) * tmppos(1) + metric(2) * tmppos(2) + metric(3) * tmppos(3))
 
147
     &+ tmppos(2) * (metric(4) * tmppos(1) + metric(5) * tmppos(2) + metric(6) * tmppos(3))
 
148
     &+ tmppos(3) * (metric(7) * tmppos(1) + metric(8) * tmppos(2) + metric(9) * tmppos(3))
 
149
 
 
150
        tmppos = positions(:, 1) - positions(:, 3)
 
151
        edgelens = edgelens +
 
152
     &  tmppos(1) * (metric(1) * tmppos(1) + metric(2) * tmppos(2) + metric(3) * tmppos(3))
 
153
     &+ tmppos(2) * (metric(4) * tmppos(1) + metric(5) * tmppos(2) + metric(6) * tmppos(3))
 
154
     &+ tmppos(3) * (metric(7) * tmppos(1) + metric(8) * tmppos(2) + metric(9) * tmppos(3))
 
155
 
 
156
        tmppos = positions(:, 2) - positions(:, 3)
 
157
        edgelens = edgelens +
 
158
     &  tmppos(1) * (metric(1) * tmppos(1) + metric(2) * tmppos(2) + metric(3) * tmppos(3))
 
159
     &+ tmppos(2) * (metric(4) * tmppos(1) + metric(5) * tmppos(2) + metric(6) * tmppos(3))
 
160
     &+ tmppos(3) * (metric(7) * tmppos(1) + metric(8) * tmppos(2) + metric(9) * tmppos(3))
 
161
 
 
162
        edgelens = sqrt(edgelens)
 
163
 
 
164
        x = edgelens / sqrt(6.0); x = min(x, 1.0/x)
 
165
        F = (x * (2-x))**3
 
166
        functional = factor * (volume / edgelens**3) * F
 
167
 
 
168
        undlen = min(max(undlen, -2.01E+3), 2.01E+3)
 
169
        intval = int(undlen * 1000000)
 
170
        undlen = float(intval) / 1000000
 
171
        if (ielm > 0) then
 
172
          biglst(7, ielm)  = intval
 
173
          biglst(7, ielm2) = int(m * functional + c)
 
174
        end if
 
175
      else
 
176
        undlen = float(cache) / 1000000
 
177
        functional = (biglst(7, ielm2) - c) / m
 
178
      end if
 
179
 
 
180
      if (negate) undlen = -undlen
 
181
      !write(0,*) "negate: ", negate, "; returning undlen: ", undlen
 
182
      assert(functional >= 0.0)
 
183
      assert(functional <= 1.0)
 
184
 
 
185
      ! And swap for adaptivity's convention
 
186
      functional = 5.0/max(functional, 5.0e-30) - 5.0
 
187
 
 
188
      end function elmfnc
 
189
#else
 
190
C Copyright (C) 2006 Imperial College London and others.
 
191
C 
 
192
C Please see the AUTHORS file in the main source directory for a full list
 
193
C of copyright holders.
 
194
C 
 
195
C Adrian Umpleby
 
196
C Applied Modelling and Computation Group
 
197
C Department of Earth Science and Engineering
 
198
C Imperial College London
 
199
C 
 
200
C adrian@Imperial.ac.uk
 
201
C 
 
202
C This library is free software; you can redistribute it and/or
 
203
C modify it under the terms of the GNU Lesser General Public
 
204
C License as published by the Free Software Foundation; either
 
205
C version 2.1 of the License.
 
206
C 
 
207
C This library is distributed in the hope that it will be useful,
 
208
C but WITHOUT ANY WARRANTY; without even the implied warranty of
 
209
C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 
210
C Lesser General Public License for more details.
 
211
C 
 
212
C You should have received a copy of the GNU Lesser General Public
 
213
C License along with this library; if not, write to the Free Software
 
214
C Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307
 
215
C USA
 
216
#include "ewrite.h"
 
217
      REAL FUNCTION ELMFNC( BIGLST, NODLST, IELM,
 
218
     :                      NOD1, NOD2, NOD3, NOD4, UNDLEN )
 
219
C-----------------------------------------------------------------------
 
220
C
 
221
C - This subroutine calculates the element functional.
 
222
C - The functional is calculated from the four node pointers (NOD1-NOD4)
 
223
C - which have to be first arranged in ascending order (this makes sure
 
224
C - that the element will have the same functional every time it is
 
225
C - calculated, by removing round-off errors because the order of nodes
 
226
C - is the same each time).
 
227
C - If IELM>0 and element IELM has not already had its functional worked
 
228
C - out then xeveral critical values (sum of edge functionals etc) are
 
229
C - worked out and stored for use next time the functional is required
 
230
C - for this element. This reduces the amount of calculation required.
 
231
C
 
232
C-----------------------------------------------------------------------
 
233
      use write_log
 
234
      IMPLICIT NONE
 
235
C
 
236
      REAL, PARAMETER :: BST = 4.898979486  ! sqrt(24)
 
237
C
 
238
      INCLUDE 'blknew.i'
 
239
C
 
240
      INCLUDE 'blkbig.i'
 
241
C
 
242
      INCLUDE 'blkerr.i'
 
243
C
 
244
      REAL ANGLEFACT
 
245
      COMMON / BLKTMP /  ANGLEFACT
 
246
C
 
247
      INTEGER IELM, NOD1, NOD2, NOD3, NOD4
 
248
C
 
249
      REAL UNDLEN
 
250
C
 
251
      integer i, j, k, nods(4), intval, nxt1, s
 
252
C
 
253
      REAL X(4), Y(4), Z(4), V1(3), V2(3), V3(3), D1, D2, D3, D, Q,
 
254
     :     M(9), AREAS(4), VOL, EDGFNC, a, b, c,xx(4),yy(4),zz(4),
 
255
     :     L(6), det, getdtm, angfnc
 
256
c
 
257
      logical negate, tstdts
 
258
c
 
259
      tstdts = .false.
 
260
c
 
261
      negate = .false.
 
262
C
 
263
      ELMFNC = 1E+30
 
264
      UNDLEN = 0.0
 
265
C
 
266
      IF( IERR .NE. 0 ) RETURN
 
267
C
 
268
C - Sort the nodes into ascending order
 
269
C
 
270
      NODS(1) = NOD1
 
271
      NODS(2) = NOD2
 
272
      NODS(3) = NOD3
 
273
      NODS(4) = NOD4
 
274
C
 
275
      K = NODS(1)
 
276
      IF( K .GT. NODS(2) ) THEN
 
277
        NODS(1) = NODS(2)
 
278
        NODS(2) = K
 
279
        NEGATE = .NOT. NEGATE
 
280
      ELSE
 
281
        K = NODS(2)
 
282
      END IF
 
283
      IF( K .GT. NODS(3) ) THEN
 
284
        NODS(2) = NODS(3)
 
285
        NODS(3) = K
 
286
        NEGATE = .NOT. NEGATE
 
287
      ELSE
 
288
        K = NODS(3)
 
289
      END IF
 
290
      IF( K .GT. NODS(4) ) THEN
 
291
        NODS(3) = NODS(4)
 
292
        NODS(4) = K
 
293
        NEGATE = .NOT. NEGATE
 
294
      END IF
 
295
      K = NODS(1)
 
296
      IF( K .GT. NODS(2) ) THEN
 
297
        NODS(1) = NODS(2)
 
298
        NODS(2) = K
 
299
        NEGATE = .NOT. NEGATE
 
300
      ELSE
 
301
        K = NODS(2)
 
302
      END IF
 
303
      IF( K .GT. NODS(3) ) THEN
 
304
        NODS(2) = NODS(3)
 
305
        NODS(3) = K
 
306
        NEGATE = .NOT. NEGATE
 
307
      END IF
 
308
      K = NODS(1)
 
309
      IF( K .GT. NODS(2) ) THEN
 
310
        NODS(1) = NODS(2)
 
311
        NODS(2) = K
 
312
        NEGATE = .NOT. NEGATE
 
313
      END IF
 
314
C
 
315
cc      if( debug ) print*,'   sorted: ',nods
 
316
C
 
317
C - Read the in-sphere radius stored in BIGLST
 
318
C
 
319
      INTVAL = 0
 
320
      NXT1   = 0
 
321
      IF( IELM .GT. 0 ) THEN
 
322
         INTVAL = BIGLST(7,IELM)
 
323
         NXT1   = BIGLST(NXTBIG,IELM)
 
324
      END IF
 
325
C
 
326
C - If there is no stored in-sphere info, then calculate it
 
327
C
 
328
      IF( INTVAL .EQ. 0 ) THEN
 
329
C
 
330
         DO J = 1, 9
 
331
            M(J) = 0.0
 
332
         END DO
 
333
C
 
334
         DO I = 1, 4
 
335
            if( tstdts .or. debug ) then
 
336
               det = getdtm( nodlst(7,nods(i)) )
 
337
               if( det .lt. 1e-20 ) then
 
338
                  print*,'+++ ELMFNC: Warning, small det:', det,nods(i)
 
339
                  print*,(nodlst(j,nods(i)),j=7,9)
 
340
                  print*,(nodlst(j,nods(i)),j=10,12)
 
341
                  print*,(nodlst(j,nods(i)),j=13,15)
 
342
               end if
 
343
            end if
 
344
            X(I) = NODLST(1,NODS(I))
 
345
            Y(I) = NODLST(2,NODS(I))
 
346
            Z(I) = NODLST(3,NODS(I))
 
347
            DO J = 1, 9
 
348
               M(J) = M(J) + NODLST(J+6,NODS(I))*0.25
 
349
            END DO
 
350
         END DO
 
351
C
 
352
C - get all the element info with respect to metric M
 
353
C - (Volume, Face areas, Edge lengths, In-sphere, Quality)
 
354
C
 
355
         CALL MTETIN( X, Y, Z, M, VOL, AREAS, L, UNDLEN, Q )
 
356
C
 
357
cc         if( debug ) print*,'   vol,rad,qly: ',vol,undlen,q
 
358
C
 
359
C - If we are using the quality of the element (instead of the in-sphere)...
 
360
C
 
361
         if( useqly ) undlen = 1.0/q/bst
 
362
C
 
363
         D = UNDLEN*BST
 
364
C
 
365
C - Don't worry about the maximum in-sphere radius, since such large
 
366
C - sizes will not contribute anything significant to the functional
 
367
C
 
368
            D = MIN(MAX(D,-2.01E+3),2.01E+3)
 
369
C
 
370
            INTVAL = INT(D*1000000)
 
371
            IF( ABS(INTVAL) .LE. 1 ) THEN
 
372
              IF( INTVAL .EQ. 0 .AND. IELM .LE. 0 ) THEN
 
373
                ELMFNC = 1E+32
 
374
                IF( NEGATE ) UNDLEN = -UNDLEN
 
375
                RETURN
 
376
              END IF
 
377
              IF( D .GT. 0.0 ) THEN
 
378
                INTVAL = 1
 
379
              ELSE IF( D .LT. 0.0 ) THEN
 
380
                INTVAL = -1
 
381
              ELSE
 
382
                write(0,*) 
 
383
     :                '*** ELMFNC: GOT ZERO IN-SPHERE!',IELM
 
384
                IF( IELM .LE. 0 ) write(0,*)
 
385
     :                '  Nodes: ',NOD1,NOD2,NOD3,NOD4
 
386
              END IF
 
387
              IF( IELM .GT. 0 .AND. INTVAL .EQ. 0 ) THEN
 
388
                write(0,*)
 
389
     :               '+++ Existing elem with zero in-sphere:',ielm
 
390
                write(0,*)
 
391
     :               '  d,v,r,q:',d,vol,undlen,q
 
392
                write(0,*)
 
393
     :               '  node 1: ',(nodlst(j,nod1),j=1,3)
 
394
                write(0,*)
 
395
     :               '  node 2: ',(nodlst(j,nod2),j=1,3)
 
396
                write(0,*)
 
397
     :               '  node 3: ',(nodlst(j,nod3),j=1,3)
 
398
                write(0,*)
 
399
     :               '  node 4: ',(nodlst(j,nod4),j=1,3)
 
400
                write(0,*)
 
401
     :               '  Element metric:'
 
402
                write(0,*)(m(i),i=1,3)
 
403
                write(0,*)(m(i),i=4,6)
 
404
                write(0,*)(m(i),i=7,9)
 
405
                do i = 1, 9
 
406
                  m(i) = 0.0
 
407
                end do
 
408
                m(1) = 1.0
 
409
                m(5) = 1.0
 
410
                m(9) = 1.0
 
411
                write(0,*)
 
412
     :               '  Calculating values in flat space...'
 
413
                call mtetin( x, y, z, m, vol, areas, l, undlen, q )
 
414
                write(0,*)
 
415
     :               '    r,v,q: ',undlen,vol,q
 
416
                IERR = -676
 
417
                RETURN
 
418
              END IF
 
419
            END IF
 
420
            IF( IELM .GT. 0 ) BIGLST(7,IELM) = INTVAL
 
421
            D = FLOAT(INTVAL)/1000000
 
422
C
 
423
      ELSE
 
424
C
 
425
         D = FLOAT(INTVAL)/1000000
 
426
C
 
427
      END IF
 
428
C
 
429
C - work out the *signed* in-sphere radius (or quality) to be returned
 
430
C - (this is based on the node numbering originally sent into the call
 
431
C - and is relative to the value calculated when the nodes are in ascending
 
432
C - order in the calculation).
 
433
C
 
434
      UNDLEN = D/BST
 
435
      IF( NEGATE ) UNDLEN = -UNDLEN
 
436
C
 
437
C - Calculate the in-sphere contribution to the element functional.
 
438
C - This is 'fiddled' so it gets large extremely quickly once the
 
439
C - in-sphere becomes very small. Note that the value is continuous,
 
440
C - even though different equations are used for different ranges of
 
441
C - the in-sphere radius.
 
442
C
 
443
      D = ABS(D)
 
444
      IF( D .GT. 9.9E-6 ) THEN
 
445
         D = 1.0/D - 1.0
 
446
         D = D*D
 
447
      ELSE IF( D .GT. 9.9E-7 ) THEN
 
448
         D = D*23713.74
 
449
         D = 1.0/D
 
450
         D = D*D*D*D
 
451
         D = D*D*D*D
 
452
      ELSE
 
453
         D = 4E+26
 
454
      END IF
 
455
C
 
456
      EDGON = .NOT. EDGON
 
457
C
 
458
      IF( EDGON ) THEN
 
459
C
 
460
         INTVAL = 0
 
461
C
 
462
C - retreive the edge functional info (if NXT1>0)
 
463
C
 
464
         IF( NXT1 .GT. 0 ) INTVAL = BIGLST(7,NXT1)
 
465
C
 
466
         IF( INTVAL .EQ. 0 ) THEN
 
467
C
 
468
C - calculate (and store, if NXT1>0) half the sum of edge functionals
 
469
C
 
470
c           L(1) = 1.0 - L(1)
 
471
c           L(2) = 1.0 - L(2)
 
472
c           L(3) = 1.0 - L(3)
 
473
c           L(4) = 1.0 - L(4)
 
474
c           L(5) = 1.0 - L(5)
 
475
c           L(6) = 1.0 - L(6)
 
476
c           D1   = L(1)*L(1) + L(2)*L(2) + L(3)*L(3)
 
477
c     :          + L(4)*L(4) + L(5)*L(5) + L(6)*L(6)
 
478
c
 
479
           D1 = 0.5 * ( EDGFNC( BIGLST, NODLST, NOD1, NOD2, D2 )
 
480
     :                + EDGFNC( BIGLST, NODLST, NOD1, NOD3, D3 )
 
481
     :                + EDGFNC( BIGLST, NODLST, NOD1, NOD4, D2 )
 
482
     :                + EDGFNC( BIGLST, NODLST, NOD2, NOD3, D3 )
 
483
     :                + EDGFNC( BIGLST, NODLST, NOD2, NOD4, D2 )
 
484
     :                + EDGFNC( BIGLST, NODLST, NOD3, NOD4, D3 ) )
 
485
C
 
486
           IF( ANGLEFACT .GT. 0.0 )
 
487
     :         D1 = D1 + ANGFNC( BIGLST, NODLST, NODS )*ANGLEFACT
 
488
C
 
489
C - large edges will be close to the functional bound, so we can cap them
 
490
C
 
491
           D1 = MIN(D1,2.01E+6)
 
492
           INTVAL = INT(D1*1000)
 
493
           IF( INTVAL .EQ. 0 ) INTVAL = 1
 
494
           IF( NXT1 .GT. 0 ) BIGLST(7,NXT1) = INTVAL
 
495
           D1 = FLOAT(INTVAL)/1000
 
496
C
 
497
         ELSE
 
498
           D1 = FLOAT(INTVAL)/1000
 
499
         END IF
 
500
C
 
501
      ELSE
 
502
         D1 = 0.0
 
503
      END IF
 
504
C
 
505
      EDGON = .NOT. EDGON
 
506
C
 
507
C - This is an attempt to fix a problem where the in-sphere radius
 
508
C - contributes too much to the functional for elements which have
 
509
C - some kind of constraint to their edges (probably through the
 
510
C - geometry), meaning the edges end up becoming far too large (as it
 
511
C - tries to increase the in-sphere radius by any means possible).
 
512
C - The basic strategy to combat this is to increase the contribution
 
513
C - from the edges as they become too large.
 
514
C - It's worth remembering that having long edges is actually far worse
 
515
C - for interpolation errors than having small/thin elements!
 
516
C - Note that the value of the functional contribution is continuous,
 
517
C - even though different equations are used for ranges of d1...
 
518
C
 
519
      if( d1 .gt. 2E+6 ) then
 
520
         d1 = 1.6E+25
 
521
c         d1 = 3.2E+26
 
522
c      else if( d1 .gt. 10 ) then
 
523
c         d1 = d1 + d1*d1*d1*d1*d1 - 91102.0
 
524
      else if( d1 .gt. 6 ) then
 
525
         d1 = d1 + d1*d1*d1*d1 - 1102.0
 
526
      else if( d1 .gt. 3 ) then
 
527
         d1 = d1 + d1*d1*d1 - 22.0
 
528
      else if( d1 .gt. 2 ) then
 
529
         d1 = d1 + d1*d1 - 4.0
 
530
      end if
 
531
C
 
532
C - Work out the entire element functional from half the sum of the edge
 
533
C - functionals plus the in-sphere functional.
 
534
C
 
535
      ELMFNC = D1 + D
 
536
C
 
537
      RETURN
 
538
      END
 
539
C
 
540
C
 
541
C
 
542
      REAL FUNCTION ANGFNC( BIGLST, NODLST, NODS )
 
543
C
 
544
      IMPLICIT NONE
 
545
C
 
546
      INTEGER NODS(4)
 
547
C
 
548
      INCLUDE 'blknew.i'
 
549
C
 
550
      INCLUDE 'blkbig.i'
 
551
C
 
552
      REAL NX1, NY1, NZ1, NX2, NY2, NZ2, NX3, NY3, NZ3, NX4, NY4, NZ4,
 
553
     :     X(6), Y(6), Z(6), V, D
 
554
C
 
555
      INTEGER I, J, K
 
556
C
 
557
      angfnc = 1E+30
 
558
C
 
559
      NX1 = NODLST(1,NODS(1))
 
560
      NY1 = NODLST(2,NODS(1))
 
561
      NZ1 = NODLST(3,NODS(1))
 
562
C
 
563
      DO I = 1, 3
 
564
         X(I) = NODLST(1,NODS(I+1)) - NX1
 
565
         Y(I) = NODLST(2,NODS(I+1)) - NY1
 
566
         Z(I) = NODLST(3,NODS(I+1)) - NZ1
 
567
      END DO
 
568
C
 
569
      NX1 = NODLST(1,NODS(2))
 
570
      NY1 = NODLST(2,NODS(2))
 
571
      NZ1 = NODLST(3,NODS(2))
 
572
C
 
573
      DO I = 4, 5
 
574
         X(I) = NODLST(1,NODS(I-1)) - NX1
 
575
         Y(I) = NODLST(2,NODS(I-1)) - NY1
 
576
         Z(I) = NODLST(3,NODS(I-1)) - NZ1
 
577
      END DO
 
578
C
 
579
      X(6) = NODLST(1,NODS(4)) - NODLST(1,NODS(3))
 
580
      Y(6) = NODLST(2,NODS(4)) - NODLST(2,NODS(3))
 
581
      Z(6) = NODLST(3,NODS(4)) - NODLST(3,NODS(3))
 
582
C
 
583
C - work out normals to faces
 
584
C
 
585
      NX1 = Y(2)*Z(4) - Y(4)*Z(2)
 
586
      NY1 = Z(2)*X(4) - Z(4)*X(2)
 
587
      NZ1 = X(2)*Y(4) - X(4)*Y(2)
 
588
C
 
589
      V = NX1*NX1 + NY1*NY1 + NZ1*NZ1
 
590
      if( v .lt. 1e-23 ) then
 
591
         return
 
592
         print*,'got small normal 1: ',v
 
593
         print*,nx1,ny1,nz1
 
594
         print*,nods
 
595
      else
 
596
         V = 1.0/SQRT(V)
 
597
         NX1 = NX1*V
 
598
         NY1 = NY1*V
 
599
         NZ1 = NZ1*V
 
600
      end if
 
601
C
 
602
      NX2 = Y(3)*Z(1) - Y(1)*Z(3)
 
603
      NY2 = Z(3)*X(1) - Z(1)*X(3)
 
604
      NZ2 = X(3)*Y(1) - X(1)*Y(3)
 
605
C
 
606
      V = NX2*NX2 + NY2*NY2 + NZ2*NZ2
 
607
      if( v .lt. 1e-23 ) then
 
608
         return
 
609
         print*,'got small normal 2: ',v
 
610
         print*,nx2,ny2,nz2
 
611
         print*,nods
 
612
      else
 
613
         V = 1.0/SQRT(V)
 
614
         NX2 = NX2*V
 
615
         NY2 = NY2*V
 
616
         NZ2 = NZ2*V
 
617
      end if
 
618
C
 
619
      NX3 = Y(5)*Z(4) - Y(4)*Z(5)
 
620
      NY3 = Z(5)*X(4) - Z(4)*X(5)
 
621
      NZ3 = X(5)*Y(4) - X(4)*Y(5)
 
622
C
 
623
      V = NX3*NX3 + NY3*NY3 + NZ3*NZ3
 
624
      if( v .lt. 1e-23 ) then
 
625
         return
 
626
         print*,'got small normal 3: ',v
 
627
         print*,nx3,ny3,nz3
 
628
         print*,nods
 
629
      else
 
630
         V = 1.0/SQRT(V)
 
631
         NX3 = NX3*V
 
632
         NY3 = NY3*V
 
633
         NZ3 = NZ3*V
 
634
      end if
 
635
C
 
636
      NX4 = Y(3)*Z(6) - Y(6)*Z(3)
 
637
      NY4 = Z(3)*X(6) - Z(6)*X(3)
 
638
      NZ4 = X(3)*Y(6) - X(6)*Y(3)
 
639
C
 
640
      V = NX4*NX4 + NY4*NY4 + NZ4*NZ4
 
641
      if( v .lt. 1e-23 ) then
 
642
         return
 
643
         print*,'got small normal 4: ',v
 
644
         print*,nx4,ny4,nz4
 
645
         print*,nods
 
646
      else
 
647
         V = 1.0/SQRT(V)
 
648
         NX4 = NX4*V
 
649
         NY4 = NY4*V
 
650
         NZ4 = NZ4*V
 
651
      end if
 
652
C
 
653
C - work out the functionals for the angle between each pair of faces
 
654
C
 
655
      D = NX1*NX2 + NY1*NY2 + NZ1*NZ2 + 1.0
 
656
      IF( D .GT. 1E-15 ) THEN
 
657
         D = ( 1.0/D - 1.0 )
 
658
      ELSE
 
659
         return
 
660
         print*,'got d = ',d,' for normals 1+2...'
 
661
         print*,nx1,ny1,nz1
 
662
         print*,nx2,ny2,nz2
 
663
         print*,x
 
664
         print*,y
 
665
         print*,z
 
666
         D = 1E+15
 
667
      END IF
 
668
C
 
669
      V = D*D
 
670
C
 
671
      D = NX1*NX3 + NY1*NY3 + NZ1*NZ3 + 1.0
 
672
      IF( D .GT. 1E-15 ) THEN
 
673
         D = ( 1.0/D - 1.0 )
 
674
      ELSE
 
675
         return
 
676
         print*,'got d = ',d,' for normals 1+3...'
 
677
         print*,nx1,ny1,nz1
 
678
         print*,nx3,ny3,nz3
 
679
         print*,x
 
680
         print*,y
 
681
         print*,z
 
682
         D = 1E+15
 
683
      END IF
 
684
C
 
685
      V = V + D*D
 
686
C
 
687
      D = NX1*NX4 + NY1*NY4 + NZ1*NZ4 + 1.0
 
688
      IF( D .GT. 1E-15 ) THEN
 
689
         D = ( 1.0/D - 1.0 )
 
690
      ELSE
 
691
         return
 
692
         print*,'got d = ',d,' for normals 1+4...'
 
693
         print*,nx1,ny1,nz1
 
694
         print*,nx4,ny4,nz4
 
695
         print*,x
 
696
         print*,y
 
697
         print*,z
 
698
         D = 1E+15
 
699
      END IF
 
700
C
 
701
      V = V + D*D
 
702
C
 
703
      D = NX2*NX3 + NY2*NY3 + NZ2*NZ3 + 1.0
 
704
      IF( D .GT. 1E-15 ) THEN
 
705
         D = ( 1.0/D - 1.0 )
 
706
      ELSE
 
707
         return
 
708
         print*,'got d = ',d,' for normals 2+3...'
 
709
         print*,nx2,ny2,nz2
 
710
         print*,nx3,ny3,nz3
 
711
         print*,x
 
712
         print*,y
 
713
         print*,z
 
714
         D = 1E+15
 
715
      END IF
 
716
C
 
717
      V = V + D*D
 
718
C
 
719
      D = NX2*NX4 + NY2*NY4 + NZ2*NZ4 + 1.0
 
720
      IF( D .GT. 1E-15 ) THEN
 
721
         D = ( 1.0/D - 1.0 )
 
722
      ELSE
 
723
         return
 
724
         print*,'got d = ',d,' for normals 2+4...'
 
725
         print*,nx2,ny2,nz2
 
726
         print*,nx4,ny4,nz4
 
727
         print*,x
 
728
         print*,y
 
729
         print*,z
 
730
         D = 1E+15
 
731
      END IF
 
732
C
 
733
      V = V + D*D
 
734
C
 
735
      D = NX3*NX4 + NY3*NY4 + NZ3*NZ4 + 1.0
 
736
      IF( D .GT. 1E-15 ) THEN
 
737
         D = ( 1.0/D - 1.0 )
 
738
      ELSE
 
739
         return
 
740
         print*,'got d = ',d,' for normals 3+4...'
 
741
         print*,nx3,ny3,nz3
 
742
         print*,nx4,ny4,nz4
 
743
         print*,x
 
744
         print*,y
 
745
         print*,z
 
746
         D = 1E+15
 
747
      END IF
 
748
C
 
749
      ANGFNC = (V + D*D)/6
 
750
C
 
751
      RETURN
 
752
      END
 
753
C
 
754
#endif