ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/cv/bsdfmesh.c
Revision: 2.25
Committed: Mon Mar 24 06:07:46 2014 UTC (10 years ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.24: +3 -4 lines
Log Message:
Minor fixes

File Contents

# Content
1 #ifndef lint
2 static const char RCSid[] = "$Id: bsdfmesh.c,v 2.24 2014/03/15 18:11:37 greg Exp $";
3 #endif
4 /*
5 * Create BSDF advection mesh from radial basis functions.
6 *
7 * G. Ward
8 */
9
10 #ifndef _WIN32
11 #include <unistd.h>
12 #include <sys/wait.h>
13 #include <sys/mman.h>
14 #endif
15 #define _USE_MATH_DEFINES
16 #include <stdio.h>
17 #include <stdlib.h>
18 #include <string.h>
19 #include <math.h>
20 #include "bsdfrep.h"
21
22 #ifndef NEIGH_FACT2
23 #define NEIGH_FACT2 0.1 /* empirical neighborhood distance weight */
24 #endif
25 /* number of processes to run */
26 int nprocs = 1;
27 /* number of children (-1 in child) */
28 static int nchild = 0;
29
30 typedef struct {
31 int nrows, ncols; /* array size (matches migration) */
32 float *price; /* migration prices */
33 short *sord; /* sort for each row, low to high */
34 float *prow; /* current price row */
35 } PRICEMAT; /* sorted pricing matrix */
36
37 #define pricerow(p,i) ((p)->price + (i)*(p)->ncols)
38 #define psortrow(p,i) ((p)->sord + (i)*(p)->ncols)
39
40 /* Create a new migration holder (sharing memory for multiprocessing) */
41 static MIGRATION *
42 new_migration(RBFNODE *from_rbf, RBFNODE *to_rbf)
43 {
44 size_t memlen = sizeof(MIGRATION) +
45 sizeof(float)*(from_rbf->nrbf*to_rbf->nrbf - 1);
46 MIGRATION *newmig;
47 #ifdef _WIN32
48 if (nprocs > 1)
49 fprintf(stderr, "%s: warning - multiprocessing not supported\n",
50 progname);
51 nprocs = 1;
52 newmig = (MIGRATION *)malloc(memlen);
53 #else
54 if (nprocs <= 1) { /* single process? */
55 newmig = (MIGRATION *)malloc(memlen);
56 } else { /* else need to share memory */
57 newmig = (MIGRATION *)mmap(NULL, memlen, PROT_READ|PROT_WRITE,
58 MAP_ANON|MAP_SHARED, -1, 0);
59 if ((void *)newmig == MAP_FAILED)
60 newmig = NULL;
61 }
62 #endif
63 if (newmig == NULL) {
64 fprintf(stderr, "%s: cannot allocate new migration\n", progname);
65 exit(1);
66 }
67 newmig->rbfv[0] = from_rbf;
68 newmig->rbfv[1] = to_rbf;
69 /* insert in edge lists */
70 newmig->enxt[0] = from_rbf->ejl;
71 from_rbf->ejl = newmig;
72 newmig->enxt[1] = to_rbf->ejl;
73 to_rbf->ejl = newmig;
74 newmig->next = mig_list; /* push onto global list */
75 return(mig_list = newmig);
76 }
77
78 #ifdef _WIN32
79 #define await_children(n) (void)(n)
80 #define run_subprocess() 0
81 #define end_subprocess() (void)0
82 #else
83
84 /* Wait for the specified number of child processes to complete */
85 static void
86 await_children(int n)
87 {
88 int exit_status = 0;
89
90 if (n > nchild)
91 n = nchild;
92 while (n-- > 0) {
93 int status;
94 if (wait(&status) < 0) {
95 fprintf(stderr, "%s: missing child(ren)!\n", progname);
96 nchild = 0;
97 break;
98 }
99 --nchild;
100 if (status) { /* something wrong */
101 if ((status = WEXITSTATUS(status)))
102 exit_status = status;
103 else
104 exit_status += !exit_status;
105 fprintf(stderr, "%s: subprocess died\n", progname);
106 n = nchild; /* wait for the rest */
107 }
108 }
109 if (exit_status)
110 exit(exit_status);
111 }
112
113 /* Start child process if multiprocessing selected */
114 static pid_t
115 run_subprocess(void)
116 {
117 int status;
118 pid_t pid;
119
120 if (nprocs <= 1) /* any children requested? */
121 return(0);
122 await_children(nchild + 1 - nprocs); /* free up child process */
123 if ((pid = fork())) {
124 if (pid < 0) {
125 fprintf(stderr, "%s: cannot fork subprocess\n",
126 progname);
127 await_children(nchild);
128 exit(1);
129 }
130 ++nchild; /* subprocess started */
131 return(pid);
132 }
133 nchild = -1;
134 return(0); /* put child to work */
135 }
136
137 /* If we are in subprocess, call exit */
138 #define end_subprocess() if (nchild < 0) _exit(0); else
139
140 #endif /* ! _WIN32 */
141
142 /* Compute normalized distribution scattering functions for comparison */
143 static void
144 compute_nDSFs(const RBFNODE *rbf0, const RBFNODE *rbf1)
145 {
146 const double nf0 = (GRIDRES*GRIDRES) / rbf0->vtotal;
147 const double nf1 = (GRIDRES*GRIDRES) / rbf1->vtotal;
148 int x, y;
149 FVECT dv;
150
151 for (x = GRIDRES; x--; )
152 for (y = GRIDRES; y--; ) {
153 ovec_from_pos(dv, x, y); /* cube root (brightness) */
154 dsf_grid[x][y].val[0] = pow(nf0*eval_rbfrep(rbf0, dv), .3333);
155 dsf_grid[x][y].val[1] = pow(nf1*eval_rbfrep(rbf1, dv), .3333);
156 }
157 }
158
159 /* Compute neighborhood distance-squared (dissimilarity) */
160 static double
161 neighborhood_dist2(int x0, int y0, int x1, int y1)
162 {
163 int rad = GRIDRES>>5;
164 double sum2 = 0.;
165 double d;
166 int p[4];
167 int i, j;
168 /* check radius */
169 p[0] = x0; p[1] = y0; p[2] = x1; p[3] = y1;
170 for (i = 4; i--; ) {
171 if (p[i] < rad) rad = p[i];
172 if (GRIDRES-1-p[i] < rad) rad = GRIDRES-1-p[i];
173 }
174 for (i = -rad; i <= rad; i++)
175 for (j = -rad; j <= rad; j++) {
176 d = dsf_grid[x0+i][y0+j].val[0] -
177 dsf_grid[x1+i][y1+j].val[1];
178 sum2 += d*d;
179 }
180 return(sum2 / (4*rad*(rad+1) + 1));
181 }
182
183 /* Comparison routine needed for sorting price row */
184 static int
185 msrt_cmp(void *b, const void *p1, const void *p2)
186 {
187 PRICEMAT *pm = (PRICEMAT *)b;
188 float c1 = pm->prow[*(const short *)p1];
189 float c2 = pm->prow[*(const short *)p2];
190
191 if (c1 > c2) return(1);
192 if (c1 < c2) return(-1);
193 return(0);
194 }
195
196 /* Compute (and allocate) migration price matrix for optimization */
197 static void
198 price_routes(PRICEMAT *pm, const RBFNODE *from_rbf, const RBFNODE *to_rbf)
199 {
200 FVECT *vto = (FVECT *)malloc(sizeof(FVECT) * to_rbf->nrbf);
201 int i, j;
202
203 compute_nDSFs(from_rbf, to_rbf);
204 pm->nrows = from_rbf->nrbf;
205 pm->ncols = to_rbf->nrbf;
206 pm->price = (float *)malloc(sizeof(float) * pm->nrows*pm->ncols);
207 pm->sord = (short *)malloc(sizeof(short) * pm->nrows*pm->ncols);
208
209 if ((pm->price == NULL) | (pm->sord == NULL) | (vto == NULL)) {
210 fprintf(stderr, "%s: Out of memory in migration_costs()\n",
211 progname);
212 exit(1);
213 }
214 for (j = to_rbf->nrbf; j--; ) /* save repetitive ops. */
215 ovec_from_pos(vto[j], to_rbf->rbfa[j].gx, to_rbf->rbfa[j].gy);
216
217 for (i = from_rbf->nrbf; i--; ) {
218 const double from_ang = R2ANG(from_rbf->rbfa[i].crad);
219 FVECT vfrom;
220 short *srow;
221 ovec_from_pos(vfrom, from_rbf->rbfa[i].gx, from_rbf->rbfa[i].gy);
222 pm->prow = pricerow(pm,i);
223 srow = psortrow(pm,i);
224 for (j = to_rbf->nrbf; j--; ) {
225 double d; /* quadratic cost function */
226 d = Acos(DOT(vfrom, vto[j]));
227 pm->prow[j] = d*d;
228 d = R2ANG(to_rbf->rbfa[j].crad) - from_ang;
229 pm->prow[j] += d*d;
230 /* neighborhood difference */
231 pm->prow[j] += NEIGH_FACT2 * neighborhood_dist2(
232 from_rbf->rbfa[i].gx, from_rbf->rbfa[i].gy,
233 to_rbf->rbfa[j].gx, to_rbf->rbfa[j].gy );
234 srow[j] = j;
235 }
236 qsort_r(srow, pm->ncols, sizeof(short), pm, &msrt_cmp);
237 }
238 free(vto);
239 }
240
241 /* Free price matrix */
242 static void
243 free_routes(PRICEMAT *pm)
244 {
245 free(pm->price); pm->price = NULL;
246 free(pm->sord); pm->sord = NULL;
247 }
248
249 /* Compute minimum (optimistic) cost for moving the given source material */
250 static double
251 min_cost(double amt2move, const double *avail, const PRICEMAT *pm, int s)
252 {
253 const short *srow = psortrow(pm,s);
254 const float *prow = pricerow(pm,s);
255 double total_cost = 0;
256 int j;
257 /* move cheapest first */
258 for (j = 0; (j < pm->ncols) & (amt2move > FTINY); j++) {
259 int d = srow[j];
260 double amt = (amt2move < avail[d]) ? amt2move : avail[d];
261
262 total_cost += amt * prow[d];
263 amt2move -= amt;
264 }
265 return(total_cost);
266 }
267
268 typedef struct {
269 short s, d; /* source and destination */
270 float dc; /* discount to push inventory */
271 } ROWSENT; /* row sort entry */
272
273 /* Compare entries by discounted moving price */
274 static int
275 rmovcmp(void *b, const void *p1, const void *p2)
276 {
277 PRICEMAT *pm = (PRICEMAT *)b;
278 const ROWSENT *re1 = (const ROWSENT *)p1;
279 const ROWSENT *re2 = (const ROWSENT *)p2;
280 double price_diff;
281
282 if (re1->d < 0) return(re2->d >= 0);
283 if (re2->d < 0) return(-1);
284 price_diff = re1->dc*pricerow(pm,re1->s)[re1->d] -
285 re2->dc*pricerow(pm,re2->s)[re2->d];
286 if (price_diff > 0) return(1);
287 if (price_diff < 0) return(-1);
288 return(0);
289 }
290
291 /* Take a step in migration by choosing reasonable bucket to transfer */
292 static double
293 migration_step(MIGRATION *mig, double *src_rem, double *dst_rem, PRICEMAT *pm)
294 {
295 const int max2check = 100;
296 const double maxamt = 1./(double)pm->ncols;
297 const double minamt = maxamt*1e-4;
298 double *src_cost;
299 ROWSENT *rord;
300 struct {
301 int s, d; /* source and destination */
302 double price; /* cost per amount moved */
303 double amt; /* amount we can move */
304 } cur, best;
305 int r2check, i, ri;
306 /*
307 * Check cheapest available routes only -- a higher adjusted
308 * destination price implies that another source is closer, so
309 * we can hold off considering more expensive options until
310 * some other (hopefully better) moves have been made.
311 * A discount based on source remaining is supposed to prioritize
312 * movement from large lobes, but it doesn't seem to do much,
313 * so we have it set to 1.0 at the moment.
314 */
315 #define discount(qr) 1.0
316 /* most promising row order */
317 rord = (ROWSENT *)malloc(sizeof(ROWSENT)*pm->nrows);
318 if (rord == NULL)
319 goto memerr;
320 for (ri = pm->nrows; ri--; ) {
321 rord[ri].s = ri;
322 rord[ri].d = -1;
323 rord[ri].dc = 1.f;
324 if (src_rem[ri] <= minamt) /* enough source material? */
325 continue;
326 for (i = 0; i < pm->ncols; i++)
327 if (dst_rem[ rord[ri].d = psortrow(pm,ri)[i] ] > minamt)
328 break;
329 if (i >= pm->ncols) { /* moved all we can? */
330 free(rord);
331 return(.0);
332 }
333 rord[ri].dc = discount(src_rem[ri]);
334 }
335 if (pm->nrows > max2check) /* sort if too many sources */
336 qsort_r(rord, pm->nrows, sizeof(ROWSENT), pm, &rmovcmp);
337 /* allocate cost array */
338 src_cost = (double *)malloc(sizeof(double)*pm->nrows);
339 if (src_cost == NULL)
340 goto memerr;
341 for (i = pm->nrows; i--; ) /* starting costs for diff. */
342 src_cost[i] = min_cost(src_rem[i], dst_rem, pm, i);
343 /* find best source & dest. */
344 best.s = best.d = -1; best.price = FHUGE; best.amt = 0;
345 if ((r2check = pm->nrows) > max2check)
346 r2check = max2check; /* put a limit on search */
347 for (ri = 0; ri < r2check; ri++) { /* check each source row */
348 double cost_others = 0;
349 cur.s = rord[ri].s;
350 if ((cur.d = rord[ri].d) < 0 ||
351 rord[ri].dc*pricerow(pm,cur.s)[cur.d] >= best.price) {
352 if (pm->nrows > max2check) break; /* sorted end */
353 continue; /* else skip this one */
354 }
355 cur.amt = (src_rem[cur.s] < dst_rem[cur.d]) ?
356 src_rem[cur.s] : dst_rem[cur.d];
357 /* don't just leave smidgen */
358 if (cur.amt > maxamt*1.02) cur.amt = maxamt;
359 dst_rem[cur.d] -= cur.amt; /* add up opportunity costs */
360 for (i = pm->nrows; i--; )
361 if (i != cur.s)
362 cost_others += min_cost(src_rem[i], dst_rem, pm, i)
363 - src_cost[i];
364 dst_rem[cur.d] += cur.amt; /* undo trial move */
365 /* discount effective price */
366 cur.price = ( pricerow(pm,cur.s)[cur.d] + cost_others/cur.amt ) *
367 rord[ri].dc;
368 if (cur.price < best.price) /* are we better than best? */
369 best = cur;
370 }
371 free(src_cost); /* clean up */
372 free(rord);
373 if ((best.s < 0) | (best.d < 0)) /* nothing left to move? */
374 return(.0);
375 /* else make the actual move */
376 mtx_coef(mig,best.s,best.d) += best.amt;
377 src_rem[best.s] -= best.amt;
378 dst_rem[best.d] -= best.amt;
379 return(best.amt);
380 memerr:
381 fprintf(stderr, "%s: Out of memory in migration_step()\n", progname);
382 exit(1);
383 #undef discount
384 }
385
386 /* Compute and insert migration along directed edge (may fork child) */
387 static MIGRATION *
388 create_migration(RBFNODE *from_rbf, RBFNODE *to_rbf)
389 {
390 const double end_thresh = 5e-6;
391 PRICEMAT pmtx;
392 MIGRATION *newmig;
393 double *src_rem, *dst_rem;
394 double total_rem = 1., move_amt;
395 int i, j;
396 /* check if exists already */
397 for (newmig = from_rbf->ejl; newmig != NULL;
398 newmig = nextedge(from_rbf,newmig))
399 if (newmig->rbfv[1] == to_rbf)
400 return(NULL);
401 /* else allocate */
402 #ifdef DEBUG
403 fprintf(stderr, "Building path from (theta,phi) (%.1f,%.1f) ",
404 get_theta180(from_rbf->invec),
405 get_phi360(from_rbf->invec));
406 fprintf(stderr, "to (%.1f,%.1f) with %d x %d matrix\n",
407 get_theta180(to_rbf->invec),
408 get_phi360(to_rbf->invec),
409 from_rbf->nrbf, to_rbf->nrbf);
410 #endif
411 newmig = new_migration(from_rbf, to_rbf);
412 if (run_subprocess())
413 return(newmig); /* child continues */
414 price_routes(&pmtx, from_rbf, to_rbf);
415 src_rem = (double *)malloc(sizeof(double)*from_rbf->nrbf);
416 dst_rem = (double *)malloc(sizeof(double)*to_rbf->nrbf);
417 if ((src_rem == NULL) | (dst_rem == NULL)) {
418 fprintf(stderr, "%s: Out of memory in create_migration()\n",
419 progname);
420 exit(1);
421 }
422 /* starting quantities */
423 memset(newmig->mtx, 0, sizeof(float)*from_rbf->nrbf*to_rbf->nrbf);
424 for (i = from_rbf->nrbf; i--; )
425 src_rem[i] = rbf_volume(&from_rbf->rbfa[i]) / from_rbf->vtotal;
426 for (j = to_rbf->nrbf; j--; )
427 dst_rem[j] = rbf_volume(&to_rbf->rbfa[j]) / to_rbf->vtotal;
428
429 do { /* move a bit at a time */
430 move_amt = migration_step(newmig, src_rem, dst_rem, &pmtx);
431 total_rem -= move_amt;
432 } while ((total_rem > end_thresh) & (move_amt > 0));
433
434 for (i = from_rbf->nrbf; i--; ) { /* normalize final matrix */
435 double nf = rbf_volume(&from_rbf->rbfa[i]);
436 if (nf <= FTINY) continue;
437 nf = from_rbf->vtotal / nf;
438 for (j = to_rbf->nrbf; j--; )
439 mtx_coef(newmig,i,j) *= nf; /* row now sums to 1.0 */
440 }
441 end_subprocess(); /* exit here if subprocess */
442 free_routes(&pmtx); /* free working arrays */
443 free(src_rem);
444 free(dst_rem);
445 return(newmig);
446 }
447
448 /* Check if prospective vertex would create overlapping triangle */
449 static int
450 overlaps_tri(const RBFNODE *bv0, const RBFNODE *bv1, const RBFNODE *pv)
451 {
452 const MIGRATION *ej;
453 RBFNODE *vother[2];
454 int im_rev;
455 /* find shared edge in mesh */
456 for (ej = pv->ejl; ej != NULL; ej = nextedge(pv,ej)) {
457 const RBFNODE *tv = opp_rbf(pv,ej);
458 if (tv == bv0) {
459 im_rev = is_rev_tri(ej->rbfv[0]->invec,
460 ej->rbfv[1]->invec, bv1->invec);
461 break;
462 }
463 if (tv == bv1) {
464 im_rev = is_rev_tri(ej->rbfv[0]->invec,
465 ej->rbfv[1]->invec, bv0->invec);
466 break;
467 }
468 }
469 if (!get_triangles(vother, ej)) /* triangle on same side? */
470 return(0);
471 return(vother[im_rev] != NULL);
472 }
473
474 /* Find convex hull vertex to complete triangle (oriented call) */
475 static RBFNODE *
476 find_chull_vert(const RBFNODE *rbf0, const RBFNODE *rbf1)
477 {
478 FVECT vmid, vejn, vp;
479 RBFNODE *rbf, *rbfbest = NULL;
480 double dprod, area2, bestarea2 = FHUGE, bestdprod = -.5;
481
482 VSUB(vejn, rbf1->invec, rbf0->invec);
483 VADD(vmid, rbf0->invec, rbf1->invec);
484 if (normalize(vejn) == 0 || normalize(vmid) == 0)
485 return(NULL);
486 /* XXX exhaustive search */
487 /* Find triangle with minimum rotation from perpendicular */
488 for (rbf = dsf_list; rbf != NULL; rbf = rbf->next) {
489 if ((rbf == rbf0) | (rbf == rbf1))
490 continue;
491 tri_orient(vp, rbf0->invec, rbf1->invec, rbf->invec);
492 if (DOT(vp, vmid) <= FTINY)
493 continue; /* wrong orientation */
494 area2 = .25*DOT(vp,vp);
495 VSUB(vp, rbf->invec, vmid);
496 dprod = -DOT(vp, vejn);
497 VSUM(vp, vp, vejn, dprod); /* above guarantees non-zero */
498 dprod = DOT(vp, vmid) / VLEN(vp);
499 if (dprod <= bestdprod + FTINY*(1 - 2*(area2 < bestarea2)))
500 continue; /* found better already */
501 if (overlaps_tri(rbf0, rbf1, rbf))
502 continue; /* overlaps another triangle */
503 rbfbest = rbf;
504 bestdprod = dprod; /* new one to beat */
505 bestarea2 = area2;
506 }
507 return(rbfbest);
508 }
509
510 /* Create new migration edge and grow mesh recursively around it */
511 static void
512 mesh_from_edge(MIGRATION *edge)
513 {
514 MIGRATION *ej0, *ej1;
515 RBFNODE *tvert[2];
516
517 if (edge == NULL)
518 return;
519 /* triangle on either side? */
520 get_triangles(tvert, edge);
521 if (tvert[0] == NULL) { /* grow mesh on right */
522 tvert[0] = find_chull_vert(edge->rbfv[0], edge->rbfv[1]);
523 if (tvert[0] != NULL) {
524 if (tvert[0]->ord > edge->rbfv[0]->ord)
525 ej0 = create_migration(edge->rbfv[0], tvert[0]);
526 else
527 ej0 = create_migration(tvert[0], edge->rbfv[0]);
528 if (tvert[0]->ord > edge->rbfv[1]->ord)
529 ej1 = create_migration(edge->rbfv[1], tvert[0]);
530 else
531 ej1 = create_migration(tvert[0], edge->rbfv[1]);
532 mesh_from_edge(ej0);
533 mesh_from_edge(ej1);
534 }
535 } else if (tvert[1] == NULL) { /* grow mesh on left */
536 tvert[1] = find_chull_vert(edge->rbfv[1], edge->rbfv[0]);
537 if (tvert[1] != NULL) {
538 if (tvert[1]->ord > edge->rbfv[0]->ord)
539 ej0 = create_migration(edge->rbfv[0], tvert[1]);
540 else
541 ej0 = create_migration(tvert[1], edge->rbfv[0]);
542 if (tvert[1]->ord > edge->rbfv[1]->ord)
543 ej1 = create_migration(edge->rbfv[1], tvert[1]);
544 else
545 ej1 = create_migration(tvert[1], edge->rbfv[1]);
546 mesh_from_edge(ej0);
547 mesh_from_edge(ej1);
548 }
549 }
550 }
551
552 /* Add normal direction if missing */
553 static void
554 check_normal_incidence(void)
555 {
556 static FVECT norm_vec = {.0, .0, 1.};
557 const int saved_nprocs = nprocs;
558 RBFNODE *near_rbf, *mir_rbf, *rbf;
559 double bestd;
560 int n;
561
562 if (dsf_list == NULL)
563 return; /* XXX should be error? */
564 near_rbf = dsf_list;
565 bestd = input_orient*near_rbf->invec[2];
566 if (single_plane_incident) { /* ordered plane incidence? */
567 if (bestd >= 1.-2.*FTINY)
568 return; /* already have normal */
569 } else {
570 switch (inp_coverage) {
571 case INP_QUAD1:
572 case INP_QUAD2:
573 case INP_QUAD3:
574 case INP_QUAD4:
575 break; /* quadrilateral symmetry? */
576 default:
577 return; /* else we can interpolate */
578 }
579 for (rbf = near_rbf->next; rbf != NULL; rbf = rbf->next) {
580 const double d = input_orient*rbf->invec[2];
581 if (d >= 1.-2.*FTINY)
582 return; /* seems we have normal */
583 if (d > bestd) {
584 near_rbf = rbf;
585 bestd = d;
586 }
587 }
588 }
589 if (mig_list != NULL) { /* need to be called first */
590 fprintf(stderr, "%s: Late call to check_normal_incidence()\n",
591 progname);
592 exit(1);
593 }
594 #ifdef DEBUG
595 fprintf(stderr, "Interpolating normal incidence by mirroring (%.1f,%.1f)\n",
596 get_theta180(near_rbf->invec), get_phi360(near_rbf->invec));
597 #endif
598 /* mirror nearest incidence */
599 n = sizeof(RBFNODE) + sizeof(RBFVAL)*(near_rbf->nrbf-1);
600 mir_rbf = (RBFNODE *)malloc(n);
601 if (mir_rbf == NULL)
602 goto memerr;
603 memcpy(mir_rbf, near_rbf, n);
604 mir_rbf->ord = near_rbf->ord - 1; /* not used, I think */
605 mir_rbf->next = NULL;
606 mir_rbf->ejl = NULL;
607 rev_rbf_symmetry(mir_rbf, MIRROR_X|MIRROR_Y);
608 nprocs = 1; /* compute migration matrix */
609 if (create_migration(mir_rbf, near_rbf) == NULL)
610 exit(1); /* XXX should never happen! */
611 norm_vec[2] = input_orient; /* interpolate normal dist. */
612 rbf = e_advect_rbf(mig_list, norm_vec, 2*near_rbf->nrbf);
613 nprocs = saved_nprocs; /* final clean-up */
614 free(mir_rbf);
615 free(mig_list);
616 mig_list = near_rbf->ejl = NULL;
617 insert_dsf(rbf); /* insert interpolated normal */
618 return;
619 memerr:
620 fprintf(stderr, "%s: Out of memory in check_normal_incidence()\n",
621 progname);
622 exit(1);
623 }
624
625 /* Build our triangle mesh from recorded RBFs */
626 void
627 build_mesh(void)
628 {
629 double best2 = M_PI*M_PI;
630 RBFNODE *shrt_edj[2];
631 RBFNODE *rbf0, *rbf1;
632 /* add normal if needed */
633 check_normal_incidence();
634 /* check if isotropic */
635 if (single_plane_incident) {
636 for (rbf0 = dsf_list; rbf0 != NULL; rbf0 = rbf0->next)
637 if (rbf0->next != NULL)
638 create_migration(rbf0, rbf0->next);
639 await_children(nchild);
640 return;
641 }
642 shrt_edj[0] = shrt_edj[1] = NULL; /* start w/ shortest edge */
643 for (rbf0 = dsf_list; rbf0 != NULL; rbf0 = rbf0->next)
644 for (rbf1 = rbf0->next; rbf1 != NULL; rbf1 = rbf1->next) {
645 double dist2 = 2. - 2.*DOT(rbf0->invec,rbf1->invec);
646 if (dist2 < best2) {
647 shrt_edj[0] = rbf0;
648 shrt_edj[1] = rbf1;
649 best2 = dist2;
650 }
651 }
652 if (shrt_edj[0] == NULL) {
653 fprintf(stderr, "%s: Cannot find shortest edge\n", progname);
654 exit(1);
655 }
656 /* build mesh from this edge */
657 if (shrt_edj[0]->ord < shrt_edj[1]->ord)
658 mesh_from_edge(create_migration(shrt_edj[0], shrt_edj[1]));
659 else
660 mesh_from_edge(create_migration(shrt_edj[1], shrt_edj[0]));
661 /* complete migrations */
662 await_children(nchild);
663 }