ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/eplus_adduvf.c
Revision: 2.17
Committed: Fri Apr 4 17:56:45 2014 UTC (9 years, 11 months ago) by greg
Content type: text/plain
Branch: MAIN
CVS Tags: rad4R2P2, rad5R0, rad4R2, rad4R2P1
Changes since 2.16: +3 -3 lines
Log Message:
Fixed major bug in subsurface calculation

File Contents

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