1 #include "mulGlobal.h"
2 #include "zbufGlobal.h"
3
4 extern char *hack_path(char*);
5
6 #define XI 0
7 #define YI 1
8 #define ZI 2
9 #define EQUIV_TOL 1.0e-9
10 /*#define PLANAR_TOL 1.0e-3*/
11 #define PLANAR_TOL 1.0e-5
12 #define MAX_REAL 1.0e+20;
13
14 #define Dot_Product(V1,V2) (V1[XI]*V2[XI]+V1[YI]*V2[YI]+V1[ZI]*V2[ZI])
15 #define DotP_Product(V1,R,S,T) ((V1[XI])*(R)+(V1[YI])*(S)+(V1[ZI])*(T))
16
17 #define ONE3 0.3333333333333
18
19 /*
20 dumps a rows x cols matrix of doubles; assumes indices from zero
21 */
dumpCorners(fp,mat,rows,cols)22 void dumpCorners(fp, mat, rows, cols)
23 int rows, cols;
24 double **mat;
25 FILE *fp;
26 {
27 int i, j;
28 for(i = 0; i < rows; i++) {
29 fprintf(fp, " corner%d ", i);
30 for(j = 0; j < cols; j++) {
31 if(mat[i][j] < 0.0) fprintf(fp, "%.5e ", mat[i][j]);
32 else fprintf(fp, " %.5e ", mat[i][j]);
33 }
34 fprintf(fp, "\n");
35 }
36 }
37
38 /*
39 dumps state of important compile flags
40 */
dumpConfig(fp,name)41 void dumpConfig(fp, name)
42 char *name;
43 FILE *fp;
44 {
45 int size = -1; /* for '#define MAXITER size' case */
46
47 fprintf(fp, "\n%s CONFIGURATION FLAGS:\n", name);
48
49 fprintf(fp, " DISCRETIZATION CONFIGURATION\n");
50
51 fprintf(fp, " WRMETH");
52 if(WRMETH == COLLOC)
53 fprintf(fp, " == COLLOC (point collocation)\n");
54 else if(WRMETH == SUBDOM)
55 fprintf(fp, " == SUBDOM (not implemented - do collocation)\n");
56 else if(WRMETH == GALKIN)
57 fprintf(fp, " == GALKIN (not implemented - do collocation)\n");
58 fprintf(fp, " ELTYPE");
59 if(ELTYPE == CONST)
60 fprintf(fp, " == CONST (constant panel densities)\n");
61 else if(ELTYPE == AFFINE)
62 fprintf(fp, " == AFFINE (not implemented - use constant)\n");
63 else if(ELTYPE == QUADRA)
64 fprintf(fp, " == QUADRA (not implemented - use constant)\n");
65
66 fprintf(fp, " MULTIPOLE CONFIGURATION\n");
67
68 fprintf(fp, " DNTYPE");
69 if(DNTYPE == NOLOCL)
70 fprintf(fp, " == NOLOCL (no locals in dwnwd pass)\n");
71 else if(DNTYPE == NOSHFT)
72 fprintf(fp, " == NOSHFT (no local2local shift dwnwd pass)\n");
73 else if(DNTYPE == GRENGD)
74 fprintf(fp, " == GRENGD (full Greengard dwnwd pass)\n");
75 fprintf(fp, " MULTI");
76 if(MULTI == ON) fprintf(fp, " == ON (include multipole part of P*q)\n");
77 else fprintf(fp, " == OFF (don't use multipole part of P*q)\n");
78 fprintf(fp, " RADINTER");
79 if(RADINTER == ON)
80 fprintf(fp," == ON (allow parent level interaction list entries)\n");
81 else
82 fprintf(fp," == OFF (use only cube level interaction list entries)\n");
83 fprintf(fp, " NNBRS == %d (max distance to a nrst neighbor)\n", NNBRS);
84 fprintf(fp, " ADAPT");
85 if(ADAPT == ON)
86 fprintf(fp, " == ON (adaptive - no expansions in exact cubes)\n");
87 else fprintf(fp, " == OFF (not adaptive - expansions in all cubes)\n");
88 fprintf(fp, " OPCNT");
89 if(OPCNT == ON)
90 fprintf(fp, " == ON (count P*q ops - exit after mat build)\n");
91 else fprintf(fp, " == OFF (no P*q op count - iterate to convergence)\n");
92
93 fprintf(fp, " MAXDEP");
94 fprintf(fp,
95 " == %d (assume no more than %d partitioning levels are needed)\n",
96 MAXDEP, MAXDEP);
97
98 fprintf(fp, " NUMDPT");
99 fprintf(fp,
100 " == %d (do %d potential evaluations for each dielectric panel)\n",
101 NUMDPT, NUMDPT);
102
103 fprintf(fp, " LINEAR SYSTEM SOLUTION CONFIGURATION\n");
104
105 fprintf(fp, " ITRTYP");
106 if(ITRTYP == GCR)
107 fprintf(fp, " == GCR (generalized conjugate residuals)\n");
108 else if(ITRTYP == GMRES)
109 fprintf(fp, " == GMRES (generalized minimum residuals)\n");
110 else fprintf(fp, " == %d (not implemented - use GCR)\n", ITRTYP);
111
112 fprintf(fp, " PRECOND");
113 if(PRECOND == BD) {
114 fprintf(fp,
115 " == BD (use block diagonal preconditioner)\n");
116 }
117 else if(PRECOND == OL) {
118 fprintf(fp,
119 " == OL (use overlap preconditioner)\n");
120 }
121 else if(PRECOND == NONE) {
122 fprintf(fp,
123 " == NONE (no preconditioner)\n");
124 }
125 else fprintf(fp, " == %d (not implemented - use BD, OL or NONE)\n", PRECOND);
126
127 fprintf(fp, " DIRSOL");
128 if(DIRSOL == ON)
129 fprintf(fp, " == ON (do the whole calculation directly)\n");
130 else fprintf(fp, " == OFF (do the calculation iteratively)\n");
131
132 fprintf(fp, " EXPGCR");
133 if(EXPGCR == ON)
134 fprintf(fp, " == ON (do all P*q's explicitly w/full matrix)\n");
135 else fprintf(fp, " == OFF (do P*q's with multipole)\n");
136
137 fprintf(fp, " MAXITER");
138 if(MAXITER < 0) {
139 fprintf(fp, " == size (for n panel system, do at most n iterations)\n");
140 }
141 else fprintf(fp, " == %d (stop after %d iterations if not converged)\n",
142 MAXITER, MAXITER);
143
144 fprintf(fp, " EXRTSH");
145 fprintf(fp,
146 " == %g (exact/ttl cubes > %g on lowest level => stop refinement)\n",
147 EXRTSH, EXRTSH);
148 }
149
150 /*
151 dump the contents of a face struct
152 */
dump_face(fp,fac)153 void dump_face(fp, fac)
154 face *fac;
155 FILE *fp;
156 {
157 int i, j;
158 face **behind = fac->behind;
159
160 fprintf(fp, "Face %d, %d sides, depth %d, mark %d, greylev %g\n",
161 fac->index, fac->numsides, fac->depth, fac->mark, fac->greylev);
162 fprintf(fp, " plane: n = (%g %g %g) rhs = %g\n",
163 fac->normal[0], fac->normal[1], fac->normal[2], fac->rhs);
164 fprintf(fp, " behind [depth(index)]:");
165 for(i = 0; i < fac->numbehind; i++) {
166 fprintf(fp, " %d(%d)", behind[i]->depth, behind[i]->index);
167 if(i % 10 == 0 && i != 0) fprintf(fp, "\n");
168 }
169 i--;
170 if(!(i % 10 && i != 0)) fprintf(fp, "\n");
171 fprintf(fp, " Corners\n");
172 dumpCorners(fp, fac->c, fac->numsides, 3);
173 }
174
initcalcp(panel_list)175 initcalcp(panel_list)
176 charge *panel_list;
177 {
178 charge *pq, *npq;
179 double vtemp[3];
180 double length, maxlength, minlength, length20, length31, sum, sum2, delta;
181 double normalize();
182 int i, j, next;
183
184 for(i=0, pq = panel_list; pq != NULL; pq = pq->next) i++;
185
186 pq = panel_list;
187 while (pq != NULL) {
188 /* Calculate edge lengths. */
189 maxlength = 0.0;
190 minlength = MAX_REAL;
191 for(i=0; i < pq->shape; i++) {
192 if(i == (pq->shape -1)) next = 0;
193 else next = i + 1;
194 for(sum= 0, j = 0; j < 3; j++) {
195 delta = pq->corner[next][j] - pq->corner[i][j];
196 sum += delta * delta;
197 }
198 pq->length[i] = length = sqrt(sum);
199 maxlength = MAX(maxlength, length);
200 minlength = MIN(minlength, length);
201 }
202
203 /* Get diags and lengths. */
204 for(sum= 0, sum2 = 0, i = 0; i < 3; i++) {
205 pq->X[i] = delta = pq->corner[2][i] - pq->corner[0][i];
206 sum += delta * delta;
207 if(pq->shape == 3) pq->Y[i] = pq->corner[0][i] - pq->corner[1][i];
208 else {
209 pq->Y[i] = delta = pq->corner[3][i] - pq->corner[1][i];
210 sum2 += delta * delta;
211 }
212 }
213 length20 = sqrt(sum);
214 length31 = sqrt(sum2);
215
216 /* Check on lengths for quad. */
217 if(pq->shape == 3) {
218 pq->max_diag = maxlength;
219 pq->min_diag = minlength;
220 }
221 else {
222 length = MAX(length20, length31);
223 pq->max_diag = length;
224 pq->min_diag = MIN(length20, length31);
225 }
226
227 /* Z-axis is normal to two diags. */
228 Cross_Product(pq->X, pq->Y, pq->Z);
229 /*#if 1 == 0*/
230 if(flip_normal(pq)) {
231 for(i = 0; i < 3; i++) {
232 pq->Z[i] = -(pq->Z[i]); /* flip the normal */
233 pq->X[i] = -(pq->X[i]); /* flip the x direction */
234 /* interchange points 0 and 2 so that corner order will be
235 consistent with X flip (note this is OK for both quads and tris) */
236 vtemp[0] = pq->corner[0][i];
237 pq->corner[0][i] = pq->corner[2][i];
238 pq->corner[2][i] = vtemp[0];
239 }
240 /* interchange length entries in length vector */
241 vtemp[0] = pq->length[0];
242 pq->length[0] = pq->length[1];
243 pq->length[1] = vtemp[0];
244 if(pq->shape == 4) {
245 vtemp[0] = pq->length[2];
246 pq->length[2] = pq->length[3];
247 pq->length[3] = vtemp[0];
248 }
249 }
250 /*#endif*/
251 pq->area = 0.5 * normalize(pq->Z);
252 normalize(pq->X);
253
254 /* Real Y-axis is normal to X and Z (resulting system is left-handed). */
255 Cross_Product(pq->X, pq->Z, pq->Y);
256
257 /* Project the corner points into the plane defined by edge midpoints. */
258 if(planarize(pq) == FALSE) {
259 /* Planarization too drastic, crack into two triangles. */
260 CALLOC(npq, 1, charge, ON, AMSC);
261 npq->next = pq->next;
262 pq->next = npq;
263 npq->shape = 3;
264 pq->shape = 3;
265 npq->cond = pq->cond;
266 npq->surf = pq->surf;
267 VCOPY(npq->corner[0], pq->corner[0]);
268 VCOPY(npq->corner[1], pq->corner[2]);
269 VCOPY(npq->corner[2], pq->corner[3]);
270 continue;
271 }
272
273 /* Calculate the centroid. */
274 centroid(pq, length20);
275
276 /* Put corners in the newly defined coord system. */
277 for(i=0; i < pq->shape; i++) {
278 pq->corner[i][XI] -= pq->x;
279 pq->corner[i][YI] -= pq->y;
280 pq->corner[i][ZI] -= pq->z;
281 }
282 for(i=0; i < pq->shape; i++) {
283 vtemp[XI] = Dot_Product(pq->corner[i], pq->X);
284 vtemp[YI] = Dot_Product(pq->corner[i], pq->Y);
285 vtemp[ZI] = Dot_Product(pq->corner[i], pq->Z);
286 pq->corner[i][XI] = vtemp[XI];
287 pq->corner[i][YI] = vtemp[YI];
288 pq->corner[i][ZI] = vtemp[ZI];
289 if(fabs(pq->corner[i][ZI]) > (EQUIV_TOL * pq->min_diag)) {
290 printf("FATAL PROGRAM ERROR: renormalized z=%g\n", pq->corner[i][ZI]);
291 exit(0);
292 }
293 pq->corner[i][ZI] = 0.0;
294 }
295
296 /* Iterate for the next panel. */
297 pq = pq->next;
298
299 }
300
301 }
302
303 /* Calculates result_vector = vector1 X vector2. */
Cross_Product(vector1,vector2,result_vector)304 Cross_Product(vector1, vector2, result_vector)
305 double vector1[], vector2[], result_vector[];
306 {
307 result_vector[XI] = vector1[YI]*vector2[ZI] - vector1[ZI]*vector2[YI];
308 result_vector[YI] = vector1[ZI]*vector2[XI] - vector1[XI]*vector2[ZI];
309 result_vector[ZI] = vector1[XI]*vector2[YI] - vector1[YI]*vector2[XI];
310 }
311
normalize(vector)312 double normalize(vector)
313 double vector[3];
314 {
315 double length;
316 int i;
317
318 length = sqrt( vector[0]*vector[0]
319 + vector[1]*vector[1]
320 + vector[2]*vector[2]);
321
322 for (i=0; i<3; i++) vector[i] = vector[i] / length;
323
324 return length;
325 }
326
327 /*
328 Determines centroid of a panel (defined as the point which make the
329 first moments vanish. Calculation begins by projection into the
330 coordinate system defined by the panel normal as the z-axis and
331 edge02 as the x-axis.
332 */
centroid(pp,x2)333 centroid(pp, x2)
334 charge *pp;
335 double x2;
336 {
337 double vertex1[3], vertex3[3];
338 double sum, dl, x1, y1, x3, y3, xc, yc;
339 int i;
340
341 /* Use vertex 0 as the origin. */
342 for(i=0; i< 3; i++) {
343 vertex1[i] = pp->corner[1][i] - pp->corner[0][i];
344 if(pp->shape == 4) vertex3[i] = pp->corner[3][i] - pp->corner[0][i];
345 else vertex3[i] = pp->corner[2][i] - pp->corner[0][i];
346 }
347
348 /* Project into the panel axes. */
349 y1 = Dot_Product(vertex1, pp->Y);
350 y3 = Dot_Product(vertex3, pp->Y);
351 x1 = Dot_Product(vertex1, pp->X);
352 x3 = Dot_Product(vertex3, pp->X);
353
354 yc = ONE3 * (y1 + y3);
355 xc = ONE3 * (x2 + ((x1 * y1 - x3 * y3)/(y1 - y3)));
356
357 pp->x = pp->corner[0][XI] + xc * pp->X[XI] + yc * pp->Y[XI];
358 pp->y = pp->corner[0][YI] + xc * pp->X[YI] + yc * pp->Y[YI];
359 pp->z = pp->corner[0][ZI] + xc * pp->X[ZI] + yc * pp->Y[ZI];
360
361 }
362
363 /*
364 Changes the corner points so that they lie in the plane defined by the
365 panel diagonals and any midpoint of an edge.
366 */
planarize(pq)367 planarize(pq)
368 charge *pq;
369 {
370 double origin[3], corner[3], delta[4][3], px, py, dx, dy, dz;
371 int i, j, numcorners = pq->shape;
372 double tolsq = PLANAR_TOL * pq->min_diag;
373
374 tolsq *= tolsq;
375
376 /* Triangular panels are just fine already. */
377 if(numcorners != 4) return(TRUE);
378
379 /* Pick edge midpoint as origin. */
380 for(i=0; i < 3; i++) origin[i] = 0.5 * (pq->corner[1][i] + pq->corner[0][i]);
381
382 for(i=0; i < numcorners; i++) {
383 for(j=0; j < 3; j++) corner[j] = pq->corner[i][j] - origin[j];
384 px = Dot_Product(corner, pq->X);
385 py = Dot_Product(corner, pq->Y);
386
387 dx = px * pq->X[XI] + py * pq->Y[XI] + origin[XI] - pq->corner[i][XI];
388 dy = px * pq->X[YI] + py * pq->Y[YI] + origin[YI] - pq->corner[i][YI];
389 dz = px * pq->X[ZI] + py * pq->Y[ZI] + origin[ZI] - pq->corner[i][ZI];
390
391 delta[i][XI] = dx;
392 delta[i][YI] = dy;
393 delta[i][ZI] = dz;
394
395 /* Moved too much, crack the panel. */
396 if((dx * dx + dy * dy + dz * dz) > tolsq) return(FALSE);
397 }
398
399 /* Still Here? Must be planarizing. */
400 for(i=0; i < numcorners; i++) {
401 for(j=0; j < 3; j++) {
402 pq->corner[i][j] += delta[i][j];
403 }
404 }
405 return(TRUE);
406 }
407
408 /*
409 determine if normal needs to be flipped to get dielectric bdry cond right
410 - this function uses 0.0 as a breakpoint when really machine precision
411 weighted checks should be done (really not an issue if ref point far)
412 */
flip_normal(panel)413 flip_normal(panel)
414 charge *panel;
415 {
416 int i;
417 double x, y, z;
418 double ctr_minus_n[3], ctr_plus_n[3], norm_minus, norm_plus, norm, norm_sq;
419 surface *surf = panel->surf;
420 int ref_inside = surf->ref_inside, flip_normal;
421 double *ref = surf->ref, *normal, angle, norm_n;
422 char *surf_name = surf->name;
423
424 if(surf->type != DIELEC && surf->type != BOTH) return(FALSE);
425
426 /* get panel corner (relative to reference point) and normal */
427 x = panel->corner[0][0] - ref[0];
428 y = panel->corner[0][1] - ref[1];
429 z = panel->corner[0][2] - ref[2];
430 norm_sq = x*x + y*y + z*z;
431 norm = sqrt(norm_sq);
432 normal = panel->Z;
433 norm_n =
434 sqrt(normal[0]*normal[0] + normal[1]*normal[1] + normal[2]*normal[2]);
435
436 /* figure the dot product between the normal and the corner-ref point
437 - ref_inside and angle <= 90 => flip
438 - ref_inside and angle > 90 => no flip
439 - !ref_inside and angle <= 90 => no flip
440 - !ref_inside and angle > 90 => flip */
441 /* figure angle */
442 angle = (x*normal[0] + y*normal[1] + z*normal[2])/(norm*norm_n);
443 if(ref_inside && angle <= 0.0) flip_normal = TRUE;
444 else if(ref_inside && angle > 0.0) flip_normal = FALSE;
445 else if(!ref_inside && angle <= 0.0) flip_normal = FALSE;
446 else if(!ref_inside && angle > 0.0) flip_normal = TRUE;
447 else {
448 fprintf(stderr,
449 "flip_normal: inconclusive test for normal flipping\n");
450 fprintf(stderr, " Surface: %s\n", hack_path(surf_name));
451 fprintf(stderr, " Translation: (%g %g %g)\n", surf->trans[0],
452 surf->trans[1], surf->trans[2]);
453 fprintf(stderr, " Reference point: (%g %g %g)\n",
454 ref[0], ref[1], ref[2]);
455 fprintf(stderr, " Panel corner: (%g %g %g)\n",
456 panel->corner[0][0], panel->corner[0][1], panel->corner[0][2]);
457 fprintf(stderr, " Normal: (%g %g %g)\n",
458 normal[0], normal[1], normal[2]);
459 exit(0);
460 }
461
462 return(flip_normal);
463 }
464