Isles.c 99.6 KB
Newer Older
1
2
3
4
5
6
/*
(c) 2005 Supratik Mukhopadhyay, Nayana Majumdar
*/

#define DEFINE_ISLESGLOBAL

7
8
#include "Isles.h"

Heinrich Schindler's avatar
Heinrich Schindler committed
9
10
#include <complex.h>
#include <math.h>
11
12
#include <stdio.h>
#include <stdlib.h>
13

14
#include <gsl/gsl_sf.h>
15
16

#define SHIFT 2.0
17
18
19
20
21
22
23
24
25
#define ARMAX 10000.0  // Maximum aspect ratio for an element
#define ARMIN 0.0001   // Minimum aspect ratio for an element

// #define XNSegApprox 10	// much less time but inaccurate results
// #define ZNSegApprox 10	// much less time but inaccurate results
#define XNSegApprox 100  // barely acceptable results
#define ZNSegApprox 100  // barely acceptable result
// #define XNSegApprox 1000	// much better results but lot more time
// #define ZNSegApprox 1000	// much better results but lot more time
26

27
#ifdef __cplusplus
Heinrich Schindler's avatar
Heinrich Schindler committed
28
#include <cmath>
29
30
#ifdef _GLIBCXX_HAVE_OBSOLETE_ISINF_ISNAN
using std::isinf;
31
using std::isnan;
32
#endif
33
34
35
namespace neBEM {
#endif

36
37
38
39
40
41
// Exact potential and flux for a rectangular element
// Potential needs all the calculations, whereas, the fluxes need part of the
// information generated during the computation for potential. This is the
// reason why we have included flux computation along with the potential
// computation.
// Expressions from:
Heinrich Schindler's avatar
Heinrich Schindler committed
42
43
int ExactRecSurf(double X, double Y, double Z, double xlo, double zlo,
                 double xhi, double zhi, double *Potential, Vector3D *Flux) {
Heinrich Schindler's avatar
Heinrich Schindler committed
44
  if (DebugISLES) printf("In ExactRecSurf ...\n");
Heinrich Schindler's avatar
Heinrich Schindler committed
45
46
47
48
49
50
51
52
53
54
55

  ++IslesCntr;
  ++ExactCntr;

  ApproxFlag = 0;

  if ((fabs(xhi - xlo) < 3.0 * MINDIST) || (fabs(zhi - zlo) < 3.0 * MINDIST)) {
    fprintf(stdout, "Element size too small! ... returning ...\n");
    return -1;
  }

56
57
  double ar = fabs((zhi - zlo) / (xhi - xlo));
  if (ar > ARMAX || ar < ARMIN) {
Heinrich Schindler's avatar
Heinrich Schindler committed
58
59
60
61
62
63
64
65
    fprintf(stdout, "Element too thin! ... returning ...\n");
    return -2;
  }

  if (fabs(X) < MINDIST) X = 0.0;
  if (fabs(Y) < MINDIST) Y = 0.0;
  if (fabs(Z) < MINDIST) Z = 0.0;

66
  double dxlo = X - xlo;  // zero at the X=xlo edge
Heinrich Schindler's avatar
Heinrich Schindler committed
67
  if (fabs(dxlo) < MINDIST) dxlo = 0.0;
68
  double dxhi = X - xhi;  // zero at the X=xhi edge
Heinrich Schindler's avatar
Heinrich Schindler committed
69
  if (fabs(dxhi) < MINDIST) dxhi = 0.0;
70
  double dzlo = Z - zlo;  // zero at the Z=zlo edge
Heinrich Schindler's avatar
Heinrich Schindler committed
71
  if (fabs(dzlo) < MINDIST) dzlo = 0.0;
72
  double dzhi = Z - zhi;  // zero at the Z=zhi edge
Heinrich Schindler's avatar
Heinrich Schindler committed
73
74
75
76
77
  if (fabs(dzhi) < MINDIST) dzhi = 0.0;

  // These four parameters can never be zero except at the four corners where
  // one of them can become zero. For example, at X=xlo, Y=0, Z=zlo, D11
  // is zero but the others are nonzero.
78
  double D11 = sqrt(dxlo * dxlo + Y * Y + dzlo * dzlo);
Heinrich Schindler's avatar
Heinrich Schindler committed
79
  if (fabs(D11) < MINDIST) D11 = 0.0;
80
  double D21 = sqrt(dxhi * dxhi + Y * Y + dzlo * dzlo);
Heinrich Schindler's avatar
Heinrich Schindler committed
81
  if (fabs(D21) < MINDIST) D21 = 0.0;
82
  double D12 = sqrt(dxlo * dxlo + Y * Y + dzhi * dzhi);
Heinrich Schindler's avatar
Heinrich Schindler committed
83
  if (fabs(D12) < MINDIST) D12 = 0.0;
84
  double D22 = sqrt(dxhi * dxhi + Y * Y + dzhi * dzhi);
Heinrich Schindler's avatar
Heinrich Schindler committed
85
86
87
  if (fabs(D22) < MINDIST) D22 = 0.0;

  // Parameters related to the Y terms
88
89
90
91
92
93
94
95
  int S1 = Sign(dzlo);
  int S2 = Sign(dzhi);
  double modY = fabs(Y);
  int SY = Sign(Y);
  double I1 = dxlo * modY;
  double I2 = dxhi * modY;
  double R1 = Y * Y + dzlo * dzlo;
  double R2 = Y * Y + dzhi * dzhi;
Heinrich Schindler's avatar
Heinrich Schindler committed
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
  if (fabs(I1) < MINDIST2) I1 = 0.0;
  if (fabs(I2) < MINDIST2) I2 = 0.0;
  if (fabs(R1) < MINDIST2) R1 = 0.0;
  if (fabs(R2) < MINDIST2) R2 = 0.0;

  if (DebugISLES) {
    fprintf(stdout, "X: %.16lg, Y: %.16lg, Z: %.16lg\n", X, Y, Z);
    fprintf(stdout, "xlo: %.16lg, zlo: %.16lg, xhi: %.16lg, zhi: %.16lg\n", xlo,
            zlo, xhi, zhi);
    fprintf(stdout, "dxlo: %.16lg, dzlo: %.16lg, dxhi: %.16lg, dzhi: %.16lg\n",
            dxlo, dzlo, dxhi, dzhi);
    fprintf(stdout, "D11: %.16lg, D12: %.16lg, D21: %.16lg, D22: %.16lg\n", D11,
            D12, D21, D22);
    fprintf(stdout, "S1: %d, S2: %d, modY: %.16lg\n", S1, S2, modY);
    fprintf(stdout, "I1: %.16lg, I2: %.16lg, R1: %.16lg, R2: %.16lg\n", I1, I2,
            R1, R2);
    fprintf(stdout, "MINDIST: %.16lg, MINDIST2: %.16lg, SHIFT: %.16lg\n",
            MINDIST, MINDIST2, SHIFT);
    fflush(stdout);
  }

  // Check for possible numerical difficuties and take care.
  // Presently the idea is to shift the field point slightly to a 'safe'
  // position. Note that a shift in Y does not work because the singularities
  // are associated with D11-s and dxlo-s and their combinations. A Y shift can
  // sometimes alleviate the problem, but it cannot gurantee a permanent1s
  // solution.
  // A better approach is to evaluate analytic expressions truly valid for
  // these critical regions, especially to take care of the >= and <=
  // comparisons. For the == case, both of the previous comparisons are true and
  // it is hard to justify one choice over the other. On the other hand, there
  // are closed form analytic solutions for the == cases, and the problem does
  // not stem from round-off errors as in the > or <, but not quite == cases.
  // Note that shifts that are exactly equal to MINDIST, can give rise to
  // un-ending recursions, leading to Segmentation fault.
  // Corners
  // Four corners
  // Averages over four point evaluations may be carried out (skipped at
  // present) This, however, may be difficult since we'll have to make sure that
  // the shift towards the element does not bring the point too close to the
  // same, or another difficult-to-evaluate situation. One of the ways to ensure
  // this is to make SHIFT large enough, but that is unreasnoable and will
  // introduce large amount of error.
139
140
  if ((fabs(D11) <= MINDIST)) {
    // close to xlo, 0, zlo
Heinrich Schindler's avatar
Heinrich Schindler committed
141
    if (DebugISLES) printf("fabs(D11) <= MINDIST ... ");
142
143
    double X1 = X;
    double Z1 = Z;
144
145
    if ((X >= xlo) && (Z >= zlo)) {
      // point on the element
Heinrich Schindler's avatar
Heinrich Schindler committed
146
      if (DebugISLES) printf("Case 1\n");
147
148
      X1 = X + SHIFT * MINDIST;
      Z1 = Z + SHIFT * MINDIST;
149
150
    } else if ((X <= xlo) && (Z >= zlo)) {
      // field point outside the element
Heinrich Schindler's avatar
Heinrich Schindler committed
151
      if (DebugISLES) printf("Case 2 ...\n");
152
153
      X1 = X - SHIFT * MINDIST;
      Z1 = Z + SHIFT * MINDIST;
154
155
    } else if ((X >= xlo) && (Z <= zlo)) {
      // field point outside the element
Heinrich Schindler's avatar
Heinrich Schindler committed
156
      if (DebugISLES) printf("Case 3 ...\n");
157
158
      X1 = X + SHIFT * MINDIST;
      Z1 = Z - SHIFT * MINDIST;
159
160
    } else if ((X <= xlo) && (Z <= zlo)) {
      // field point outside the element
Heinrich Schindler's avatar
Heinrich Schindler committed
161
      if (DebugISLES) printf("Case 4 ...\n");
162
163
      X1 = X - SHIFT * MINDIST;
      Z1 = Z - SHIFT * MINDIST;
Heinrich Schindler's avatar
Heinrich Schindler committed
164
    }
165
166
167
168
169
170
171
172
    double Pot1;
    Vector3D Flux1;
    ExactRecSurf(X1, Y, Z1, xlo, zlo, xhi, zhi, &Pot1, &Flux1);
    *Potential = Pot1;
    Flux->X = Flux1.X;
    Flux->Y = Flux1.Y;
    Flux->Z = Flux1.Z;
    return 0;
Heinrich Schindler's avatar
Heinrich Schindler committed
173
  }
174
175
  if ((fabs(D21) <= MINDIST)) {
    // close to xhi, 0, zlo
Heinrich Schindler's avatar
Heinrich Schindler committed
176
    if (DebugISLES) printf("fabs(D21) <= MINDIST ... ");
177
178
    double X1 = X;
    double Z1 = Z;
179
180
    if ((X >= xhi) && (Z >= zlo)) {
      // point outside the element
Heinrich Schindler's avatar
Heinrich Schindler committed
181
      if (DebugISLES) printf("Case 1 ...\n");
182
183
      X1 = X + SHIFT * MINDIST;
      Z1 = Z + SHIFT * MINDIST;
184
185
    } else if ((X <= xhi) && (Z >= zlo)) {
      // point on the element
Heinrich Schindler's avatar
Heinrich Schindler committed
186
      if (DebugISLES) printf("Case 2 ...\n");
187
188
      X1 = X - SHIFT * MINDIST;
      Z1 = Z + SHIFT * MINDIST;
189
190
    } else if ((X >= xhi) && (Z <= zlo)) {
      // field point outside the element
Heinrich Schindler's avatar
Heinrich Schindler committed
191
      if (DebugISLES) printf("Case 3 ...\n");
192
193
      X1 = X + SHIFT * MINDIST;
      Z1 = Z - SHIFT * MINDIST;
194
195
    } else if ((X <= xhi) && (Z <= zlo)) {
      // field point outside the element
Heinrich Schindler's avatar
Heinrich Schindler committed
196
      if (DebugISLES) printf("Case 4 ...\n");
197
198
      X1 = X - SHIFT * MINDIST;
      Z1 = Z - SHIFT * MINDIST;
Heinrich Schindler's avatar
Heinrich Schindler committed
199
    }
200
201
202
203
204
205
206
207
    double Pot1;
    Vector3D Flux1;
    ExactRecSurf(X1, Y, Z1, xlo, zlo, xhi, zhi, &Pot1, &Flux1);
    *Potential = Pot1;
    Flux->X = Flux1.X;
    Flux->Y = Flux1.Y;
    Flux->Z = Flux1.Z;
    return 0;
Heinrich Schindler's avatar
Heinrich Schindler committed
208
  }
209
  if ((fabs(D12) <= MINDIST)) {
210
    // close to xlo, 0, zhi
Heinrich Schindler's avatar
Heinrich Schindler committed
211
    if (DebugISLES) printf("fabs(D12) <= MINDIST ... ");
212
213
    double X1 = X;
    double Z1 = Z;
214
215
    if ((X >= xlo) && (Z >= zhi)) {
      // point outside the element
Heinrich Schindler's avatar
Heinrich Schindler committed
216
      if (DebugISLES) printf("Case 1 ...\n");
217
218
      X1 = X + SHIFT * MINDIST;
      Z1 = Z + SHIFT * MINDIST;
219
220
    } else if ((X <= xlo) && (Z >= zhi)) {
      // field point outside the element
Heinrich Schindler's avatar
Heinrich Schindler committed
221
      if (DebugISLES) printf("Case 2 ...\n");
222
223
      X1 = X - SHIFT * MINDIST;
      Z1 = Z + SHIFT * MINDIST;
224
225
    } else if ((X >= xlo) && (Z <= zhi)) {
      // field point on the element
Heinrich Schindler's avatar
Heinrich Schindler committed
226
      if (DebugISLES) printf("Case 3 ...\n");
227
228
      X1 = X + SHIFT * MINDIST;
      Z1 = Z - SHIFT * MINDIST;
229
230
    } else if ((X <= xlo) && (Z <= zhi)) {
      // field point outside the element
Heinrich Schindler's avatar
Heinrich Schindler committed
231
      if (DebugISLES) printf("Case 4 ...\n");
232
233
      X1 = X - SHIFT * MINDIST;
      Z1 = Z - SHIFT * MINDIST;
Heinrich Schindler's avatar
Heinrich Schindler committed
234
    }
235
236
237
238
239
240
241
242
    double Pot1;
    Vector3D Flux1;
    ExactRecSurf(X1, Y, Z1, xlo, zlo, xhi, zhi, &Pot1, &Flux1);
    *Potential = Pot1;
    Flux->X = Flux1.X;
    Flux->Y = Flux1.Y;
    Flux->Z = Flux1.Z;
    return 0;
Heinrich Schindler's avatar
Heinrich Schindler committed
243
  }
244
245
  if ((fabs(D22) <= MINDIST)) {
    // close to xhi, 0, zhi
Heinrich Schindler's avatar
Heinrich Schindler committed
246
    if (DebugISLES) printf("fabs(D22) <= MINDIST ... ");
247
248
    double X1 = X;
    double Z1 = Z;
249
250
    if ((X >= xhi) && (Z >= zhi)) {
      // point outside the element
Heinrich Schindler's avatar
Heinrich Schindler committed
251
      if (DebugISLES) printf("Case 1 ...\n");
252
253
      X1 = X + SHIFT * MINDIST;
      Z1 = Z + SHIFT * MINDIST;
254
255
    } else if ((X <= xhi) && (Z >= zhi)) {
      // field point outside the element
Heinrich Schindler's avatar
Heinrich Schindler committed
256
      if (DebugISLES) printf("Case 2 ...\n");
257
258
      X1 = X - SHIFT * MINDIST;
      Z1 = Z + SHIFT * MINDIST;
259
260
    } else if ((X >= xhi) && (Z <= zhi)) {
      // field point outside the element
Heinrich Schindler's avatar
Heinrich Schindler committed
261
      if (DebugISLES) printf("Case 3 ...\n");
262
263
      X1 = X + SHIFT * MINDIST;
      Z1 = Z - SHIFT * MINDIST;
264
265
    } else if ((X <= xhi) && (Z <= zhi)) {
      // field point on the element
Heinrich Schindler's avatar
Heinrich Schindler committed
266
      if (DebugISLES) printf("Case 4 ...\n");
267
268
      X1 = X - SHIFT * MINDIST;
      Z1 = Z - SHIFT * MINDIST;
Heinrich Schindler's avatar
Heinrich Schindler committed
269
    }
270
271
272
273
274
275
276
277
    double Pot1;
    Vector3D Flux1;
    ExactRecSurf(X1, Y, Z1, xlo, zlo, xhi, zhi, &Pot1, &Flux1);
    *Potential = Pot1;
    Flux->X = Flux1.X;
    Flux->Y = Flux1.Y;
    Flux->Z = Flux1.Z;
    return 0;
Heinrich Schindler's avatar
Heinrich Schindler committed
278
279
280
281
282
  }
  // Four edges	- average over two points on two sides of the edge may be ok.
  // Here also we'll have to make sure that the shift
  // towards the element does not bring the point too close to the same, or
  // another difficult-to-evaluate situation. One of the ways to ensure this
283
  // is to make SHIFT large enough, but that is unreasonable and will introduce
Heinrich Schindler's avatar
Heinrich Schindler committed
284
  // large amount of error.
285
286
  if (fabs(dxlo) < MINDIST) {
    // edge at x=xlo || to Z - axis
Heinrich Schindler's avatar
Heinrich Schindler committed
287
    if (DebugISLES) printf("fabs(dxlo) < MINDIST ... ");
288
289
290
291
    double X1 = X;
    double X2 = X;
    if (X >= xlo) {
      // field point on +ve side of YZ plane
Heinrich Schindler's avatar
Heinrich Schindler committed
292
      if (DebugISLES) printf("Case 1 ...\n");
293
294
295
296
      X1 += SHIFT * MINDIST;
      X2 -= SHIFT * MINDIST;
    } else {
      // field point on -ve side of YZ plane
Heinrich Schindler's avatar
Heinrich Schindler committed
297
      if (DebugISLES) printf("Case 2 ...\n");
298
299
      X1 -= SHIFT * MINDIST;
      X2 += SHIFT * MINDIST;
Heinrich Schindler's avatar
Heinrich Schindler committed
300
    }
301
302
303
304
305
306
307
308
309
    double Pot1, Pot2;
    Vector3D Flux1, Flux2;
    ExactRecSurf(X1, Y, Z, xlo, zlo, xhi, zhi, &Pot1, &Flux1);
    ExactRecSurf(X2, Y, Z, xlo, zlo, xhi, zhi, &Pot2, &Flux2);
    *Potential = 0.5 * (Pot1 + Pot2);
    Flux->X = 0.5 * (Flux1.X + Flux2.X);
    Flux->Y = 0.5 * (Flux1.Y + Flux2.Y);
    Flux->Z = 0.5 * (Flux1.Z + Flux2.Z);
    return 0;
Heinrich Schindler's avatar
Heinrich Schindler committed
310
  }
311
312
  if (fabs(dzlo) < MINDIST) {
    // edge at z=zlo, || to X axis
Heinrich Schindler's avatar
Heinrich Schindler committed
313
    if (DebugISLES) printf("fabs(dzlo) < MINDIST ... ");
314
315
316
317
    double Z1 = Z;
    double Z2 = Z;
    if (Z >= zlo) {
      // field point on +ve side of XY plane
Heinrich Schindler's avatar
Heinrich Schindler committed
318
      if (DebugISLES) printf("Case 1 ...\n");
319
320
321
322
      Z1 += SHIFT * MINDIST;
      Z2 -= SHIFT * MINDIST;
    } else {
      // field point on -ve side of XY plane
Heinrich Schindler's avatar
Heinrich Schindler committed
323
      if (DebugISLES) printf("Case 2 ...\n");
324
325
      Z1 -= SHIFT * MINDIST;
      Z2 += SHIFT * MINDIST;
Heinrich Schindler's avatar
Heinrich Schindler committed
326
    }
327
328
329
330
331
332
333
334
335
    double Pot1, Pot2;
    Vector3D Flux1, Flux2;
    ExactRecSurf(X, Y, Z1, xlo, zlo, xhi, zhi, &Pot1, &Flux1);
    ExactRecSurf(X, Y, Z2, xlo, zlo, xhi, zhi, &Pot2, &Flux2);
    *Potential = 0.5 * (Pot1 + Pot2);
    Flux->X = 0.5 * (Flux1.X + Flux2.X);
    Flux->Y = 0.5 * (Flux1.Y + Flux2.Y);
    Flux->Z = 0.5 * (Flux1.Z + Flux2.Z);
    return 0;
Heinrich Schindler's avatar
Heinrich Schindler committed
336
  }
337
338
  if (fabs(dxhi) < MINDIST) {
    // edge at x=xhi, || to Z axis
Heinrich Schindler's avatar
Heinrich Schindler committed
339
    if (DebugISLES) printf("fabs(dxhi) < MINDIST ... ");
340
341
342
343
    double X1 = X;
    double X2 = X;
    if (X >= xhi) {
      // field point on +ve side of YZ plane
Heinrich Schindler's avatar
Heinrich Schindler committed
344
      if (DebugISLES) printf("Case 1 ...\n");
345
346
347
348
      X1 += SHIFT * MINDIST;
      X2 -= SHIFT * MINDIST;
    } else {
      // field point on -ve side of YZ plane
Heinrich Schindler's avatar
Heinrich Schindler committed
349
      if (DebugISLES) printf("Case 2 ...\n");
350
351
      X1 -= SHIFT * MINDIST;
      X2 += SHIFT * MINDIST;
Heinrich Schindler's avatar
Heinrich Schindler committed
352
    }
353
354
355
356
357
358
359
360
361
    double Pot1, Pot2;
    Vector3D Flux1, Flux2;
    ExactRecSurf(X1, Y, Z, xlo, zlo, xhi, zhi, &Pot1, &Flux1);
    ExactRecSurf(X2, Y, Z, xlo, zlo, xhi, zhi, &Pot2, &Flux2);
    *Potential = 0.5 * (Pot1 + Pot2);
    Flux->X = 0.5 * (Flux1.X + Flux2.X);
    Flux->Y = 0.5 * (Flux1.Y + Flux2.Y);
    Flux->Z = 0.5 * (Flux1.Z + Flux2.Z);
    return 0;
362
  }
363
364
  if (fabs(dzhi) < MINDIST) {
    // edge at z=zhi || to X axis
Heinrich Schindler's avatar
Heinrich Schindler committed
365
    if (DebugISLES) printf("fabs(dzhi) < MINDIST ... ");
366
367
368
369
    double Z1 = Z;
    double Z2 = Z;
    if (Z >= zhi) {
      // field point on +ve side of XY plane
Heinrich Schindler's avatar
Heinrich Schindler committed
370
      if (DebugISLES) printf("Case 1 ...\n");
371
372
373
374
      Z1 += SHIFT * MINDIST;
      Z2 -= SHIFT * MINDIST;
    } else {
      // field point on -ve side of XY plane
Heinrich Schindler's avatar
Heinrich Schindler committed
375
      if (DebugISLES) printf("Case 2 ...\n");
376
377
      Z1 -= SHIFT * MINDIST;
      Z2 += SHIFT * MINDIST;
Heinrich Schindler's avatar
Heinrich Schindler committed
378
    }
379
380
381
382
383
384
385
386
387
    double Pot1, Pot2;
    Vector3D Flux1, Flux2;
    ExactRecSurf(X, Y, Z1, xlo, zlo, xhi, zhi, &Pot1, &Flux1);
    ExactRecSurf(X, Y, Z2, xlo, zlo, xhi, zhi, &Pot2, &Flux2);
    *Potential = 0.5 * (Pot1 + Pot2);
    Flux->X = 0.5 * (Flux1.X + Flux2.X);
    Flux->Y = 0.5 * (Flux1.Y + Flux2.Y);
    Flux->Z = 0.5 * (Flux1.Z + Flux2.Z);
    return 0;
Heinrich Schindler's avatar
Heinrich Schindler committed
388
  }
389
390

  // Logarithmic weak singularities are possible.
391
  // Checks to be performed for 0 or -ve denominators and also
392
393
394
395
396
397
398
399
400
401
402
403
404
  // 0 and +ve numerators.
  // Interestingly, 0/0 does not cause a problem.
  double DZTerm1 = log((D11 - dzlo) / (D12 - dzhi));
  double DZTerm2 = log((D21 - dzlo) / (D22 - dzhi));
  double DXTerm1 = log((D11 - dxlo) / (D21 - dxhi));
  double DXTerm2 = log((D12 - dxlo) / (D22 - dxhi));

  if (DebugISLES) {
    fprintf(
        stdout,
        "DZTerm1: %.16lg, DZTerm2: %.16lg, DXTerm1: %.16lg, DXTerm2: %.16lg\n",
        DZTerm1, DZTerm2, DXTerm1, DXTerm2);
  }
Heinrich Schindler's avatar
Heinrich Schindler committed
405
  // Four conditions based on the logarithmic terms
406
  if (isnan(DZTerm1) || isinf(DZTerm1)) {
Heinrich Schindler's avatar
Heinrich Schindler committed
407
408
409
410
411
412
413
414
415
    ++FailureCntr;
    --ExactCntr;
    ApproxFlag = 1;
    fprintf(fIsles, "DZTerm1 problem ... approximating: %d.\n", ApproxFlag);
    if (DebugISLES)
      fprintf(stdout, "DZTerm1 problem ... approximating: %d.\n", ApproxFlag);
    return (ApproxRecSurf(X, Y, Z, xlo, zlo, xhi, zhi, XNSegApprox, ZNSegApprox,
                          Potential, Flux));
  }
416
  if (isnan(DZTerm2) || isinf(DZTerm2)) {
Heinrich Schindler's avatar
Heinrich Schindler committed
417
418
419
420
421
422
423
424
425
    ++FailureCntr;
    --ExactCntr;
    ApproxFlag = 2;
    fprintf(fIsles, "DZTerm2 problem ... approximating: %d.\n", ApproxFlag);
    if (DebugISLES)
      fprintf(stdout, "DZTerm2 problem ... approximating: %d.\n", ApproxFlag);
    return (ApproxRecSurf(X, Y, Z, xlo, zlo, xhi, zhi, XNSegApprox, ZNSegApprox,
                          Potential, Flux));
  }
426
  if (isnan(DXTerm1) || isinf(DXTerm1)) {
Heinrich Schindler's avatar
Heinrich Schindler committed
427
428
429
430
431
432
433
434
435
    ++FailureCntr;
    --ExactCntr;
    ApproxFlag = 3;
    fprintf(fIsles, "DXTerm1 problem ... approximating: %d.\n", ApproxFlag);
    if (DebugISLES)
      fprintf(stdout, "DXTerm1 problem ... approximating: %d.\n", ApproxFlag);
    return (ApproxRecSurf(X, Y, Z, xlo, zlo, xhi, zhi, XNSegApprox, ZNSegApprox,
                          Potential, Flux));
  }
436
  if (isnan(DXTerm2) || isinf(DXTerm2)) {
Heinrich Schindler's avatar
Heinrich Schindler committed
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
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
494
495
496
497
498
499
500
501
502
    ++FailureCntr;
    --ExactCntr;
    ApproxFlag = 4;
    fprintf(fIsles, "DXTerm2 problem ... approximating: %d.\n", ApproxFlag);
    if (DebugISLES)
      fprintf(stdout, "DXTerm2 problem ... approximating: %d.\n", ApproxFlag);
    return (ApproxRecSurf(X, Y, Z, xlo, zlo, xhi, zhi, XNSegApprox, ZNSegApprox,
                          Potential, Flux));
  }
  // Four criticalities based on the tanhyperbolic terms
  // Since gsl_complex_arctanh_real can have any real number as its argument,
  // all the criticalities are related to gsl_complex_arctanh where the
  // imaginary component is present. So, fabs(I1) > mindist and
  // fabs(I2) > mindist are only being tested.
  if (S1 != 0) {
    if (fabs(I1) > MINDIST2) {
      if (D11 * fabs(dzlo) < MINDIST2) {
        ++FailureCntr;
        --ExactCntr;
        ApproxFlag = 5;
        fprintf(fIsles, "S1-I1 problem ... approximating: %d.\n", ApproxFlag);
        if (DebugISLES)
          fprintf(stdout, "S1-I1 problem ... approximating: %d.\n", ApproxFlag);
        return (ApproxRecSurf(X, Y, Z, xlo, zlo, xhi, zhi, XNSegApprox,
                              ZNSegApprox, Potential, Flux));
      }
    }
    if (fabs(I2) > MINDIST2) {
      if (D21 * fabs(dzlo) < MINDIST2) {
        ++FailureCntr;
        --ExactCntr;
        ApproxFlag = 6;
        fprintf(fIsles, "S1-I2 problem ... approximating: %d.\n", ApproxFlag);
        if (DebugISLES)
          fprintf(stdout, "S1-I2 problem ... approximating: %d.\n", ApproxFlag);
        return (ApproxRecSurf(X, Y, Z, xlo, zlo, xhi, zhi, XNSegApprox,
                              ZNSegApprox, Potential, Flux));
      }
    }
  }
  if (S2 != 0) {
    if (fabs(I1) > MINDIST2) {
      if (D12 * fabs(dzhi) < MINDIST2) {
        ++FailureCntr;
        --ExactCntr;
        ApproxFlag = 7;
        fprintf(fIsles, "S2-I1 problem ... approximating: %d.\n", ApproxFlag);
        if (DebugISLES)
          fprintf(stdout, "S2-I1 problem ... approximating: %d.\n", ApproxFlag);
        return (ApproxRecSurf(X, Y, Z, xlo, zlo, xhi, zhi, XNSegApprox,
                              ZNSegApprox, Potential, Flux));
      }
    }
    if (fabs(I2) > MINDIST2) {
      if (D22 * fabs(dzhi) < MINDIST2) {
        ++FailureCntr;
        --ExactCntr;
        ApproxFlag = 8;
        fprintf(fIsles, "S2-I2 problem ... approximating: %d.\n", ApproxFlag);
        if (DebugISLES)
          fprintf(stdout, "S2-I2 problem ... approximating: %d.\n", ApproxFlag);
        return (ApproxRecSurf(X, Y, Z, xlo, zlo, xhi, zhi, XNSegApprox,
                              ZNSegApprox, Potential, Flux));
      }
    }
  }
503

504
  double sumTanTerms = 0.;
505
506
  // The possibility of singularities for dzhi or dzlo (division by zero)
  // is overridden by the fact that S1 or S2 becomes zero in such cases
Heinrich Schindler's avatar
Heinrich Schindler committed
507
  // and the singularity is avoided.
508
509
  if (S1 != 0) {
    if (fabs(I1) > MINDIST2) {
510
511
512
513
      double a = D11 * D11 * dzlo * dzlo;
      double b = R1 * R1 + I1 * I1;
      double tmp = -S1 * atan(2 * I1 * D11 * fabs(dzlo) / (a - b));
      if (b > a) {
514
        if ((X > xlo && Z > zlo) || (X < xlo && Z < zlo)) {
515
          tmp -= ST_PI;
516
        } else if ((X < xlo && Z > zlo) || (X > xlo && Z < zlo)) {
517
          tmp += ST_PI;
518
519
        }
      }
520
521
522
523
      sumTanTerms += tmp;
    }

    if (fabs(I2) > MINDIST2) {
524
525
526
527
      double a = D21 * D21 * dzlo * dzlo;
      double b = R1 * R1 + I2 * I2;
      double tmp = -S1 * atan(2 * I2 * D21 * fabs(dzlo) / (a - b));
      if (b > a) {
528
        if ((X > xhi && Z > zlo) || (X < xhi && Z < zlo)) {
529
          tmp -= ST_PI;
530
        } else if ((X < xhi && Z > zlo) || (X > xhi && Z < zlo)) {
531
          tmp += ST_PI;
532
533
        }
      }
534
535
536
      sumTanTerms -= tmp;
    }
  }
537

538
539
  if (S2 != 0) {
    if (fabs(I1) > MINDIST2) {
540
      double a = D12 * D12 * dzhi * dzhi;
541
      double b = R2 * R2 + I1 * I1;
542
543
      double tmp = -S2 * atan(2 * I1 * D12 * fabs(dzhi) / (a - b));
      if (b > a) {
544
        if ((X > xlo && Z > zhi) || (X < xlo && Z < zhi)) {
545
          tmp -= ST_PI;
546
        } else if ((X < xlo && Z > zhi) || (X > xlo && Z < zhi)) {
547
          tmp += ST_PI;
548
        }
549
      }
550
      sumTanTerms -= tmp;
Heinrich Schindler's avatar
Heinrich Schindler committed
551
    }
552

553
    if (fabs(I2) > MINDIST2) {
554
555
556
557
      double a = D22 * D22 * dzhi * dzhi;
      double b = R2 * R2 + I2 * I2;
      double tmp = -S2 * atan(2 * I2 * D22 * fabs(dzhi) / (a - b));
      if (b > a) {
558
559
560
561
562
563
564
        if ((X > xhi && Z > zhi) || (X < xhi && Z < zhi)) {
          tmp -= ST_PI;
        } else if ((X < xhi && Z > zhi) || (X > xhi && Z < zhi)) {
          tmp += ST_PI;
        }
      }
      sumTanTerms += tmp;
Heinrich Schindler's avatar
Heinrich Schindler committed
565
    }
566
  }
567

568
  sumTanTerms *= -0.5;
Heinrich Schindler's avatar
Heinrich Schindler committed
569
  double Pot = -dxlo * DZTerm1 + dxhi * DZTerm2 + modY * sumTanTerms -
570
               dzlo * DXTerm1 + dzhi * DXTerm2;
Heinrich Schindler's avatar
Heinrich Schindler committed
571
572
573
  double Fx = DZTerm1 - DZTerm2;
  double Fy = -SY * sumTanTerms;
  double Fz = DXTerm1 - DXTerm2;
574
575
576
577
578
579
580
  if (DebugISLES) {
    printf("XTerms: %.16lg, YTerms: %.16lg, ZTerms: %.16lg\n",
           -dxlo * DZTerm1 + dxhi * DZTerm2, modY * sumTanTerms,
           -dzlo * DXTerm1 + dzhi * DXTerm2);
    printf("Pot: %lf, Fx: %lf, Fy: %lf, Fz: %lf\n", Pot, Fx, Fy, Fz);
    fflush(stdout);
  }
Heinrich Schindler's avatar
Heinrich Schindler committed
581
582

  // constants of integration
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
  // The only logic for the Fy constant seems to be the fact that the
  // potential has a negative of this constant
  if (((X > (xlo + MINDIST)) && (X < (xhi - MINDIST))) &&
      ((Z > (zlo + MINDIST)) && (Z < (zhi - MINDIST)))) {
    Pot -= 2.0 * modY * ST_PI;
    if (SY != 0)
      Fy += 2.0 * (double)SY * ST_PI;
    else
      Fy = 2.0 * ST_PI;
  }
  if (DebugISLES) {
    printf("Constants of integration added for potential and Fy.\n");
    printf("Pot: %lf, Fx: %lf, Fy: %lf, Fz: %lf\n", Pot, Fx, Fy, Fz);
    fflush(stdout);
  }
Heinrich Schindler's avatar
Heinrich Schindler committed
598
599

  // Error situations handled before returning the values
600
  if ((Pot < 0.0) || (isnan(Pot) || isinf(Pot))) {
Heinrich Schindler's avatar
Heinrich Schindler committed
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
    fprintf(fIsles, "\n--- Approximation in ExactRecSurf ---\n");
    fprintf(fIsles, "Negative, nan or inf potential ... approximating!\n");
    if (DebugISLES) {
      fprintf(stdout, "\n--- Approximation in ExactRecSurf ---\n");
      fprintf(stdout, "Negative, nan or inf potential ... approximating!\n");
    }
    fprintf(fIsles, "X: %.16lg, Y: %.16lg, Z: %.16lg\n", X, Y, Z);
    fprintf(fIsles, "xlo: %.16lg, zlo: %.16lg, xhi: %.16lg, zhi: %.16lg\n", xlo,
            zlo, xhi, zhi);
    fprintf(fIsles, "dxlo: %.16lg, dzlo: %.16lg, dxhi: %.16lg, dzhi: %.16lg\n",
            dxlo, dzlo, dxhi, dzhi);
    fprintf(fIsles, "D11: %.16lg, D12: %.16lg, D21: %.16lg, D22: %.16lg\n", D11,
            D12, D21, D22);
    fprintf(fIsles, "modY: %.16lg\n", modY);
    fprintf(fIsles, "I1: %.16lg, I2: %.16lg, R1: %.16lg, R2: %.16lg\n", I1, I2,
            R1, R2);
    fprintf(fIsles, "S1: %d, S2: %d\n", S1, S2);
    fprintf(fIsles, "Computed Pot: %.16lg\n", Pot);
    ApproxFlag = 13;
    ++FailureCntr;
    --ExactCntr;
    fprintf(fIsles, "Approximating: %d.\n", ApproxFlag);
    return (ApproxRecSurf(X, Y, Z, xlo, zlo, xhi, zhi, XNSegApprox, ZNSegApprox,
                          Potential, Flux));
  }
626
  if (isnan(Fx) || isinf(Fx)) {
Heinrich Schindler's avatar
Heinrich Schindler committed
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
    fprintf(fIsles, "\n--- Approximation in ExactRecSurf ---\n");
    fprintf(fIsles, "Nan or inf Fx ... approximating!\n");
    if (DebugISLES) {
      fprintf(stdout, "\n--- Approximation in ExactRecSurf ---\n");
      fprintf(stdout, "Nan or inf Fx ... approximating!\n");
    }
    fprintf(fIsles, "X: %.16lg, Y: %.16lg, Z: %.16lg\n", X, Y, Z);
    fprintf(fIsles, "xlo: %.16lg, zlo: %.16lg, xhi: %.16lg, zhi: %.16lg\n", xlo,
            zlo, xhi, zhi);
    fprintf(fIsles, "dxlo: %.16lg, dzlo: %.16lg, dxhi: %.16lg, dzhi: %.16lg\n",
            dxlo, dzlo, dxhi, dzhi);
    fprintf(fIsles, "D11: %.16lg, D12: %.16lg, D21: %.16lg, D22: %.16lg\n", D11,
            D12, D21, D22);
    fprintf(fIsles, "Computed Fx: %.16lg\n", Fx);
    ApproxFlag = 14;
    ++FailureCntr;
    --ExactCntr;
    fprintf(fIsles, "Approximating: %d.\n", ApproxFlag);
    return (ApproxRecSurf(X, Y, Z, xlo, zlo, xhi, zhi, XNSegApprox, ZNSegApprox,
                          Potential, Flux));
  }
648
  if (isnan(Fy) || isinf(Fy)) {
Heinrich Schindler's avatar
Heinrich Schindler committed
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
    fprintf(fIsles, "\n--- Approximation in ExactRecSurf ---\n");
    fprintf(fIsles, "Nan or inf Fy ... approximating!\n");
    if (DebugISLES) {
      fprintf(stdout, "\n--- Approximation in ExactRecSurf ---\n");
      fprintf(stdout, "Nan or inf Fy ... approximating!\n");
    }
    fprintf(fIsles, "X: %.16lg, Y: %.16lg, Z: %.16lg\n", X, Y, Z);
    fprintf(fIsles, "xlo: %.16lg, zlo: %.16lg, xhi: %.16lg, zhi: %.16lg\n", xlo,
            zlo, xhi, zhi);
    fprintf(fIsles, "dxlo: %.16lg, dzlo: %.16lg, dxhi: %.16lg, dzhi: %.16lg\n",
            dxlo, dzlo, dxhi, dzhi);
    fprintf(fIsles, "D11: %.16lg, D12: %.16lg, D21: %.16lg, D22: %.16lg\n", D11,
            D12, D21, D22);
    fprintf(fIsles, "S1: %d, S2: %d, SY: %d, modY: %.16lg\n", S1, S2, SY, modY);
    fprintf(fIsles, "I1: %.16lg, I2: %.16lg, R1: %.16lg, R2: %.16lg\n", I1, I2,
            R1, R2);
    fprintf(fIsles, "Computed Fy: %.16lg\n", Fy);
    ApproxFlag = 15;
    ++FailureCntr;
    --ExactCntr;
    fprintf(fIsles, "Approximating: %d.\n", ApproxFlag);
    return (ApproxRecSurf(X, Y, Z, xlo, zlo, xhi, zhi, XNSegApprox, ZNSegApprox,
                          Potential, Flux));
  }
673
  if (isnan(Fz) || isinf(Fz)) {
Heinrich Schindler's avatar
Heinrich Schindler committed
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
    fprintf(fIsles, "\n--- Approximation in ExactRecSurf ---\n");
    fprintf(fIsles, "Nan or inf Fz ... approximating!\n");
    if (DebugISLES) {
      fprintf(stdout, "\n--- Approximation in ExactRecSurf ---\n");
      fprintf(stdout, "Nan or inf Fz ... approximating!\n");
    }
    fprintf(fIsles, "X: %.16lg, Y: %.16lg, Z: %.16lg\n", X, Y, Z);
    fprintf(fIsles, "xlo: %.16lg, zlo: %.16lg, xhi: %.16lg, zhi: %.16lg\n", xlo,
            zlo, xhi, zhi);
    fprintf(fIsles, "dxlo: %.16lg, dzlo: %.16lg, dxhi: %.16lg, dzhi: %.16lg\n",
            dxlo, dzlo, dxhi, dzhi);
    fprintf(fIsles, "D11: %.16lg, D12: %.16lg, D21: %.16lg, D22: %.16lg\n", D11,
            D12, D21, D22);
    fprintf(fIsles, "Computed Fz: %.16lg\n", Fz);
    ApproxFlag = 16;
    ++FailureCntr;
    --ExactCntr;
    fprintf(fIsles, "Approximating: %d.\n", ApproxFlag);
    return (ApproxRecSurf(X, Y, Z, xlo, zlo, xhi, zhi, XNSegApprox, ZNSegApprox,
                          Potential, Flux));
  }

  *Potential = Pot;
  Flux->X = Fx;
  Flux->Y = Fy;
  Flux->Z = Fz;
700

Heinrich Schindler's avatar
Heinrich Schindler committed
701
  if (DebugISLES) printf("Going out of ExactRecSurf ...\n");
702

Heinrich Schindler's avatar
Heinrich Schindler committed
703
704
  return 0;
}  // ExactRecSurf ends
705

Heinrich Schindler's avatar
Heinrich Schindler committed
706
707
708
709
710
711
int ApproxRecSurf(double X, double Y, double Z, double xlo, double zlo,
                  double xhi, double zhi, int xseg, int zseg, double *Potential,
                  Vector3D *Flux) {
  if (DebugISLES) {
    printf("In ApproxRecSurf ...\n");
  }
712

Heinrich Schindler's avatar
Heinrich Schindler committed
713
  ++ApproxCntr;
714

715
716
717
718
  double dx = (xhi - xlo) / xseg;
  double dz = (zhi - zlo) / zseg;
  double xel = (xhi - xlo) / xseg;
  double zel = (zhi - zlo) / zseg;
Heinrich Schindler's avatar
Heinrich Schindler committed
719
720
  double diag = sqrt(dx * dx + dz * dz);
  double area = xel * zel;
721

722
  double Pot = 0., XFlux = 0., YFlux = 0., ZFlux = 0.;
Heinrich Schindler's avatar
Heinrich Schindler committed
723

724
  if (area > MINDIST2) {  // else not necessary
725
726
727
728
    for (int i = 1; i <= xseg; ++i) {
      double xi = xlo + (dx / 2.0) + (i - 1) * dx;
      for (int k = 1; k <= zseg; ++k) {
        double zk = zlo + (dz / 2.0) + (k - 1) * dz;
Heinrich Schindler's avatar
Heinrich Schindler committed
729

730
        double dist = sqrt((X - xi) * (X - xi) + Y * Y + (Z - zk) * (Z - zk));
Heinrich Schindler's avatar
Heinrich Schindler committed
731
732
733
        if (DebugISLES) printf("dist: %lg\n", dist);
        if (dist >= diag) {
          Pot += area / dist;
734
735
        } else if (dist <= MINDIST) {
          // Self influence
Heinrich Schindler's avatar
Heinrich Schindler committed
736
737
          Pot += 2.0 * (xel * log((zel + sqrt(xel * xel + zel * zel)) / xel) +
                        zel * log((xel + sqrt(xel * xel + zel * zel)) / zel));
738
739
        } else {
          // in the intermediate region where diag > dist > MINDIST
Heinrich Schindler's avatar
Heinrich Schindler committed
740
741
          Pot += area / diag;  // replace by expression of self-influence
          if (DebugISLES) printf("Special Pot: %lg\n", area / diag);
742
743
        }

Heinrich Schindler's avatar
Heinrich Schindler committed
744
        if (dist >= diag) {
745
746
747
748
749
750
751
752
753
          double f = area / (dist * dist * dist);
          XFlux += f * (X - xi);
          YFlux += f * Y;
          ZFlux += f * (Z - zk);
        } else {
          double f = area / (diag * diag * diag);
          XFlux += f * (X - xi);
          YFlux += f * Y;
          ZFlux += f * (Z - zk);
Heinrich Schindler's avatar
Heinrich Schindler committed
754
          if (DebugISLES) {
755
756
            printf("Special XFlux: %lg, YFlux: %lg, ZFlux: %lg\n", f * (X - xi),
                   f * Y, f * (Z - zk));
757
          }
Heinrich Schindler's avatar
Heinrich Schindler committed
758
759
760
761
        }  // else dist >= diag
      }    // zseg
    }      // xseg
  }        // if area > MINDIST2
762

Heinrich Schindler's avatar
Heinrich Schindler committed
763
764
765
766
  *Potential = Pot;
  Flux->X = XFlux;
  Flux->Y = YFlux;
  Flux->Z = ZFlux;
767

Heinrich Schindler's avatar
Heinrich Schindler committed
768
769
  return 0;
}  // ApproxRecSurf ends
770
771
772
773
774
775
776
777
778
779
780
781
782

// Exact potential for a triangular element (a right-angled
// triangle having (0,0), (1,0), (0,zMax) as its vertices).
// It is assumed that the affine transformation has been carried out separately
// and necessary adjustments to X,Y,Z have been carried out and supplied to
// the following function.
// Potential needs all the calculations, whereas, the fluxes need part of the
// information generated during the computation for potential. This is the
// reason why we have included flux computation along with the potential
// computation.
// Parameters like dxlo (X), dxhi (X-1.0), dzlo (Z), dzhi (Z-zMax) can be used
// as done in ExactRectSurf
// Expressions from:
Heinrich Schindler's avatar
Heinrich Schindler committed
783
784
int ExactTriSurf(double zMax, double X, double Y, double Z, double *Potential,
                 Vector3D *Flux) {
Heinrich Schindler's avatar
Heinrich Schindler committed
785
  if (DebugISLES) printf("In ExactTriSurf ...\n");
786

Heinrich Schindler's avatar
Heinrich Schindler committed
787
788
  ++IslesCntr;
  ++ExactCntr;
Heinrich Schindler's avatar
Heinrich Schindler committed
789
  // Reset the flag to indicate approximate evaluation of potential.
790
  ApproxFlag = 0;
Heinrich Schindler's avatar
Heinrich Schindler committed
791

792
  // We do not need to check for X, since element extent is always 0 - 1.
Heinrich Schindler's avatar
Heinrich Schindler committed
793
794
795
  if (zMax < 3.0 * SHIFT * MINDIST) {
    // should allow enough space for Z corrections
    // One SHIFT should not lead to another criticality
Heinrich Schindler's avatar
Heinrich Schindler committed
796
797
    fprintf(stdout, "Element size too small! ... returning ...\n");
    return -1;
798
  }
Heinrich Schindler's avatar
Heinrich Schindler committed
799

800
  if (zMax > ARMAX || zMax < ARMIN) {
Heinrich Schindler's avatar
Heinrich Schindler committed
801
802
    fprintf(stdout, "Element too thin! ... returning ...\n");
    return -1;
803
804
  }

Heinrich Schindler's avatar
Heinrich Schindler committed
805
806
807
808
809
810
  // These three parameters can never be zero except at the three corners where
  // one of them can become zero. For example, at X=0, Y=0, Z=0, D11
  // is zero but the others are nonzero.
  if (fabs(X) < MINDIST) X = 0.0;
  if (fabs(Y) < MINDIST) Y = 0.0;
  if (fabs(Z) < MINDIST) Z = 0.0;
811
  double modY = fabs(Y);
Heinrich Schindler's avatar
Heinrich Schindler committed
812
  if (modY < MINDIST) modY = 0.0;
813
814
  int S1 = Sign(Z);
  int SY = Sign(Y);
Heinrich Schindler's avatar
Heinrich Schindler committed
815
816

  // distances from corners (0,0,0), (1,0,0) and (0,0,zMax)
817
  double D11 = sqrt(X * X + Y * Y + Z * Z);
Heinrich Schindler's avatar
Heinrich Schindler committed
818
  if (D11 < MINDIST) D11 = 0.0;
819
  double D21 = sqrt((X - 1.0) * (X - 1.0) + Y * Y + Z * Z);
Heinrich Schindler's avatar
Heinrich Schindler committed
820
  if (D21 < MINDIST) D21 = 0.0;
821
  double D12 = sqrt(X * X + Y * Y + (Z - zMax) * (Z - zMax));
Heinrich Schindler's avatar
Heinrich Schindler committed
822
823
  if (D12 < MINDIST) D12 = 0.0;

824
  double G = zMax * (X - 1.0) + Z;
Heinrich Schindler's avatar
Heinrich Schindler committed
825
  if (fabs(G) < MINDIST) G = 0.0;
826
  double E1 = (X - zMax * (Z - zMax));
Heinrich Schindler's avatar
Heinrich Schindler committed
827
  if (fabs(E1) < MINDIST) E1 = 0.0;
828
  double E2 = (X - 1.0 - zMax * Z);
Heinrich Schindler's avatar
Heinrich Schindler committed
829
  if (fabs(E2) < MINDIST) E2 = 0.0;
830
  double H1 = Y * Y + (Z - zMax) * G;
Heinrich Schindler's avatar
Heinrich Schindler committed
831
  if (fabs(H1) < MINDIST2) H1 = 0.0;
832
  double H2 = Y * Y + Z * G;
Heinrich Schindler's avatar
Heinrich Schindler committed
833
  if (fabs(H2) < MINDIST2) H2 = 0.0;
834
  double R1 = Y * Y + Z * Z;
Heinrich Schindler's avatar
Heinrich Schindler committed
835
  if (fabs(R1) < MINDIST2) R1 = 0.0;
836
  double I1 = modY * X;
Heinrich Schindler's avatar
Heinrich Schindler committed
837
  if (fabs(I1) < MINDIST2) I1 = 0.0;
838
  double I2 = modY * (X - 1.0);
Heinrich Schindler's avatar
Heinrich Schindler committed
839
  if (fabs(I2) < MINDIST2) I2 = 0.0;
840
  double Hypot = sqrt(1.0 + zMax * zMax);
841
  if (Hypot < MINDIST) {  // superfluous
Heinrich Schindler's avatar
Heinrich Schindler committed
842
843
844
    fprintf(stdout, "Hypotenuse too small! ... returning ...\n");
    return -1;
  }
845
  double Zhyp = zMax * (1.0 - X);  // Z on hypotenuse (extended) for given X
Heinrich Schindler's avatar
Heinrich Schindler committed
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867

  if (DebugISLES) {
    printf("\n\nzMax: %.16lg, X: %.16lg, Y: %.16lg, Z: %.16lg\n", zMax, X, Y,
           Z);
    printf("D11: %.16lg, D21: %.16lg, D12: %.16lg, Hypot: %.16lg\n", D11, D21,
           D12, Hypot);
    printf("modY: %.16lg, G: %.16lg\n", modY, G);
    printf("E1: %.16lg, E2: %.16lg, H1: %.16lg, H2: %.16lg\n", E1, E2, H1, H2);
    printf("S1: %d, SY: %d, R1: %.16lg, I1: %.16lg, I2: %.16lg\n", S1, SY, R1,
           I1, I2);
    printf("H1+G*D12: %.16lg, E1-zMax*D12: %.16lg\n", H1 + G * D12,
           E1 - zMax * D12);
    printf("H2+G*D21: %.16lg, E2-zMax*D21: %.16lg\n", H2 + G * D21,
           E2 - zMax * D21);
    printf("D11*fabs(Z): %.16lg, D21*fabs(Z): %.16lg\n", D11 * fabs(Z),
           D21 * fabs(Z));
    printf("Hypot*D12 - E1: %.16lg\n", Hypot * D12 - E1);
    printf("Hypot*D21 - E2: %.16lg\n\n", Hypot * D21 - E2);
    fprintf(stdout, "MINDIST: %.16lg, MINDIST2: %.16lg, SHIFT: %.16lg\n",
            MINDIST, MINDIST2, SHIFT);
    fflush(stdout);
  }
868

869
  // Check for possible numerical difficulties and take care.
Heinrich Schindler's avatar
Heinrich Schindler committed
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
  // Presently the idea is to shift the field point slightly to a 'safe'
  // position. Note that a shift in Y does not work because the singularities
  // are associated with D11-s and dxlo-s and their combinations. A Y shift can
  // sometimes alleviate the problem, but it cannot gurantee a permanent1s
  // solution.
  // A better approach is to evaluate analytic expressions truly valid for
  // these critical regions, especially to take care of the >= and <=
  // comparisons. For the == case, both of the previous comparisons are true and
  // it is hard to justify one choice over the other. On the other hand, there
  // are closed form analytic solutions for the == cases, and the problem does
  // not stem from round-off errors as in the > or <, but not quite == cases.
  // Possible singularity at D21 corner where X=1, Z=0
  // Possible singularity at D12 corner where X=0, Z=zMax
  // Check for possible numerical difficuties and take care
  // Note that modY = 0 cannot be a condition of criticality since considering
  // that would mean omitting even the barycenter properties.
  // Also note that shifts that are exactly equal to MINDIST, can give rise to
  // un-ending recursions, leading to Segmentation fault.
  // Special points and combinations

  // Three corners
  if ((fabs(D11) <= MINDIST)) {
    if (DebugISLES) printf("D11 <= MINDIST\n");
893
894
895
896
    double X1 = X;
    double Z1 = Z;
    if ((X >= 0.0) && (Z >= 0.0)) {
      // point on the element
Heinrich Schindler's avatar
Heinrich Schindler committed
897
898
899
      if (DebugISLES) printf("Case1=>\n");
      X1 = X + SHIFT * MINDIST;
      Z1 = Z + SHIFT * MINDIST;
900
901
    } else if ((X <= 0.0) && (Z >= 0.0)) {
      // field point outside the element
Heinrich Schindler's avatar
Heinrich Schindler committed
902
903
904
      if (DebugISLES) printf("Case2=>\n");
      X1 = X - SHIFT * MINDIST;
      Z1 = Z + SHIFT * MINDIST;
905
906
    } else if ((X >= 0.0) && (Z <= 0.0)) {
      // field point outside the element
Heinrich Schindler's avatar
Heinrich Schindler committed
907
908
909
      if (DebugISLES) printf("Case3=>\n");
      X1 = X + SHIFT * MINDIST;
      Z1 = Z - SHIFT * MINDIST;
910
911
    } else if ((X <= 0.0) && (Z <= 0.0)) {
      // field point outside the element
Heinrich Schindler's avatar
Heinrich Schindler committed
912
913
914
915
      if (DebugISLES) printf("Case4=>\n");
      X1 = X - SHIFT * MINDIST;
      Z1 = Z - SHIFT * MINDIST;
    }
916
917
918
919
920
921
922
923
    double Pot1;
    Vector3D Flux1;
    ExactTriSurf(zMax, X1, Y, Z1, &Pot1, &Flux1);
    *Potential = Pot1;
    Flux->X = Flux1.X;
    Flux->Y = Flux1.Y;
    Flux->Z = Flux1.Z;
    return 0;
Heinrich Schindler's avatar
Heinrich Schindler committed
924
925
926
  }
  if ((fabs(D21) <= MINDIST)) {
    if (DebugISLES) printf("D21 <= MINDIST\n");
927
928
929
930
    double X1 = X;
    double Z1 = Z;
    if ((X >= 1.0) && (Z >= 0.0)) {
      // point outside the element
Heinrich Schindler's avatar
Heinrich Schindler committed
931
932
933
      if (DebugISLES) printf("Case1=>\n");
      X1 = X + SHIFT * MINDIST;
      Z1 = Z + SHIFT * MINDIST;
934
935
936
    } else if ((X <= 1.0) && (Z >= 0.0)) {
      // field point on the element
      // difficult to decide (chk figure) - chk whether Z is beyond Zhyp
Heinrich Schindler's avatar
Heinrich Schindler committed
937
938
939
      if (DebugISLES) printf("Case2=>\n");
      X1 = X - SHIFT * MINDIST;
      Z1 = Z + SHIFT * MINDIST;
940
941
    } else if ((X >= 1.0) && (Z <= 0.0)) {
      // field point outside the element
Heinrich Schindler's avatar
Heinrich Schindler committed
942
943
944
      if (DebugISLES) printf("Case3=>\n");
      X1 = X + SHIFT * MINDIST;
      Z1 = Z - SHIFT * MINDIST;
945
946
    } else if ((X <= 1.0) && (Z <= 0.0)) {
      // field point outside the element
Heinrich Schindler's avatar
Heinrich Schindler committed
947
948
949
950
      if (DebugISLES) printf("Case4=>\n");
      X1 = X - SHIFT * MINDIST;
      Z1 = Z - SHIFT * MINDIST;
    }
951
952
953
954
955
956
957
958
    double Pot1;
    Vector3D Flux1;
    ExactTriSurf(zMax, X1, Y, Z1, &Pot1, &Flux1);
    *Potential = Pot1;
    Flux->X = Flux1.X;
    Flux->Y = Flux1.Y;
    Flux->Z = Flux1.Z;
    return 0;
Heinrich Schindler's avatar
Heinrich Schindler committed
959
960
961
  }
  if ((fabs(D12) <= MINDIST)) {
    if (DebugISLES) printf("D12 <= MINDIST\n");
962
963
964
965
    double X1 = X;
    double Z1 = Z;
    if ((X >= 0.0) && (Z >= zMax)) {
      // point outside the element
Heinrich Schindler's avatar
Heinrich Schindler committed
966
967
968
      if (DebugISLES) printf("Case1=>\n");
      X1 = X + SHIFT * MINDIST;
      Z1 = Z + SHIFT * MINDIST;
969
970
    } else if ((X <= 0.0) && (Z >= zMax)) {
      // field point outside the element
Heinrich Schindler's avatar
Heinrich Schindler committed
971
972
973
      if (DebugISLES) printf("Case2=>\n");
      X1 = X - SHIFT * MINDIST;
      Z1 = Z + SHIFT * MINDIST;
974
975
976
    } else if ((X >= 0.0) && (Z <= zMax)) {
      // field point on the element
      // can create problem for small element - chk whether Z is beyond Zhyp
Heinrich Schindler's avatar
Heinrich Schindler committed
977
978
979
      if (DebugISLES) printf("Case3=>\n");
      X1 = X + SHIFT * MINDIST;
      Z1 = Z - SHIFT * MINDIST;
980
981
    } else if ((X <= 0.0) && (Z <= zMax)) {
      // field point outside the element
Heinrich Schindler's avatar
Heinrich Schindler committed
982
983
984
985
      if (DebugISLES) printf("Case4=>\n");
      X1 = X - SHIFT * MINDIST;
      Z1 = Z - SHIFT * MINDIST;
    }
986
987
988
989
990
991
992
993
    double Pot1;
    Vector3D Flux1;
    ExactTriSurf(zMax, X1, Y, Z1, &Pot1, &Flux1);
    *Potential = Pot1;
    Flux->X = Flux1.X;
    Flux->Y = Flux1.Y;
    Flux->Z = Flux1.Z;
    return 0;
Heinrich Schindler's avatar
Heinrich Schindler committed
994
995
  }
  // Three Y lines at corners
996
997
  if ((fabs(X) < MINDIST) && (fabs(Z) < MINDIST)) {
    // Y line at (0,0,0) corner
Heinrich Schindler's avatar
Heinrich Schindler committed
998
    if (DebugISLES) printf("Y line at (0,0,0) corner\n");
999
1000
1001
1002
    double X1 = X;
    double Z1 = Z;
    if ((X >= 0.0) && (Z >= 0.0)) {
      // point on the element
Heinrich Schindler's avatar
Heinrich Schindler committed
1003
1004
1005
      if (DebugISLES) printf("Case1=>\n");
      X1 = X + SHIFT * MINDIST;
      Z1 = Z + SHIFT * MINDIST;
1006
1007
    } else if ((X <= 0.0) && (Z >= 0.0)) {
      // field point outside the element
Heinrich Schindler's avatar
Heinrich Schindler committed
1008
1009
1010
      if (DebugISLES) printf("Case2=>\n");
      X1 = X - SHIFT * MINDIST;
      Z1 = Z + SHIFT * MINDIST;
1011
1012
    } else if ((X >= 0.0) && (Z <= 0.0)) {
      // field point outside the element
Heinrich Schindler's avatar
Heinrich Schindler committed
1013
1014
1015
      if (DebugISLES) printf("Case3=>\n");
      X1 = X + SHIFT * MINDIST;
      Z1 = Z - SHIFT * MINDIST;
1016
1017
    } else if ((X <= 0.0) && (Z <= 0.0)) {
      // field point outside the element
Heinrich Schindler's avatar
Heinrich Schindler committed
1018
1019
1020
1021
      if (DebugISLES) printf("Case4=>\n");
      X1 = X - SHIFT * MINDIST;
      Z1 = Z - SHIFT * MINDIST;
    }
1022
1023
1024
1025
1026
1027
1028
1029
    double Pot1;
    Vector3D Flux1;
    ExactTriSurf(zMax, X1, Y, Z1, &Pot1, &Flux1);
    *Potential = Pot1;
    Flux->X = Flux1.X;
    Flux->Y = Flux1.Y;
    Flux->Z = Flux1.Z;
    return 0;
Heinrich Schindler's avatar
Heinrich Schindler committed
1030
  }
1031
1032
  if ((fabs(X - 1.0) < MINDIST) && (fabs(Z) < MINDIST)) {
    // Y line at (1,0) corner
Heinrich Schindler's avatar
Heinrich Schindler committed
1033
    if (DebugISLES) printf("Y line at (1,0,0) corner\n");
1034
1035
1036
1037
    double X1 = X;
    double Z1 = Z;
    if ((X >= 1.0) && (Z >= 0.0)) {
      // point on the element
Heinrich Schindler's avatar
Heinrich Schindler committed
1038
1039
1040
      if (DebugISLES) printf("Case1=>\n");
      X1 = X + SHIFT * MINDIST;
      Z1 = Z + SHIFT * MINDIST;
1041
1042
    } else if ((X <= 1.0) && (Z >= 0.0)) {
      // field point outside the element
Heinrich Schindler's avatar
Heinrich Schindler committed
1043
1044
1045
      if (DebugISLES) printf("Case2=>\n");
      X1 = X - SHIFT * MINDIST;
      Z1 = Z + SHIFT * MINDIST;
1046
1047
    } else if ((X >= 1.0) && (Z <= 0.0)) {
      // field point outside the element
Heinrich Schindler's avatar
Heinrich Schindler committed
1048
1049
1050
      if (DebugISLES) printf("Case3=>\n");
      X1 = X + SHIFT * MINDIST;
      Z1 = Z - SHIFT * MINDIST;
1051
1052
    } else if ((X <= 1.0) && (Z <= 0.0)) {
      // field point outside the element
Heinrich Schindler's avatar
Heinrich Schindler committed
1053
1054
1055
1056
      if (DebugISLES) printf("Case4=>\n");
      X1 = X - SHIFT * MINDIST;
      Z1 = Z - SHIFT * MINDIST;
    }
1057
1058
1059
1060
1061
1062
1063
1064
    double Pot1;
    Vector3D Flux1;
    ExactTriSurf(zMax, X1, Y, Z1, &Pot1, &Flux1);
    *Potential = Pot1;
    Flux->X = Flux1.X;
    Flux->Y = Flux1.Y;
    Flux->Z = Flux1.Z;
    return 0;
Heinrich Schindler's avatar
Heinrich Schindler committed
1065
  }
1066
1067
  if ((fabs(X) < MINDIST) && (fabs(Z - zMax) < MINDIST)) {
    // Y line at (0,zMax)
Heinrich Schindler's avatar
Heinrich Schindler committed
1068
    if (DebugISLES) printf("Y line at (0,0,zMax) corner\n");
1069