1 /*
2 Copyright (c) 1990 Massachusetts Institute of Technology, Cambridge, MA.
3 All rights reserved.
4
5 This Agreement gives you, the LICENSEE, certain rights and obligations.
6 By using the software, you indicate that you have read, understood, and
7 will comply with the terms.
8
9 Permission to use, copy and modify for internal, noncommercial purposes
10 is hereby granted. Any distribution of this program or any part thereof
11 is strictly prohibited without prior written consent of M.I.T.
12
13 Title to copyright to this software and to any associated documentation
14 shall at all times remain with M.I.T. and LICENSEE agrees to preserve
15 same. LICENSEE agrees not to make any copies except for LICENSEE'S
16 internal noncommercial use, or to use separately any portion of this
17 software without prior written consent of M.I.T. LICENSEE agrees to
18 place the appropriate copyright notice on any such copies.
19
20 Nothing in this Agreement shall be construed as conferring rights to use
21 in advertising, publicity or otherwise any trademark or the name of
22 "Massachusetts Institute of Technology" or "M.I.T."
23
24 M.I.T. MAKES NO REPRESENTATIONS OR WARRANTIES, EXPRESS OR IMPLIED. By
25 way of example, but not limitation, M.I.T. MAKES NO REPRESENTATIONS OR
26 WARRANTIES OF MERCHANTABILITY OR FITNESS FOR ANY PARTICULAR PURPOSE OR
27 THAT THE USE OF THE LICENSED SOFTWARE COMPONENTS OR DOCUMENTATION WILL
28 NOT INFRINGE ANY PATENTS, COPYRIGHTS, TRADEMARKS OR OTHER RIGHTS.
29 M.I.T. shall not be held liable for any liability nor for any direct,
30 indirect or consequential damages with respect to any claim by LICENSEE
31 or any third party on account of or arising from this Agreement or use
32 of this software.
33 */
34
35 #include "mulGlobal.h"
36
37 #ifndef FALSE
38 #define FALSE 0
39 #endif
40 #ifndef TRUE
41 #define TRUE 1
42 #endif
43 #define XI 0
44 #define YI 1
45 #define ZI 2
46 #define EQUIV_TOL 1.0e-9
47 /*#define PLANAR_TOL 1.0e-3*/
48 #define PLANAR_TOL 1.0e-5
49 #define MAX_REAL 1.0e+20;
50
51 /* Obvious Constants. */
52 #define PI 3.1415927
53 #define TWOPI 6.2831853
54 #define HALFPI 1.5707963E+00
55
56 /* Constants to save typing. */
57 #define FIVE3 1.666666666667
58 #define SEVEN3 2.3333333333333
59 #define ONE6 0.16666666666667
60 #define ONE3 0.3333333333333
61 #define FT3 4.66666667
62
63
64 /* Defines breakpoints panel multipoles. */
65 #define LIMITFOURTH 9.0
66 #define LIMITSECOND 36.0
67
68 /* Constants used for Hart arctan approximation. */
69 #define B1 0.24091197
70 #define B2 3.7851122
71 #define B3 5.6770721
72 #define B4 5.6772854
73 #define B5 5.6770747
74
75 #define Dot_Product(V1,V2) (V1[XI]*V2[XI]+V1[YI]*V2[YI]+V1[ZI]*V2[ZI])
76 #define DotP_Product(V1,R,S,T) ((V1[XI])*(R)+(V1[YI])*(S)+(V1[ZI])*(T))
77
78
79 static int num2nd=0, num4th=0, numexact=0;
80 static int num2ndsav=0, num4thsav=0, numexactsav=0;
81
initcalcp(panel_list)82 initcalcp(panel_list)
83 charge *panel_list;
84 {
85 charge *pq, *npq;
86 double vtemp[3];
87 double length, maxlength, minlength, length20, length31, sum, sum2, delta;
88 double normalize();
89 int i, j, next;
90
91 #if JACDBG == ON
92 FILE *foo;
93
94 foo = fopen("corners.txt", "w");
95 #endif
96
97 for(i=0, pq = panel_list; pq != NULL; pq = pq->next) i++;
98 #if JACDBG == ON
99 printf("Initial number of panels = %d\n", i);
100 #endif
101
102 pq = panel_list;
103 while (pq != NULL) {
104 /* Calculate edge lengths. */
105 maxlength = 0.0;
106 minlength = MAX_REAL;
107 for(i=0; i < pq->shape; i++) {
108 if(i == (pq->shape -1)) next = 0;
109 else next = i + 1;
110 for(sum= 0, j = 0; j < 3; j++) {
111 delta = pq->corner[next][j] - pq->corner[i][j];
112 sum += delta * delta;
113 }
114 pq->length[i] = length = sqrt(sum);
115 maxlength = MAX(maxlength, length);
116 minlength = MIN(minlength, length);
117 }
118
119 /* Get diags and lengths. */
120 for(sum= 0, sum2 = 0, i = 0; i < 3; i++) {
121 pq->X[i] = delta = pq->corner[2][i] - pq->corner[0][i];
122 sum += delta * delta;
123 if(pq->shape == 3) pq->Y[i] = pq->corner[0][i] - pq->corner[1][i];
124 else {
125 pq->Y[i] = delta = pq->corner[3][i] - pq->corner[1][i];
126 sum2 += delta * delta;
127 }
128 }
129 length20 = sqrt(sum);
130 length31 = sqrt(sum2);
131
132 /* Check on lengths for quad. */
133 if(pq->shape == 3) {
134 pq->max_diag = maxlength;
135 pq->min_diag = minlength;
136 }
137 else {
138 length = MAX(length20, length31);
139 pq->max_diag = length;
140 pq->min_diag = MIN(length20, length31);
141 }
142
143 /* Z-axis is normal to two diags. */
144 Cross_Product(pq->X, pq->Y, pq->Z);
145 /*#if 1 == 0*/
146 if(flip_normal(pq)) {
147 for(i = 0; i < 3; i++) {
148 pq->Z[i] = -(pq->Z[i]); /* flip the normal */
149 pq->X[i] = -(pq->X[i]); /* flip the x direction */
150 /* interchange points 0 and 2 so that corner order will be
151 consistent with X flip (note this is OK for both quads and tris) */
152 vtemp[0] = pq->corner[0][i];
153 pq->corner[0][i] = pq->corner[2][i];
154 pq->corner[2][i] = vtemp[0];
155 }
156 /* interchange length entries in length vector */
157 vtemp[0] = pq->length[0];
158 pq->length[0] = pq->length[1];
159 pq->length[1] = vtemp[0];
160 if(pq->shape == 4) {
161 vtemp[0] = pq->length[2];
162 pq->length[2] = pq->length[3];
163 pq->length[3] = vtemp[0];
164 }
165 }
166 /*#endif*/
167 pq->area = 0.5 * normalize(pq->Z);
168 normalize(pq->X);
169
170 /* Real Y-axis is normal to X and Z (resulting system is left-handed). */
171 Cross_Product(pq->X, pq->Z, pq->Y);
172
173 /* Project the corner points into the plane defined by edge midpoints. */
174 if(planarize(pq) == FALSE) {
175 /* Planarization too drastic, crack into two triangles. */
176 CALLOC(npq, 1, charge, ON, AMSC);
177 npq->next = pq->next;
178 pq->next = npq;
179 npq->shape = 3;
180 pq->shape = 3;
181 npq->cond = pq->cond;
182 npq->surf = pq->surf;
183 VCOPY(npq->corner[0], pq->corner[0]);
184 VCOPY(npq->corner[1], pq->corner[2]);
185 VCOPY(npq->corner[2], pq->corner[3]);
186 continue;
187 }
188
189 /* Calculate the centroid. */
190 centroid(pq, length20);
191
192 /* Put corners in the newly defined coord system. */
193 for(i=0; i < pq->shape; i++) {
194 pq->corner[i][XI] -= pq->x;
195 pq->corner[i][YI] -= pq->y;
196 pq->corner[i][ZI] -= pq->z;
197 }
198 for(i=0; i < pq->shape; i++) {
199 vtemp[XI] = Dot_Product(pq->corner[i], pq->X);
200 vtemp[YI] = Dot_Product(pq->corner[i], pq->Y);
201 vtemp[ZI] = Dot_Product(pq->corner[i], pq->Z);
202 pq->corner[i][XI] = vtemp[XI];
203 pq->corner[i][YI] = vtemp[YI];
204 pq->corner[i][ZI] = vtemp[ZI];
205 if(fabs(pq->corner[i][ZI]) > (EQUIV_TOL * pq->min_diag)) {
206 printf("FATAL PROGRAM ERROR: renormalized z=%g\n", pq->corner[i][ZI]);
207 exit(1);
208 }
209 pq->corner[i][ZI] = 0.0;
210 }
211
212 #if JACDBG == ON /* dump corners, center to file */
213 fprintf(foo, "%g %g %g: ", pq->x, pq->y, pq->z);
214 for(i = 0; i < pq->shape; i++) {
215 fprintf(foo, "%g %g %g ", pq->corner[i][XI], pq->corner[i][YI],
216 pq->corner[i][ZI]);
217 }
218 fprintf(foo, "\n");
219 #endif
220
221 /* Finally, calculate the moments. */
222 ComputeMoments(pq);
223
224 /* Iterate for the next panel. */
225 pq = pq->next;
226
227 }
228
229 }
230
231 /*
232 determine if normal needs to be flipped to get dielectric bdry cond right
233 */
oldflip_normal(panel)234 oldflip_normal(panel)
235 charge *panel;
236 {
237 int i;
238 double x, y, z;
239 double ctr_minus_n[3], ctr_plus_n[3], norm_minus, norm_plus, norm, norm_sq;
240 surface *surf = panel->surf;
241 int ref_inside = surf->ref_inside, flip_normal;
242 double *ref = surf->ref, *normal;
243 char *surf_name = surf->name;
244
245 if(surf->type != DIELEC && surf->type != BOTH) return(FALSE);
246
247 /* get panel corner (relative to reference point) and normal */
248 x = panel->corner[0][0] - ref[0];
249 y = panel->corner[0][1] - ref[1];
250 z = panel->corner[0][2] - ref[2];
251 norm_sq = x*x + y*y + z*z;
252 norm = sqrt(norm_sq);
253 normal = panel->Z;
254
255 /* add the (scaled) normal and negative normal to the panel corner */
256 /* negative normal result should be closer to ref point if(ref_inside) */
257 ctr_minus_n[0] = x - 0.1*norm*normal[0];
258 ctr_minus_n[1] = y - 0.1*norm*normal[1];
259 ctr_minus_n[2] = z - 0.1*norm*normal[2];
260 ctr_plus_n[0] = x + 0.1*norm*normal[0];
261 ctr_plus_n[1] = y + 0.1*norm*normal[1];
262 ctr_plus_n[2] = z + 0.1*norm*normal[2];
263
264 /* get norms of test points, one inside (minus) other out (plus) */
265 norm_minus = ctr_minus_n[0]*ctr_minus_n[0];
266 norm_plus = ctr_plus_n[0]*ctr_plus_n[0];
267 for(i = 1; i < 3; i++) {
268 norm_minus += ctr_minus_n[i]*ctr_minus_n[i];
269 norm_plus += ctr_plus_n[i]*ctr_plus_n[i];
270 }
271
272 flip_normal = FALSE;
273 if(norm_minus > norm_sq) {
274 if(norm_plus > norm_sq) {
275 fprintf(stderr,
276 "flip_normal: both test points on non-reference side\n");
277 fprintf(stderr, " Surface: %s\n", hack_path(surf_name));
278 fprintf(stderr, " Translation: (%g %g %g)\n", surf->trans[0],
279 surf->trans[1], surf->trans[2]);
280 fprintf(stderr, " Reference point: (%g %g %g)\n",
281 ref[0], ref[1], ref[2]);
282 fprintf(stderr, " Panel corner: (%g %g %g)\n",
283 panel->corner[0][0], panel->corner[0][1], panel->corner[0][2]);
284 fprintf(stderr, " Normal: (%g %g %g)\n",
285 normal[0], normal[1], normal[2]);
286 exit(1);
287 }
288 if(ref_inside) flip_normal = TRUE;
289 }
290 else if(norm_plus < norm_sq) {
291 if(norm_minus < norm_sq) {
292 fprintf(stderr,
293 "flip_normal: both test points on reference point side\n");
294 fprintf(stderr, " Surface: %s\n", hack_path(surf_name));
295 fprintf(stderr, " Translation: (%g %g %g)\n", surf->trans[0],
296 surf->trans[1], surf->trans[2]);
297 fprintf(stderr, " Reference point: (%g %g %g)\n",
298 ref[0], ref[1], ref[2]);
299 fprintf(stderr, " Panel corner: (%g %g %g)\n",
300 panel->corner[0][0], panel->corner[0][1], panel->corner[0][2]);
301 fprintf(stderr, " Normal: (%g %g %g)\n",
302 normal[0], normal[1], normal[2]);
303 exit(1);
304 }
305 if(!ref_inside) flip_normal = TRUE;
306 }
307
308 return(flip_normal);
309 }
310
311 /*
312 determine if normal needs to be flipped to get dielectric bdry cond right
313 - this function uses 0.0 as a breakpoint when really machine precision
314 weighted checks should be done (really not an issue if ref point far)
315 */
flip_normal(panel)316 flip_normal(panel)
317 charge *panel;
318 {
319 int i;
320 double x, y, z;
321 double ctr_minus_n[3], ctr_plus_n[3], norm_minus, norm_plus, norm, norm_sq;
322 surface *surf = panel->surf;
323 int ref_inside = surf->ref_inside, flip_normal;
324 double *ref = surf->ref, *normal, angle, norm_n;
325 char *surf_name = surf->name;
326
327 if(surf->type != DIELEC && surf->type != BOTH) return(FALSE);
328
329 /* get panel corner (relative to reference point) and normal */
330 x = panel->corner[0][0] - ref[0];
331 y = panel->corner[0][1] - ref[1];
332 z = panel->corner[0][2] - ref[2];
333 norm_sq = x*x + y*y + z*z;
334 norm = sqrt(norm_sq);
335 normal = panel->Z;
336 norm_n =
337 sqrt(normal[0]*normal[0] + normal[1]*normal[1] + normal[2]*normal[2]);
338
339 /* figure the dot product between the normal and the corner-ref point
340 - ref_inside and angle <= 90 => flip
341 - ref_inside and angle > 90 => no flip
342 - !ref_inside and angle <= 90 => no flip
343 - !ref_inside and angle > 90 => flip */
344 /* figure angle */
345 angle = (x*normal[0] + y*normal[1] + z*normal[2])/(norm*norm_n);
346 if(ref_inside && angle <= 0.0) flip_normal = TRUE;
347 else if(ref_inside && angle > 0.0) flip_normal = FALSE;
348 else if(!ref_inside && angle <= 0.0) flip_normal = FALSE;
349 else if(!ref_inside && angle > 0.0) flip_normal = TRUE;
350 else {
351 fprintf(stderr,
352 "flip_normal: inconclusive test for normal flipping\n");
353 fprintf(stderr, " Surface: %s\n", hack_path(surf_name));
354 fprintf(stderr, " Translation: (%g %g %g)\n", surf->trans[0],
355 surf->trans[1], surf->trans[2]);
356 fprintf(stderr, " Reference point: (%g %g %g)\n",
357 ref[0], ref[1], ref[2]);
358 fprintf(stderr, " Panel corner: (%g %g %g)\n",
359 panel->corner[0][0], panel->corner[0][1], panel->corner[0][2]);
360 fprintf(stderr, " Normal: (%g %g %g)\n",
361 normal[0], normal[1], normal[2]);
362 exit(1);
363 }
364
365 return(flip_normal);
366 }
367
368
369 /*
370 Changes the corner points so that they lie in the plane defined by the
371 panel diagonals and any midpoint of an edge.
372 */
planarize(pq)373 planarize(pq)
374 charge *pq;
375 {
376 double origin[3], corner[3], delta[4][3], px, py, dx, dy, dz;
377 int i, j, numcorners = pq->shape;
378 double tolsq = PLANAR_TOL * pq->min_diag;
379
380 tolsq *= tolsq;
381
382 /* Triangular panels are just fine already. */
383 if(numcorners != 4) return(TRUE);
384
385 /* Pick edge midpoint as origin. */
386 for(i=0; i < 3; i++) origin[i] = 0.5 * (pq->corner[1][i] + pq->corner[0][i]);
387
388 for(i=0; i < numcorners; i++) {
389 for(j=0; j < 3; j++) corner[j] = pq->corner[i][j] - origin[j];
390 px = Dot_Product(corner, pq->X);
391 py = Dot_Product(corner, pq->Y);
392
393 dx = px * pq->X[XI] + py * pq->Y[XI] + origin[XI] - pq->corner[i][XI];
394 dy = px * pq->X[YI] + py * pq->Y[YI] + origin[YI] - pq->corner[i][YI];
395 dz = px * pq->X[ZI] + py * pq->Y[ZI] + origin[ZI] - pq->corner[i][ZI];
396
397 delta[i][XI] = dx;
398 delta[i][YI] = dy;
399 delta[i][ZI] = dz;
400
401 /* Moved too much, crack the panel. */
402 if((dx * dx + dy * dy + dz * dz) > tolsq) return(FALSE);
403 }
404
405 /* Still Here? Must be planarizing. */
406 for(i=0; i < numcorners; i++) {
407 for(j=0; j < 3; j++) {
408 pq->corner[i][j] += delta[i][j];
409 }
410 }
411 return(TRUE);
412 }
413
414
415 /*
416 Determines centroid of a panel (defined as the point which make the
417 first moments vanish. Calculation begins by projection into the
418 coordinate system defined by the panel normal as the z-axis and
419 edge02 as the x-axis.
420 */
centroid(pp,x2)421 centroid(pp, x2)
422 charge *pp;
423 double x2;
424 {
425 double vertex1[3], vertex3[3];
426 double sum, dl, x1, y1, x3, y3, xc, yc;
427 int i;
428
429 /* Use vertex 0 as the origin. */
430 for(i=0; i< 3; i++) {
431 vertex1[i] = pp->corner[1][i] - pp->corner[0][i];
432 if(pp->shape == 4) vertex3[i] = pp->corner[3][i] - pp->corner[0][i];
433 else vertex3[i] = pp->corner[2][i] - pp->corner[0][i];
434 }
435
436 /* Project into the panel axes. */
437 y1 = Dot_Product(vertex1, pp->Y);
438 y3 = Dot_Product(vertex3, pp->Y);
439 x1 = Dot_Product(vertex1, pp->X);
440 x3 = Dot_Product(vertex3, pp->X);
441
442 yc = ONE3 * (y1 + y3);
443 xc = ONE3 * (x2 + ((x1 * y1 - x3 * y3)/(y1 - y3)));
444
445 pp->x = pp->corner[0][XI] + xc * pp->X[XI] + yc * pp->Y[XI];
446 pp->y = pp->corner[0][YI] + xc * pp->X[YI] + yc * pp->Y[YI];
447 pp->z = pp->corner[0][ZI] + xc * pp->X[ZI] + yc * pp->Y[ZI];
448
449 }
450
normalize(vector)451 double normalize(vector)
452 double vector[3];
453 {
454 double length;
455 int i;
456
457 length = sqrt( vector[0]*vector[0]
458 + vector[1]*vector[1]
459 + vector[2]*vector[2]);
460
461 for (i=0; i<3; i++) vector[i] = vector[i] / length;
462
463 return length;
464 }
465
466 /* Assumes the vectors are normalized. */
If_Equal(vector1,vector2)467 int If_Equal(vector1, vector2)
468 double vector1[3], vector2[3];
469 {
470 int i;
471
472 for (i=0; i<3; i++)
473 if (fabs(vector1[i] - vector2[i]) > EQUIV_TOL) return FALSE;
474 return TRUE;
475 }
476
477 /* Calculates result_vector = vector1 X vector2. */
Cross_Product(vector1,vector2,result_vector)478 Cross_Product(vector1, vector2, result_vector)
479 double vector1[], vector2[], result_vector[];
480 {
481 result_vector[XI] = vector1[YI]*vector2[ZI] - vector1[ZI]*vector2[YI];
482 result_vector[YI] = vector1[ZI]*vector2[XI] - vector1[XI]*vector2[ZI];
483 result_vector[ZI] = vector1[XI]*vector2[YI] - vector1[YI]*vector2[XI];
484 }
485
486 /*
487 Computes the potential at x, y, z due to a unit source on panel
488 and due to a dipole is returned in pfd.
489
490 Note, the code is subtle because there are 5 cases depending on
491 the placement of the collocation point
492
493 CASE1: evaluation point projection strictly outside of panel (happens most
494 of the time).
495 CASE2: eval pnt proj. strictly inside panel (happens >= #panels times,
496 at least for the self terms)
497 CASE3: eval pnt proj. on a panel side, not on a corner (happens
498 when paneled faces meet at right angles, also possible other ways)
499 CASE4: eval pnt proj. on a panel corner (happens rarely to never)
500 CASE5: eval pnt proj. on side extension (happens when paneled
501 faces meet at right angles, also possible other ways).
502 */
calcp(panel,x,y,z,pfd)503 double calcp(panel, x, y, z, pfd)
504 charge *panel;
505 double x, y, z, *pfd;
506 {
507 double r[4], fe[4], xmxv[4], ymyv[4];
508 double xc, yc, zc, zsq, xn, yn, zn, znabs, xsq, ysq, rsq, diagsq, dtol;
509 double v, arg, st, ct, length, s1, c1, s2, c2, s12, c12, val;
510 double *s;
511 double rInv, r2Inv, r3Inv, r5Inv, r7Inv, r9Inv, zr2Inv;
512 double ss1, ss3, ss5, ss7, ss9;
513 double s914, s813, s411, s512, s1215;
514 double fs, fd, fdsum;
515 int okay, i, next;
516 struct edge *edge;
517 double *corner;
518
519 /* Put the evaluation point into this panel's coordinates. */
520 xc = x - panel->x;
521 yc = y - panel->y;
522 zc = z - panel->z;
523
524 xn = DotP_Product(panel->X, xc, yc, zc);
525 yn = DotP_Product(panel->Y, xc, yc, zc);
526 zn = DotP_Product(panel->Z, xc, yc, zc);
527
528 zsq = zn * zn;
529 xsq = xn * xn;
530 ysq = yn * yn;
531 rsq = zsq + xsq + ysq;
532 diagsq = panel->max_diag * panel->max_diag;
533
534 if(rsq > (LIMITFOURTH * diagsq)) {
535 fs = 0.0; fd = 0.0;
536 s = panel->moments;
537 /* First, second moments. */
538 r2Inv = 1.0 / rsq;
539 rInv = sqrt(r2Inv);
540 r3Inv = r2Inv * rInv;
541 r5Inv = r3Inv * r2Inv;
542 zr2Inv = zn * r2Inv;
543 ss1 = s[1] * rInv;
544 ss3 = -(s[3] + s[10]) * r3Inv;
545 ss5 = (xsq * s[10] + (xn * yn * s[7]) + ysq * s[3]) * r5Inv;
546 fs = ss1 + ONE3 * ss3 + ss5;
547 fdsum = ss1 + ss3 + 5.0 * ss5;
548 fd = zr2Inv * fdsum;
549 if(rsq < (LIMITSECOND * diagsq)) {
550 /* Third and fourth moments added for diagsq/r2 between 40 and 150. */
551 s914 = s[9] + s[14];
552 s813 = s[8] + s[13];
553 s411 = s[4] + s[11];
554 s512 = s[5] + s[12];
555 s1215 = s[12] + s[15];
556 r7Inv = r5Inv * r2Inv;
557 r9Inv = r7Inv * r2Inv;
558 ss5 = (-xn * s813 - yn * s411 + 0.1 * (s512 + s1215)) * r5Inv;
559
560 ss7 = (FIVE3 *((xn * xsq * s[13] + yn * ysq * s[4])
561 + 3.0 * xn * yn * (xn * s[11] + yn * s[8]))
562 - xsq * s1215 - ysq * s512 - xn * yn * s914) * r7Inv;
563
564 ss9 = (7.0 * (ONE6 * (xsq * xsq * s[15] + ysq * ysq * s[5])
565 + xsq * ysq * s[12])
566 + SEVEN3 * xn * yn * (xsq * s[14] + ysq * s[9])) * r9Inv;
567
568 fs += ss5 + ss7 + ss9;
569 fdsum = 5.0 * ss5 + 7.0 * ss7 + 9.0 * ss9;
570 fd += zr2Inv * fdsum;
571 num4th++;
572 }
573 else num2nd++;
574 }
575 else {
576 dtol = EQUIV_TOL * panel->min_diag;
577 znabs = fabs(zn);
578
579 /* Always move the evaluation point a little bit off the panel. */
580 if(znabs < dtol) {
581 zn = 0.5 * dtol; /* Half of dtol insures detection for zero dipole. */
582 znabs = 0.5 * dtol;
583 }
584
585 /* Once per corner computations. */
586 for(okay = TRUE, i=0; i < panel->shape; i++) {
587 corner = panel->corner[i];
588 xmxv[i] = xc = xn - corner[XI];
589 ymyv[i] = yc = yn - corner[YI];
590 zc = zn - corner[ZI];
591 fe[i] = xc * xc + zc * zc;
592 r[i] = sqrt(yc * yc + fe[i]);
593 if(r[i] < (1.005 * znabs)) { /* If r almost z, on vertex normal. */
594 okay = FALSE;
595 }
596 }
597
598 /* Once per edge computations. */
599 fs = 0.0; fd = 0.0;
600 for(i=0; i < panel->shape; i++) {
601 if(i == (panel->shape - 1)) next = 0;
602 else next = i + 1;
603
604 /* Now calculate the edge contributions to a panel. */
605 length = panel->length[i];
606 ct = (panel->corner[next][XI] - panel->corner[i][XI]) / length;
607 st = (panel->corner[next][YI] - panel->corner[i][YI]) / length;
608
609 /* v is projection of eval-i edge onto perpend to next-i edge. */
610 /* Exploits the fact that corner points in panel coordinates. */
611 v = xmxv[i] * st - ymyv[i] * ct;
612
613 /* arg == zero if eval on next-i edge, but then v = 0. */
614 arg = (r[i] + r[next] - length)/(r[i] + r[next] + length);
615 if(arg > 0.0) fs -= v * log(arg);
616
617 /* Okay means eval not near a vertex normal, Use Hess-Smith. */
618 if(okay) {
619 s1 = v * r[i];
620 c1 = znabs * (xmxv[i] * ct + ymyv[i] * st);
621 s2 = v * r[next];
622 c2 = znabs * (xmxv[next] * ct + ymyv[next] * st);
623 }
624 /* Near a vertex normal, use Newman. */
625 else {
626 s1 = (fe[i] * st) - (xmxv[i] * ymyv[i] * ct);
627 c1 = znabs * r[i] * ct;
628 s2 = (fe[next] * st) - (xmxv[next] * ymyv[next] * ct);
629 c2 = znabs * r[next] * ct;
630 }
631
632 s12 = (s1 * c2) - (s2 * c1);
633 c12 = (c1 * c2) + (s1 * s2);
634 val = atan2(s12, c12);
635 fd += val;
636 }
637 /* Adjust the computed values. */
638
639 if(fd < 0.0) fd += TWOPI;
640 if(zn < 0.0) fd *= -1.0;
641 if(znabs < dtol) fd = 0.0;
642
643 fs -= zn * fd;
644 numexact++;
645 }
646
647 /* Return values of the source and dipole, normalized by area. */
648 fs /= panel->area;
649 fd /= panel->area;
650 if(pfd != NULL) *pfd = fd;
651
652
653 if(fs < 0.0) {
654 fprintf(stderr,
655 "\ncalcp: panel potential coeff. less than zero = %g\n", fs);
656 fprintf(stderr, "Okay = %d Evaluation Point = %g %g %g\n", okay, x, y, z);
657 fprintf(stderr, "Evaluation Point in local coords = %g %g %g\n",xn,yn, zn);
658 fprintf(stderr, "Panel Description Follows\n");
659 dp(panel);
660 /*exit(1);*/
661 }
662
663
664 return (fs);
665 }
666
667
dumpnums(flag,size)668 dumpnums(flag, size)
669 int flag, size;
670 {
671 double total;
672
673 if(flag == ON) { /* if first call */
674 num2ndsav = num2nd;
675 num4thsav = num4th;
676 numexactsav = numexact;
677 }
678 else {
679 total = num2ndsav + num4thsav + numexactsav;
680 #if MULDAT == ON
681 fprintf(stdout, "Potential coefficient counts\n multipole only:\n");
682 fprintf(stdout,
683 " 2nd order: %d %.3g%%; 4th: %d %.3g%%; Integral: %d %.3g%%\n",
684 num2nd, 100*(num2ndsav/total), num4th, 100*(num4thsav/total),
685 numexact, 100*(numexactsav/total));
686 #endif
687 total = num2nd + num4th + numexact;
688 #if MULDAT == ON
689 fprintf(stdout, " multipole plus adaptive:\n");
690 fprintf(stdout,
691 " 2nd order: %d %.3g%%; 4th: %d %.3g%%; Integral: %d %.3g%%\n",
692 num2nd, 100*(num2nd/total), num4th, 100*(num4th/total),
693 numexact, 100*(numexact/total));
694 #endif
695 fprintf(stdout, "Percentage of multiplies done by multipole: %.3g%%\n",
696 100*(size*size - total)/(size*size));
697 if(size*size == total)
698 fprintf(stdout, "Warning: no multipole acceleration\n");
699 }
700 }
701
tilelength(nq)702 double tilelength(nq)
703 charge *nq;
704 {
705 return nq->max_diag;
706 }
707
708
709
710 /*
711 Evaluation of moments of quadrilateral surface relative to
712 local system, array S(15). First initialize array
713 Note that S(2)=S(6)=0 due to transfer above
714 */
715
ComputeMoments(pp)716 ComputeMoments(pp)
717 charge *pp;
718 {
719 int order=MAXORDER;
720 int i, j, nside, N, M, N1, M1, M2, MN1, MN2;
721 double dx, dy, dxdy, dydx, x, y, z, SI, *xp, *yp, *xpn, *ypn;
722 double ypp[3];
723 static double *XP[4], *YP[4], **I;
724 static int maxorder = 0;
725 static double CS[16] = { 0.0, 1.0, 1.0, 1.5, 1.5, 3.75, 1.0, 3.0,
726 1.5, 7.5, 1.5, 1.5, 3.75, 1.5, 7.5, 3.75 };
727 double *multi, sumc, sums, sign, **createBinom(), **createPmn();
728 int m, n, r, halfn, flrm, ceilm, numterms, rterms;
729
730 /* Allocate temporary storage and initialize arrays. */
731 if(order > maxorder) {
732 for(i = 0; i < 4; i++) {
733 CALLOC(XP[i], (order+3), double, ON, AQ2P);
734 CALLOC(YP[i], (order+3), double, ON, AQ2P);
735 }
736 /* Allocate the euclidean moments matrix, Imn. */
737 CALLOC(I, (order+1), double*, ON, AQ2P);
738 for(i = 0; i <= order; i++) CALLOC(I[i], (order+1), double, ON, AQ2P);
739
740 maxorder = order;
741 }
742
743 /* First zero out the Moments matrix. */
744 for(i = 0; i <= order; i++) {
745 for(j = 0; j <= order; j++) {
746 I[i][j] = 0.0;
747 }
748 }
749
750 /* Compute powers of x and y at corner pts. */
751 for(i = 0; i < pp->shape; i++) {
752 xp = XP[i];
753 yp = YP[i];
754 xp[1] = pp->corner[i][XI];
755 yp[1] = pp->corner[i][YI];
756 for(j = 2; j <= order+2; j++) {
757 xp[j] = xp[j-1] * xp[1];
758 yp[j] = yp[j-1] * yp[1];
759 }
760 }
761
762 /* First moment, easy, just the panel area. */
763 I[0][0] = pp->area;
764
765 /* By using centroid, (1,0) and (0,1) are zero, so begin with (2,0). */
766 for(nside = 0; nside < pp->shape; nside++) {
767 xp = XP[nside];
768 yp = YP[nside];
769 if(nside == (pp->shape - 1)) {
770 xpn = XP[0];
771 ypn = YP[0];
772 }
773 else {
774 xpn = XP[nside + 1];
775 ypn = YP[nside + 1];
776 }
777
778 dx = xpn[1] - xp[1];
779 dy = ypn[1] - yp[1];
780
781 if(fabs(dx) >= fabs(dy)) {
782 dydx = dy/dx;
783 for(M = 2; M <= order; M++) {
784 M1 = M + 1;
785 M2 = M + 2;
786
787 SI = ((xpn[M1] * ypn[1]) - (xp[M1] * yp[1])) / M1
788 + dydx * (xp[M2] - xpn[M2]) / (M1 * M2);
789 I[M][0] += SI;
790
791 for(N = 1; N <= M; N++) {
792 N1 = N + 1;
793 MN1 = M - N + 1;
794 SI = (xpn[MN1] * ypn[N1] - xp[MN1] * yp[N1]) / (MN1 * N1)
795 - (dydx * N * SI) / MN1;
796 I[M-N][N] += SI;
797 }
798 }
799 }
800 else {
801 dxdy = dx/dy;
802 for(M = 2; M <= order; M++) {
803 M1 = M + 1;
804 M2 = M + 2;
805 SI = (dxdy / (M1 * M2)) * (ypn[M2] - yp[M2]);
806 I[0][M] += SI;
807 for(N = 1; N <= M; N++) {
808 MN1 = M - N + 1;
809 MN2 = MN1 + 1;
810 SI = dxdy * ((xpn[N] * ypn[MN2] - xp[N] * yp[MN2]) / (MN1 * MN2)
811 - (N * SI / MN1));
812 I[N][M-N] += SI;
813 }
814 }
815 }
816 }
817
818 /* Now Create the S vector for calcp. */
819 for(i = 0, M = 0; M <= 4; M++) {
820 for(N = 0; N <= (4 - M); N++) {
821 i++;
822 pp->moments[i] = I[M][N] * CS[i];
823 }
824 }
825
826 }
827
828 /* Debugging Print Routines follow. */
829
dp(panel)830 dp(panel)
831 charge *panel;
832 {
833 int i;
834 double c[4][3];
835
836 printf("shape=%d maxdiag=%g mindiag=%g area=%g\n",
837 panel->shape,
838 panel->max_diag, panel->min_diag, panel->area);
839
840 fprintf(stdout, "surface: `%s'\n", panel->surf->name);
841
842 printf("x=%g y=%g z=%g\n", panel->x, panel->y, panel->z);
843 printf("X= %g %g %g\n", panel->X[0], panel->X[1], panel->X[2]);
844 printf("Y= %g %g %g\n", panel->Y[0], panel->Y[1], panel->Y[2]);
845 printf("Z= %g %g %g\n", panel->Z[0], panel->Z[1], panel->Z[2]);
846
847 for(i=0; i < panel->shape; i++)
848 printf("corner%d = %g %g %g\n",
849 i, panel->corner[i][0], panel->corner[i][1], panel->corner[i][2]);
850
851 for(i = 0; i < panel->shape; i++) {
852 c[i][0] = panel->x + panel->corner[i][0]*panel->X[0]
853 + panel->corner[i][0]*panel->X[1] + panel->corner[i][0]*panel->X[2];
854 c[i][1] = panel->y + panel->corner[i][1]*panel->Y[0]
855 + panel->corner[i][1]*panel->Y[1] + panel->corner[i][1]*panel->Y[2];
856 c[i][2] = panel->z + panel->corner[i][2]*panel->Z[0]
857 + panel->corner[i][2]*panel->Z[1] + panel->corner[i][2]*panel->Z[2];
858 printf("absolute corner%d = %g %g %g\n", i, c[i][0], c[i][1], c[i][2]);
859 }
860
861 for(i=0; i < panel->shape; i++)
862 printf("length%d = %g\n", i, panel->length[i]);
863
864 printf("multipole coeffs: ");
865 for(i=0; i < 16; i++) {
866 printf("%g ", panel->moments[i]);
867 if( (i % 6) == 0) printf("\n");
868 }
869 printf("\n");
870 }
871
872
873
874 #define DIS 2
875 #define SCALE 5
876
testCalcp(pp)877 testCalcp(pp)
878 charge *pp;
879 {
880
881 double offx, offy, offz, x, y, z, mult;
882 int i, j, k;
883
884 offx = pp->x;
885 offy = pp->y;
886 offz = pp->z;
887
888 mult = 0.5 * pp->min_diag;
889
890 printf("\n\nCenter Point %g %g %g\n", offx, offy, offz);
891 for(i=0; i < DIS; i++) {
892 for(j=0; j < DIS; j++) {
893 for(k=0; k < DIS; k++) {
894 x = offx + i * mult * SCALE;
895 y = offy + j * mult * SCALE;
896 z = offz + k * mult * SCALE;
897 printf("Eval pt = %g %g %g\n", x, y, z);
898 calcp(pp, x, y, z, NULL);
899 }
900 }
901 }
902 }
903
904
fileCorners(pp,f)905 fileCorners(pp, f)
906 FILE *f;
907 charge *pp;
908
909 {
910 int i;
911
912 for(i=0; i < pp->shape; i++)
913 fprintf(f, "%g %g\n", pp->corner[i][0], pp->corner[i][1]);
914 }
915
916
917 /* Test the moment code. */
calcpm(multi,x,y,z,origorder,order)918 calcpm(multi, x, y, z, origorder, order)
919 double *multi, x, y, z;
920 int order;
921 {
922 charge panel, *ppanel;
923 double **mat, **mulMulti2P(), potential;
924 int i, numterms;
925
926 /* Create a temporary panel which has evaluation point as centroid. */
927 panel.x = x;
928 panel.y = y;
929 panel.z = z;
930 ppanel = &panel;
931
932 mat = mulMulti2P(0.0, 0.0, 0.0, &ppanel, 1, origorder);
933 numterms = multerms(origorder);
934
935 /* Calculate the potential. */
936 for(potential = 0.0, i = 0; i < numterms; i++) {
937 potential += mat[0][i] * multi[i];
938 }
939 }
940
941
942
943 /*
944
945 dismat(mat, size)
946 double **mat;
947 int size;
948 {
949 int i, j;
950
951 for(i=0; i <= size; i++) {
952 printf("\nrow i = ");
953 for(j=0; j <= size; j++) {
954 printf("%g ", mat[i][j]);
955 }
956 }
957 printf("\n");
958 }
959 */
960