ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/eplus_adduvf.c
Revision: 2.11
Committed: Sat Feb 15 01:31:36 2014 UTC (10 years, 2 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.10: +242 -84 lines
Log Message:
Trial implementation of subsurface sampling -- needs work

File Contents

# User Rev Content
1 greg 2.1 #ifndef lint
2 greg 2.11 static const char RCSid[] = "$Id: eplus_adduvf.c,v 2.10 2014/02/13 17:33:37 greg Exp $";
3 greg 2.1 #endif
4     /*
5     * Add User View Factors to EnergyPlus Input Data File
6     *
7     * G.Ward for LBNL
8     */
9    
10     #include <stdlib.h>
11 greg 2.3 #include "rtio.h"
12     #include "rtmath.h"
13     #include "random.h"
14 greg 2.1 #include "eplus_idf.h"
15     #include "triangulate.h"
16 greg 2.2 #include "rtprocess.h"
17 greg 2.1
18     #ifndef NSAMPLES
19 greg 2.7 #define NSAMPLES 80000 /* default number of samples */
20 greg 2.1 #endif
21    
22 greg 2.11 #define SURF_EPS 0.0005 /* surface testing epsilon */
23    
24 greg 2.2 char *progname; /* global argv[0] */
25    
26 greg 2.7 int nsamps = NSAMPLES; /* number of samples to use */
27    
28 greg 2.5 char temp_octree[128]; /* temporary octree */
29 greg 2.2
30     const char UVF_PNAME[] =
31 greg 2.10 "ZoneProperty:UserViewFactors:bySurfaceName";
32 greg 2.1
33 greg 2.11 const char SUBSURF_PNAME[] =
34     "FenestrationSurface:Detailed";
35    
36 greg 2.1 const char ADD_HEADER[] =
37 greg 2.5 "\n!+++ User View Factors computed by Radiance +++!\n\n";
38 greg 2.1
39     #define NAME_FLD 1 /* name field always first? */
40    
41 greg 2.11 #define SS_BASE_FLD 4 /* subsurface base surface */
42     #define SS_VERT_FLD 10 /* subsurface vertex count */
43    
44 greg 2.1 typedef struct {
45 greg 2.9 const char *pname; /* object type name */
46 greg 2.1 short zone_fld; /* zone field index */
47     short vert_fld; /* vertex field index */
48     } SURF_PTYPE; /* surface type we're interested in */
49    
50     const SURF_PTYPE surf_type[] = {
51     {"BuildingSurface:Detailed", 4, 10},
52 greg 2.5 {"Floor:Detailed", 3, 9},
53     {"RoofCeiling:Detailed", 3, 9},
54     {"Wall:Detailed", 3, 9},
55 greg 2.1 {NULL}
56 greg 2.7 }; /* IDF surface types */
57 greg 2.1
58 greg 2.2 typedef struct s_zone {
59 greg 2.1 const char *zname; /* zone name */
60     struct s_zone *next; /* next zone in list */
61 greg 2.2 int nsurf; /* surface count */
62 greg 2.11 int ntotal; /* surfaces+subsurfaces */
63 greg 2.9 IDF_OBJECT *pfirst; /* first matching object */
64 greg 2.11 IDF_OBJECT *plast; /* last before subsurfaces */
65 greg 2.2 } ZONE; /* a list of collected zone surfaces */
66    
67     ZONE *zone_list = NULL; /* our list of zones */
68 greg 2.1
69 greg 2.11 LUTAB zonesurf_lut = /* zone surface lookup table */
70     LU_SINIT(NULL,NULL);
71    
72 greg 2.1 IDF_LOADED *our_idf = NULL; /* loaded/modified IDF */
73    
74 greg 2.7 typedef struct {
75 greg 2.11 FVECT norm; /* surface normal */
76     double area; /* surface area */
77     int nv; /* number of vertices */
78     FVECT vl[3]; /* vertex list (extends struct) */
79     } SURFACE; /* a polygonal surface */
80    
81     typedef struct {
82 greg 2.7 FVECT sdir[3]; /* UVW unit sampling vectors */
83     double poff; /* W-offset for plane of polygon */
84     double area_left; /* area left to sample */
85     int samp_left; /* remaining samples */
86     int wd; /* output file descriptor */
87     } POLYSAMP; /* structure for polygon sampling */
88    
89 greg 2.1 /* Create a new zone and push to top of our list */
90 greg 2.2 static ZONE *
91 greg 2.9 new_zone(const char *zname, IDF_OBJECT *param)
92 greg 2.1 {
93 greg 2.2 ZONE *znew = (ZONE *)malloc(sizeof(ZONE));
94 greg 2.1
95     if (znew == NULL)
96     return(NULL);
97     znew->zname = zname; /* assumes static */
98     znew->next = zone_list;
99     znew->pfirst = znew->plast = param;
100 greg 2.11 znew->ntotal = znew->nsurf = 1;
101 greg 2.1 return(zone_list = znew);
102     }
103    
104     /* Add the detailed surface (polygon) to the named zone */
105 greg 2.2 static ZONE *
106 greg 2.9 add2zone(IDF_OBJECT *param, const char *zname)
107 greg 2.1 {
108 greg 2.11 IDF_FIELD *nfp = idf_getfield(param, NAME_FLD);
109     ZONE *zptr;
110     LUENT *lep;
111 greg 2.1
112 greg 2.11 if (nfp == NULL) {
113     fputs(progname, stderr);
114     fputs(": surface missing name field!\n", stderr);
115     return(NULL);
116     }
117 greg 2.1 for (zptr = zone_list; zptr != NULL; zptr = zptr->next)
118     if (!strcmp(zptr->zname, zname))
119     break;
120 greg 2.11 if (zptr == NULL) {
121     zptr = new_zone(zname, param);
122     } else { /* keep surfaces together */
123     if (!idf_movobject(our_idf, param, zptr->plast))
124     return(NULL);
125     zptr->plast = param;
126     zptr->nsurf++;
127     zptr->ntotal++;
128     }
129     /* add to lookup table */
130     lep = lu_find(&zonesurf_lut, nfp->val);
131     if (lep == NULL)
132     return(NULL);
133     if (lep->data != NULL) {
134     fputs(progname, stderr);
135     fputs(": duplicate surface name '", stderr);
136     fputs(nfp->val, stderr);
137     fputs("'\n", stderr);
138     return(NULL);
139     }
140     lep->key = nfp->val;
141     lep->data = (char *)zptr;
142     return(zptr);
143     }
144    
145     /* Add a subsurface by finding its base surface and the corresponding zone */
146     static ZONE *
147     add_subsurf(IDF_OBJECT *param)
148     {
149     IDF_FIELD *bfp = idf_getfield(param, SS_BASE_FLD);
150     ZONE *zptr;
151     LUENT *lep;
152    
153     if (bfp == NULL) {
154     fputs(progname, stderr);
155     fputs(": missing base field name in subsurface!\n", stderr);
156     return(NULL);
157     }
158     lep = lu_find(&zonesurf_lut, bfp->val); /* find base zone */
159     if (lep == NULL || lep->data == NULL) {
160     fputs(progname, stderr);
161     fputs(": cannot find referenced base surface '", stderr);
162     fputs(bfp->val, stderr);
163     fputs("'\n", stderr);
164     return(NULL);
165     }
166     zptr = (ZONE *)lep->data; /* add this subsurface */
167 greg 2.9 if (!idf_movobject(our_idf, param, zptr->plast))
168 greg 2.1 return(NULL);
169 greg 2.11 zptr->ntotal++;
170 greg 2.1 return(zptr);
171     }
172    
173 greg 2.9 /* Return field for vertices in the given object */
174 greg 2.3 static IDF_FIELD *
175 greg 2.9 get_vlist(IDF_OBJECT *param, const char *zname)
176 greg 2.2 {
177     int i = 0;
178 greg 2.11 /* check if subsurface */
179     if (!strcmp(param->pname, SUBSURF_PNAME)) {
180     if (zname != NULL) {
181     LUENT *lep = lu_find(&zonesurf_lut,
182     idf_getfield(param,SS_BASE_FLD)->val);
183     if (strcmp((*(ZONE *)lep->data).zname, zname))
184     return(NULL);
185     }
186     return(idf_getfield(param, SS_VERT_FLD));
187     }
188 greg 2.2 /* check for surface type */
189     while (strcmp(surf_type[i].pname, param->pname))
190     if (surf_type[++i].pname == NULL)
191 greg 2.3 return(NULL);
192    
193     if (zname != NULL) { /* matches specified zone? */
194 greg 2.11 IDF_FIELD *fptr = idf_getfield(param, surf_type[i].zone_fld);
195 greg 2.3 if (fptr == NULL || strcmp(fptr->val, zname))
196     return(NULL);
197     }
198 greg 2.2 /* return field for #verts */
199 greg 2.3 return(idf_getfield(param, surf_type[i].vert_fld));
200 greg 2.2 }
201    
202 greg 2.11 /* Get/allocate surface polygon */
203     static SURFACE *
204     get_surface(IDF_FIELD *fptr)
205     {
206     SURFACE *surf;
207     int nv;
208     FVECT e1, e2, vc;
209     int i, j;
210     /* get number of vertices */
211     if (fptr == NULL || (nv = atoi(fptr->val)) < 3) {
212     fputs(progname, stderr);
213     fputs(": bad vertex count for surface!\n", stderr);
214     return(NULL);
215     }
216     surf = (SURFACE *)malloc(sizeof(SURFACE)+sizeof(FVECT)*(nv-3));
217     if (surf == NULL)
218     return(NULL);
219     surf->nv = nv;
220     for (i = nv; i--; ) /* reverse vertex order/orientation */
221     for (j = 0; j < 3; j++) {
222     fptr = fptr->next;
223     if (fptr == NULL) {
224     fputs(progname, stderr);
225     fputs(": missing vertex in get_surface()\n", stderr);
226     free(surf);
227     return(NULL);
228     }
229     if ((surf->vl[i][j] = atof(fptr->val)) == 0 &&
230     !isflt(fptr->val)) {
231     fputs(progname, stderr);
232     fputs(": bad vertex in get_surface()\n", stderr);
233     free(surf);
234     return(NULL);
235     }
236     }
237     /* compute area and normal */
238     surf->norm[0] = surf->norm[1] = surf->norm[2] = 0;
239     VSUB(e1, surf->vl[1], surf->vl[0]);
240     for (i = 2; i < nv; i++) {
241     VSUB(e2, surf->vl[i], surf->vl[0]);
242     fcross(vc, e1, e2);
243     surf->norm[0] += vc[0];
244     surf->norm[1] += vc[1];
245     surf->norm[2] += vc[2];
246     VCOPY(e1, e2);
247     }
248     surf->area = .5 * normalize(surf->norm);
249     if (surf->area == 0) {
250     fputs(progname, stderr);
251     fputs(": degenerate polygon in get_surface()\n", stderr);
252     free(surf);
253     return(NULL);
254     }
255     return(surf);
256     }
257    
258 greg 2.2 /* Convert surface to Radiance with modifier based on unique name */
259     static int
260 greg 2.9 rad_surface(IDF_OBJECT *param, FILE *ofp)
261 greg 2.2 {
262 greg 2.3 const char *sname = idf_getfield(param,NAME_FLD)->val;
263     IDF_FIELD *fptr = get_vlist(param, NULL);
264 greg 2.11 const char **fargs;
265     int nvert, i, j;
266 greg 2.2
267 greg 2.3 if (fptr == NULL || (nvert = atoi(fptr->val)) < 3) {
268 greg 2.2 fprintf(stderr, "%s: bad surface '%s'\n", progname, sname);
269     return(0);
270     }
271 greg 2.11 fargs = (const char **)malloc(sizeof(const char *)*3*nvert);
272     if (fargs == NULL)
273     return(0);
274 greg 2.2 fprintf(ofp, "\nvoid glow '%s'\n0\n0\n4 1 1 1 0\n", sname);
275     fprintf(ofp, "\n'%s' polygon 's_%s'\n0\n0\n%d\n", sname, sname, 3*nvert);
276 greg 2.11 for (j = nvert; j--; ) { /* get arguments in reverse */
277     for (i = 0; i < 3; i++) {
278 greg 2.2 fptr = fptr->next;
279 greg 2.3 if (fptr == NULL || !isflt(fptr->val)) {
280 greg 2.2 fprintf(stderr,
281 greg 2.3 "%s: missing/bad vertex for %s '%s'\n",
282     progname, param->pname, sname);
283 greg 2.2 return(0);
284     }
285 greg 2.11 fargs[3*j + i] = fptr->val;
286     }
287     }
288     for (j = 0; j < nvert; j++) { /* output reversed verts */
289     for (i = 0; i < 3; i++) {
290 greg 2.2 fputc('\t', ofp);
291 greg 2.11 fputs(fargs[3*j + i], ofp);
292 greg 2.2 }
293     fputc('\n', ofp);
294     }
295 greg 2.11 free(fargs);
296     return(!ferror(ofp));
297     }
298    
299     /* Convert subsurface to Radiance with modifier based on unique name */
300     static int
301     rad_subsurface(IDF_OBJECT *param, FILE *ofp)
302     {
303     const char *sname = idf_getfield(param,NAME_FLD)->val;
304     SURFACE *surf = get_surface(idf_getfield(param,SS_VERT_FLD));
305     int i;
306    
307     if (surf == NULL) {
308     fprintf(stderr, "%s: bad subsurface '%s'\n", progname, sname);
309     return(0);
310     }
311     fprintf(ofp, "\nvoid glow '%s'\n0\n0\n4 1 1 1 0\n", sname);
312     fprintf(ofp, "\n'%s' polygon 'ss_%s'\n0\n0\n%d\n",
313     sname, sname, 3*surf->nv);
314     for (i = 0; i < surf->nv; i++) { /* offset vertices */
315     FVECT vert;
316     VSUM(vert, surf->vl[i], surf->norm, 2.*SURF_EPS);
317     fprintf(ofp, "\t%.12f %.12f %.12f\n", vert[0], vert[1], vert[2]);
318     }
319     free(surf);
320 greg 2.2 return(!ferror(ofp));
321     }
322    
323 greg 2.3 /* Start rcontrib process */
324 greg 2.2 static int
325     start_rcontrib(SUBPROC *pd, ZONE *zp)
326     {
327 greg 2.11 #define BASE_AC 6
328 greg 2.2 static char *base_av[BASE_AC] = {
329 greg 2.11 "rcontrib", "-V+", "-fff", "-h", "-x", "1"
330 greg 2.2 };
331     char cbuf[300];
332     char **av;
333     FILE *ofp;
334 greg 2.9 IDF_OBJECT *pptr;
335 greg 2.11 IDF_FIELD *fptr;
336 greg 2.2 int i, n;
337     /* start oconv command */
338     sprintf(cbuf, "oconv - > '%s'", temp_octree);
339     if ((ofp = popen(cbuf, "w")) == NULL) {
340     fputs(progname, stderr);
341     fputs(": cannot open oconv process\n", stderr);
342     return(0);
343     }
344     /* allocate argument list */
345 greg 2.11 av = (char **)malloc(sizeof(char *)*(BASE_AC+4+2*(zp->ntotal)));
346 greg 2.2 if (av == NULL)
347     return(0);
348     for (i = 0; i < BASE_AC; i++)
349     av[i] = base_av[i];
350 greg 2.7 sprintf(cbuf, "%d", nsamps);
351 greg 2.2 av[i++] = "-c";
352 greg 2.11 av[i++] = cbuf; /* add modifiers & surfaces */
353     for (n = 0, pptr = zp->pfirst; n < zp->nsurf; n++, pptr = pptr->dnext) {
354     fptr = idf_getfield(pptr,NAME_FLD);
355 greg 2.3 if (fptr == NULL || !fptr->val[0]) {
356 greg 2.2 fputs(progname, stderr);
357 greg 2.9 fputs(": missing name for surface object\n", stderr);
358 greg 2.2 return(0);
359     }
360 greg 2.3 if (!rad_surface(pptr, ofp)) /* add surface to octree */
361     return(0);
362 greg 2.2 av[i++] = "-m";
363 greg 2.3 av[i++] = fptr->val;
364 greg 2.2 }
365 greg 2.11 /* now subsurfaces */
366     for ( ; n < zp->ntotal; n++, pptr = pptr->dnext) {
367     fptr = idf_getfield(pptr,NAME_FLD);
368     if (fptr == NULL || !fptr->val[0]) {
369     fputs(progname, stderr);
370     fputs(": missing name for subsurface object\n", stderr);
371     return(0);
372     }
373     if (!rad_subsurface(pptr, ofp)) /* add surface to octree */
374     return(0);
375     av[i++] = "-m";
376     av[i++] = fptr->val;
377     }
378 greg 2.2 if (pclose(ofp) != 0) { /* finish oconv */
379     fputs(progname, stderr);
380     fputs(": error running oconv process\n", stderr);
381     return(0);
382     }
383     av[i++] = temp_octree; /* add final octree argument */
384     av[i] = NULL;
385 greg 2.8 if (open_process(pd, av) <= 0) { /* start process */
386 greg 2.2 fputs(progname, stderr);
387     fputs(": cannot start rcontrib process\n", stderr);
388     return(0);
389     }
390     free(av); /* all done -- clean up */
391     return(1);
392     #undef BASE_AC
393     }
394    
395 greg 2.3 /* Initialize polygon sampling */
396     static Vert2_list *
397 greg 2.11 init_poly(POLYSAMP *ps, IDF_FIELD *f0)
398 greg 2.3 {
399 greg 2.11 SURFACE *surf;
400     Vert2_list *vl2;
401     int i;
402     /* get 3-D polygon vertices */
403     if ((surf = get_surface(f0)) == NULL)
404     return(NULL);
405     vl2 = polyAlloc(surf->nv);
406 greg 2.3 if (vl2 == NULL)
407     return(NULL);
408 greg 2.11 /* create X & Y axes */
409     VCOPY(ps->sdir[2], surf->norm);
410     VSUB(ps->sdir[0], surf->vl[1], surf->vl[0]);
411     if (normalize(ps->sdir[0]) == 0)
412 greg 2.3 return(NULL);
413     fcross(ps->sdir[1], ps->sdir[2], ps->sdir[0]);
414     /* compute plane offset */
415 greg 2.11 ps->poff = DOT(surf->vl[0], ps->sdir[2]);
416 greg 2.3 /* assign 2-D vertices */
417 greg 2.11 for (i = 0; i < surf->nv; i++) {
418     vl2->v[i].mX = DOT(surf->vl[i], ps->sdir[0]);
419     vl2->v[i].mY = DOT(surf->vl[i], ps->sdir[1]);
420     }
421     ps->area_left = surf->area;
422     free(surf); /* it's ready! */
423     vl2->p = (void *)ps;
424 greg 2.3 return(vl2);
425     }
426    
427     /* Generate samples on 2-D triangle */
428     static int
429     sample_triangle(const Vert2_list *vl2, int a, int b, int c)
430     {
431     POLYSAMP *ps = (POLYSAMP *)vl2->p;
432     float *samp;
433     FVECT orig;
434     FVECT ab, ac;
435     double area;
436     int i, j, ns;
437     /* compute sampling axes */
438     for (i = 3; i--; ) {
439     orig[i] = vl2->v[a].mX*ps->sdir[0][i] +
440     vl2->v[a].mY*ps->sdir[1][i] +
441 greg 2.11 (ps->poff+SURF_EPS)*ps->sdir[2][i];
442 greg 2.3 ab[i] = (vl2->v[b].mX - vl2->v[a].mX)*ps->sdir[0][i] +
443     (vl2->v[b].mY - vl2->v[a].mY)*ps->sdir[1][i];
444     ac[i] = (vl2->v[c].mX - vl2->v[a].mX)*ps->sdir[0][i] +
445     (vl2->v[c].mY - vl2->v[a].mY)*ps->sdir[1][i];
446     }
447     /* compute number of samples to take */
448     area = .5*(vl2->v[a].mX*vl2->v[b].mY - vl2->v[b].mX*vl2->v[a].mY +
449     vl2->v[b].mX*vl2->v[c].mY - vl2->v[c].mX*vl2->v[b].mY +
450     vl2->v[c].mX*vl2->v[a].mY - vl2->v[a].mX*vl2->v[c].mY);
451     if (area < .0) {
452     fputs(progname, stderr);
453     fputs(": negative triangle area in sample_triangle()\n", stderr);
454     return(0);
455     }
456     if (area >= ps->area_left) {
457     ns = ps->samp_left;
458     ps->area_left = 0;
459     } else {
460     ns = (ps->samp_left*area/ps->area_left + .5);
461     ps->samp_left -= ns;
462     ps->area_left -= area;
463     }
464     if (ns <= 0) /* XXX should be error? */
465     return(1);
466     /* buffer sample rays */
467     samp = (float *)malloc(sizeof(float)*6*ns);
468     if (samp == NULL)
469     return(0);
470 greg 2.4 for (i = ns; i--; ) { /* stratified Monte Carlo sampling */
471 greg 2.3 double sv[4];
472 greg 2.5 FVECT dv;
473 greg 2.3 multisamp(sv, 4, (i+frandom())/(double)ns);
474     sv[0] *= sv[1] = sqrt(sv[1]);
475     sv[1] = 1. - sv[1];
476     for (j = 3; j--; )
477 greg 2.5 samp[i*6 + j] = orig[j] + sv[0]*ab[j] + sv[1]*ac[j];
478     sv[2] = sqrt(sv[2]);
479     sv[3] *= 2.*PI;
480     dv[0] = tcos(sv[3]) * sv[2];
481     dv[1] = tsin(sv[3]) * sv[2];
482     dv[2] = sqrt(1. - sv[2]*sv[2]);
483     for (j = 3; j--; )
484     samp[i*6 + 3 + j] = dv[0]*ps->sdir[0][j] +
485     dv[1]*ps->sdir[1][j] +
486     dv[2]*ps->sdir[2][j] ;
487 greg 2.3 }
488     /* send to our process */
489     writebuf(ps->wd, (char *)samp, sizeof(float)*6*ns);
490     free(samp); /* that's it! */
491     return(1);
492     }
493    
494     /* Sample the given surface */
495     static int
496 greg 2.9 sample_surface(IDF_OBJECT *param, int wd)
497 greg 2.3 {
498     POLYSAMP psamp;
499     int nv;
500     Vert2_list *vlist2;
501     /* set up our polygon sampler */
502 greg 2.11 if ((vlist2 = init_poly(&psamp, get_vlist(param, NULL))) == NULL) {
503 greg 2.5 fprintf(stderr, "%s: bad polygon %s '%s'\n",
504 greg 2.3 progname, param->pname,
505     idf_getfield(param,NAME_FLD)->val);
506     return(0);
507     }
508 greg 2.7 psamp.samp_left = nsamps; /* assign samples & destination */
509 greg 2.3 psamp.wd = wd;
510 greg 2.11 /* hack for subsurface sampling */
511     psamp.poff += 2.*SURF_EPS * !strcmp(param->pname, SUBSURF_PNAME);
512 greg 2.3 /* sample each subtriangle */
513     if (!polyTriangulate(vlist2, &sample_triangle))
514     return(0);
515     polyFree(vlist2); /* clean up and return */
516     return(1);
517     }
518    
519 greg 2.2 /* Compute User View Factors using open rcontrib process */
520     static int
521 greg 2.3 compute_uvfs(SUBPROC *pd, ZONE *zp)
522 greg 2.2 {
523 greg 2.9 IDF_OBJECT *pptr, *pout, *pptr1;
524 greg 2.3 float *uvfa;
525     char uvfbuf[24];
526     int n, m;
527 greg 2.9 /* create output object */
528     pout = idf_newobject(our_idf, UVF_PNAME,
529 greg 2.11 " ! computed by Radiance\n ", our_idf->plast);
530 greg 2.3 if (pout == NULL) {
531     fputs(progname, stderr);
532 greg 2.9 fputs(": cannot create new IDF object\n", stderr);
533 greg 2.3 return(0);
534     }
535     if (!idf_addfield(pout, zp->zname,
536     " !- Zone Name\n ")) {
537     fputs(progname, stderr);
538     fputs(": cannot add zone name field\n", stderr);
539     return(0);
540     }
541     /* allocate read buffer */
542 greg 2.11 uvfa = (float *)malloc(sizeof(float)*3*zp->ntotal);
543 greg 2.3 if (uvfa == NULL)
544     return(0);
545     /* UVFs from each surface */
546 greg 2.11 for (n = 0, pptr = zp->pfirst; n < zp->ntotal; n++, pptr = pptr->dnext) {
547 greg 2.4 double vfsum = 0;
548 greg 2.3 /* send samples to rcontrib */
549     if (!sample_surface(pptr, pd->w))
550     return(0);
551     /* read results */
552 greg 2.11 if (readbuf(pd->r, (char *)uvfa, sizeof(float)*3*zp->ntotal) !=
553     sizeof(float)*3*zp->ntotal) {
554 greg 2.3 fputs(progname, stderr);
555     fputs(": read error from rcontrib process\n", stderr);
556     return(0);
557     }
558     /* append UVF fields */
559     for (m = 0, pptr1 = zp->pfirst;
560 greg 2.11 m < zp->ntotal; m++, pptr1 = pptr1->dnext) {
561 greg 2.4 vfsum += uvfa[3*m + 1];
562 greg 2.11 if (pptr1 == pptr && uvfa[3*m + 1] > .001)
563     fprintf(stderr,
564 greg 2.5 "%s: warning - non-zero self-VF (%.1f%%) for surface '%s'\n",
565     progname, 100.*uvfa[3*m + 1],
566 greg 2.3 idf_getfield(pptr,NAME_FLD)->val);
567 greg 2.6 sprintf(uvfbuf, "%.4f", uvfa[3*m + 1]);
568 greg 2.3 if (!idf_addfield(pout,
569     idf_getfield(pptr,NAME_FLD)->val, NULL) ||
570     !idf_addfield(pout,
571     idf_getfield(pptr1,NAME_FLD)->val, NULL) ||
572     !idf_addfield(pout, uvfbuf,
573 greg 2.11 (n+m < 2*zp->ntotal-2) ?
574     "\n " : "\n\n")) {
575 greg 2.3 fputs(progname, stderr);
576     fputs(": error adding UVF fields\n", stderr);
577     return(0);
578     }
579     }
580 greg 2.4 if (vfsum < 0.95)
581     fprintf(stderr,
582     "%s: warning - missing %.1f%% of energy from surface '%s'\n",
583     progname, 100.*(1.-vfsum),
584     idf_getfield(pptr,NAME_FLD)->val);
585 greg 2.3 }
586     free(uvfa); /* clean up and return */
587     return(1);
588 greg 2.2 }
589    
590 greg 2.1 /* Compute zone User View Factors */
591     static int
592     compute_zones(void)
593     {
594 greg 2.2 ZONE *zptr;
595     /* temporary octree name */
596 greg 2.5 mktemp(strcpy(temp_octree, TEMPLATE));
597 greg 2.2 /* compute each zone */
598     for (zptr = zone_list; zptr != NULL; zptr = zptr->next) {
599     SUBPROC rcproc;
600     /* start rcontrib process */
601     if (!start_rcontrib(&rcproc, zptr))
602     return(0);
603     /* compute+add view factors */
604     if (!compute_uvfs(&rcproc, zptr))
605     return(0);
606     if (close_process(&rcproc) != 0) {
607     fputs(progname, stderr);
608     fputs(": bad return status from rcontrib\n", stderr);
609     return(0);
610     }
611     }
612     unlink(temp_octree); /* remove octree file */
613     return(1);
614 greg 2.1 }
615    
616     /* Load IDF and compute User View Factors */
617     int
618     main(int argc, char *argv[])
619     {
620 greg 2.2 int incl_comments = 1;
621 greg 2.1 char *origIDF, *revIDF;
622 greg 2.9 IDF_OBJECT *pptr;
623 greg 2.1 int i;
624    
625 greg 2.7 progname = *argv++; argc--; /* get options if any */
626     while (argc > 1 && argv[0][0] == '-')
627     switch (argv[0][1]) {
628     case 'c': /* elide comments */
629     incl_comments = -1; /* header only */
630     argv++; argc--;
631     continue;
632     case 's': /* samples */
633     nsamps = 1000*atoi(*++argv);
634     argv++; argc -= 2;
635     continue;
636     default:
637     fputs(progname, stderr);
638     fputs(": unknown option '", stderr);
639     fputs(argv[0], stderr);
640     fputs("'\n", stderr);
641     goto userr;
642     }
643     if ((argc < 1) | (argc > 2))
644     goto userr;
645     origIDF = argv[0];
646     revIDF = (argc == 1) ? argv[0] : argv[1];
647 greg 2.1 /* load Input Data File */
648     our_idf = idf_load(origIDF);
649     if (our_idf == NULL) {
650 greg 2.2 fputs(progname, stderr);
651 greg 2.1 fputs(": cannot load IDF '", stderr);
652     fputs(origIDF, stderr);
653     fputs("'\n", stderr);
654     return(1);
655     }
656     /* remove existing UVFs */
657 greg 2.9 if ((pptr = idf_getobject(our_idf, UVF_PNAME)) != NULL) {
658     IDF_OBJECT *pnext;
659 greg 2.2 fputs(progname, stderr);
660 greg 2.1 fputs(": removing previous User View Factors\n", stderr);
661     do {
662     pnext = pptr->pnext;
663 greg 2.9 idf_delobject(our_idf, pptr);
664 greg 2.1 } while (pnext != NULL);
665     }
666     /* add to header */
667     if (our_idf->hrem == NULL ||
668     (i = strlen(our_idf->hrem)-strlen(ADD_HEADER)) < 0 ||
669     strcmp(our_idf->hrem+i, ADD_HEADER))
670     idf_add2hdr(our_idf, ADD_HEADER);
671     /* gather zone surfaces */
672     for (i = 0; surf_type[i].pname != NULL; i++)
673 greg 2.9 for (pptr = idf_getobject(our_idf, surf_type[i].pname);
674 greg 2.1 pptr != NULL; pptr = pptr->pnext) {
675     IDF_FIELD *fptr = idf_getfield(pptr,
676     surf_type[i].zone_fld);
677     if (fptr == NULL) {
678 greg 2.2 fputs(progname, stderr);
679     fputs(": warning - missing zone field\n", stderr);
680 greg 2.1 continue;
681     }
682 greg 2.3 if (add2zone(pptr, fptr->val) == NULL)
683 greg 2.1 return(1);
684     }
685 greg 2.11 /* add subsurfaces */
686     for (pptr = idf_getobject(our_idf, SUBSURF_PNAME);
687     pptr != NULL; pptr = pptr->pnext)
688     if (add_subsurf(pptr) == NULL)
689     return(1);
690 greg 2.1 /* run rcontrib on each zone */
691     if (!compute_zones())
692     return(1);
693     /* write out modified IDF */
694 greg 2.2 if (!idf_write(our_idf, revIDF, incl_comments)) {
695     fputs(progname, stderr);
696 greg 2.1 fputs(": error writing IDF '", stderr);
697     fputs(revIDF, stderr);
698     fputs("'\n", stderr);
699     return(1);
700     }
701     return(0); /* finito! */
702 greg 2.7 userr:
703     fputs("Usage: ", stderr);
704     fputs(progname, stderr);
705     fputs(" [-c][-s Ksamps] Model.idf [Revised.idf]\n", stderr);
706     return(1);
707 greg 2.1 }