ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/eplus_adduvf.c
Revision: 2.9
Committed: Wed Feb 12 17:40:07 2014 UTC (10 years, 1 month ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.8: +22 -22 lines
Log Message:
Changed all "parameter" references to "object"

File Contents

# User Rev Content
1 greg 2.1 #ifndef lint
2 greg 2.9 static const char RCSid[] = "$Id: eplus_adduvf.c,v 2.8 2014/02/12 00:12:38 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.2 char *progname; /* global argv[0] */
23    
24 greg 2.7 int nsamps = NSAMPLES; /* number of samples to use */
25    
26 greg 2.5 char temp_octree[128]; /* temporary octree */
27 greg 2.2
28     const char UVF_PNAME[] =
29     "ZoneProperty:UserViewFactor:bySurfaceName";
30 greg 2.1
31     const char ADD_HEADER[] =
32 greg 2.5 "\n!+++ User View Factors computed by Radiance +++!\n\n";
33 greg 2.1
34     #define NAME_FLD 1 /* name field always first? */
35    
36     typedef struct {
37 greg 2.9 const char *pname; /* object type name */
38 greg 2.1 short zone_fld; /* zone field index */
39     short vert_fld; /* vertex field index */
40     } SURF_PTYPE; /* surface type we're interested in */
41    
42     const SURF_PTYPE surf_type[] = {
43     {"BuildingSurface:Detailed", 4, 10},
44 greg 2.5 {"Floor:Detailed", 3, 9},
45     {"RoofCeiling:Detailed", 3, 9},
46     {"Wall:Detailed", 3, 9},
47 greg 2.1 {NULL}
48 greg 2.7 }; /* IDF surface types */
49 greg 2.1
50 greg 2.2 typedef struct s_zone {
51 greg 2.1 const char *zname; /* zone name */
52     struct s_zone *next; /* next zone in list */
53 greg 2.2 int nsurf; /* surface count */
54 greg 2.9 IDF_OBJECT *pfirst; /* first matching object */
55     IDF_OBJECT *plast; /* last matching object */
56 greg 2.2 } ZONE; /* a list of collected zone surfaces */
57    
58     ZONE *zone_list = NULL; /* our list of zones */
59 greg 2.1
60     IDF_LOADED *our_idf = NULL; /* loaded/modified IDF */
61    
62 greg 2.7 typedef struct {
63     FVECT sdir[3]; /* UVW unit sampling vectors */
64     double poff; /* W-offset for plane of polygon */
65     double area_left; /* area left to sample */
66     int samp_left; /* remaining samples */
67     int wd; /* output file descriptor */
68     } POLYSAMP; /* structure for polygon sampling */
69    
70 greg 2.1 /* Create a new zone and push to top of our list */
71 greg 2.2 static ZONE *
72 greg 2.9 new_zone(const char *zname, IDF_OBJECT *param)
73 greg 2.1 {
74 greg 2.2 ZONE *znew = (ZONE *)malloc(sizeof(ZONE));
75 greg 2.1
76     if (znew == NULL)
77     return(NULL);
78     znew->zname = zname; /* assumes static */
79     znew->next = zone_list;
80     znew->pfirst = znew->plast = param;
81 greg 2.2 znew->nsurf = 1;
82 greg 2.1 return(zone_list = znew);
83     }
84    
85     /* Add the detailed surface (polygon) to the named zone */
86 greg 2.2 static ZONE *
87 greg 2.9 add2zone(IDF_OBJECT *param, const char *zname)
88 greg 2.1 {
89 greg 2.2 ZONE *zptr;
90 greg 2.1
91     for (zptr = zone_list; zptr != NULL; zptr = zptr->next)
92     if (!strcmp(zptr->zname, zname))
93     break;
94     if (zptr == NULL)
95     return(new_zone(zname, param));
96     /* keep surfaces together */
97 greg 2.9 if (!idf_movobject(our_idf, param, zptr->plast))
98 greg 2.1 return(NULL);
99     zptr->plast = param;
100 greg 2.2 zptr->nsurf++;
101 greg 2.1 return(zptr);
102     }
103    
104 greg 2.9 /* Return field for vertices in the given object */
105 greg 2.3 static IDF_FIELD *
106 greg 2.9 get_vlist(IDF_OBJECT *param, const char *zname)
107 greg 2.2 {
108     int i = 0;
109     IDF_FIELD *fptr;
110     /* check for surface type */
111     while (strcmp(surf_type[i].pname, param->pname))
112     if (surf_type[++i].pname == NULL)
113 greg 2.3 return(NULL);
114    
115     if (zname != NULL) { /* matches specified zone? */
116     fptr = idf_getfield(param, surf_type[i].zone_fld);
117     if (fptr == NULL || strcmp(fptr->val, zname))
118     return(NULL);
119     }
120 greg 2.2 /* return field for #verts */
121 greg 2.3 return(idf_getfield(param, surf_type[i].vert_fld));
122 greg 2.2 }
123    
124     /* Convert surface to Radiance with modifier based on unique name */
125     static int
126 greg 2.9 rad_surface(IDF_OBJECT *param, FILE *ofp)
127 greg 2.2 {
128 greg 2.3 const char *sname = idf_getfield(param,NAME_FLD)->val;
129     IDF_FIELD *fptr = get_vlist(param, NULL);
130 greg 2.2 int nvert, i;
131    
132 greg 2.3 if (fptr == NULL || (nvert = atoi(fptr->val)) < 3) {
133 greg 2.2 fprintf(stderr, "%s: bad surface '%s'\n", progname, sname);
134     return(0);
135     }
136     fprintf(ofp, "\nvoid glow '%s'\n0\n0\n4 1 1 1 0\n", sname);
137     fprintf(ofp, "\n'%s' polygon 's_%s'\n0\n0\n%d\n", sname, sname, 3*nvert);
138     while (nvert--) {
139     for (i = 3; i--; ) {
140     fptr = fptr->next;
141 greg 2.3 if (fptr == NULL || !isflt(fptr->val)) {
142 greg 2.2 fprintf(stderr,
143 greg 2.3 "%s: missing/bad vertex for %s '%s'\n",
144     progname, param->pname, sname);
145 greg 2.2 return(0);
146     }
147     fputc('\t', ofp);
148 greg 2.3 fputs(fptr->val, ofp);
149 greg 2.2 }
150     fputc('\n', ofp);
151     }
152     return(!ferror(ofp));
153     }
154    
155 greg 2.3 /* Start rcontrib process */
156 greg 2.2 static int
157     start_rcontrib(SUBPROC *pd, ZONE *zp)
158     {
159 greg 2.5 #define BASE_AC 5
160 greg 2.2 static char *base_av[BASE_AC] = {
161 greg 2.8 "rcontrib", "-fff", "-h", "-x", "1"
162 greg 2.2 };
163     char cbuf[300];
164     char **av;
165     FILE *ofp;
166 greg 2.9 IDF_OBJECT *pptr;
167 greg 2.2 int i, n;
168     /* start oconv command */
169     sprintf(cbuf, "oconv - > '%s'", temp_octree);
170     if ((ofp = popen(cbuf, "w")) == NULL) {
171     fputs(progname, stderr);
172     fputs(": cannot open oconv process\n", stderr);
173     return(0);
174     }
175     /* allocate argument list */
176     av = (char **)malloc(sizeof(char *)*(BASE_AC+4+2*zp->nsurf));
177     if (av == NULL)
178     return(0);
179     for (i = 0; i < BASE_AC; i++)
180     av[i] = base_av[i];
181 greg 2.7 sprintf(cbuf, "%d", nsamps);
182 greg 2.2 av[i++] = "-c";
183     av[i++] = cbuf; /* add modifier arguments */
184     for (n = zp->nsurf, pptr = zp->pfirst; n--; pptr = pptr->dnext) {
185 greg 2.3 IDF_FIELD *fptr = idf_getfield(pptr,NAME_FLD);
186     if (fptr == NULL || !fptr->val[0]) {
187 greg 2.2 fputs(progname, stderr);
188 greg 2.9 fputs(": missing name for surface object\n", stderr);
189 greg 2.2 return(0);
190     }
191 greg 2.3 if (!rad_surface(pptr, ofp)) /* add surface to octree */
192     return(0);
193 greg 2.2 av[i++] = "-m";
194 greg 2.3 av[i++] = fptr->val;
195 greg 2.2 }
196     if (pclose(ofp) != 0) { /* finish oconv */
197     fputs(progname, stderr);
198     fputs(": error running oconv process\n", stderr);
199     return(0);
200     }
201     av[i++] = temp_octree; /* add final octree argument */
202     av[i] = NULL;
203 greg 2.8 if (open_process(pd, av) <= 0) { /* start process */
204 greg 2.2 fputs(progname, stderr);
205     fputs(": cannot start rcontrib process\n", stderr);
206     return(0);
207     }
208     free(av); /* all done -- clean up */
209     return(1);
210     #undef BASE_AC
211     }
212    
213 greg 2.3 /* Initialize polygon sampling */
214     static Vert2_list *
215     init_poly(POLYSAMP *ps, IDF_FIELD *f0, int nv)
216     {
217     IDF_FIELD *fptr = f0;
218     int i, j;
219     FVECT *vl3, e1, e2, vc;
220     Vert2_list *vl2 = polyAlloc(nv);
221    
222     if (vl2 == NULL)
223     return(NULL);
224     vl2->p = ps;
225     /* get 3-D vertices */
226     vl3 = (FVECT *)malloc(sizeof(FVECT)*nv);
227     if (vl3 == NULL)
228     return(NULL);
229 greg 2.5 for (i = nv; i--; ) /* reverse vertex ordering */
230 greg 2.3 for (j = 0; j < 3; j++) {
231     if (fptr == NULL) {
232     fputs(progname, stderr);
233 greg 2.5 fputs(": missing vertex in init_poly()\n", stderr);
234 greg 2.3 return(NULL);
235     }
236     vl3[i][j] = atof(fptr->val);
237 greg 2.5 fptr = fptr->next;
238 greg 2.3 }
239     /* compute area and normal */
240     ps->sdir[2][0] = ps->sdir[2][1] = ps->sdir[2][2] = 0;
241     VSUB(e1, vl3[1], vl3[0]);
242     for (i = 2; i < nv; i++) {
243     VSUB(e2, vl3[i], vl3[0]);
244     fcross(vc, e1, e2);
245     ps->sdir[2][0] += vc[0];
246     ps->sdir[2][1] += vc[1];
247     ps->sdir[2][2] += vc[2];
248     VCOPY(e1, e2);
249     }
250     ps->area_left = .5 * normalize(ps->sdir[2]);
251     if (ps->area_left == .0) {
252     fputs(progname, stderr);
253     fputs(": degenerate polygon in init_poly()\n", stderr);
254     return(0);
255     }
256     /* create X & Y axes */
257     VCOPY(ps->sdir[0], e1);
258     normalize(ps->sdir[0]);
259     fcross(ps->sdir[1], ps->sdir[2], ps->sdir[0]);
260     /* compute plane offset */
261 greg 2.5 ps->poff = DOT(vl3[0], ps->sdir[2]);
262 greg 2.3 /* assign 2-D vertices */
263     for (i = 0; i < nv; i++) {
264     vl2->v[i].mX = DOT(vl3[i], ps->sdir[0]);
265     vl2->v[i].mY = DOT(vl3[i], ps->sdir[1]);
266     }
267     free(vl3); /* it's ready! */
268     return(vl2);
269     }
270    
271     /* Generate samples on 2-D triangle */
272     static int
273     sample_triangle(const Vert2_list *vl2, int a, int b, int c)
274     {
275     POLYSAMP *ps = (POLYSAMP *)vl2->p;
276     float *samp;
277     FVECT orig;
278     FVECT ab, ac;
279     double area;
280     int i, j, ns;
281     /* compute sampling axes */
282     for (i = 3; i--; ) {
283     orig[i] = vl2->v[a].mX*ps->sdir[0][i] +
284     vl2->v[a].mY*ps->sdir[1][i] +
285 greg 2.5 (ps->poff+.001)*ps->sdir[2][i];
286 greg 2.3 ab[i] = (vl2->v[b].mX - vl2->v[a].mX)*ps->sdir[0][i] +
287     (vl2->v[b].mY - vl2->v[a].mY)*ps->sdir[1][i];
288     ac[i] = (vl2->v[c].mX - vl2->v[a].mX)*ps->sdir[0][i] +
289     (vl2->v[c].mY - vl2->v[a].mY)*ps->sdir[1][i];
290     }
291     /* compute number of samples to take */
292     area = .5*(vl2->v[a].mX*vl2->v[b].mY - vl2->v[b].mX*vl2->v[a].mY +
293     vl2->v[b].mX*vl2->v[c].mY - vl2->v[c].mX*vl2->v[b].mY +
294     vl2->v[c].mX*vl2->v[a].mY - vl2->v[a].mX*vl2->v[c].mY);
295     if (area < .0) {
296     fputs(progname, stderr);
297     fputs(": negative triangle area in sample_triangle()\n", stderr);
298     return(0);
299     }
300     if (area >= ps->area_left) {
301     ns = ps->samp_left;
302     ps->area_left = 0;
303     } else {
304     ns = (ps->samp_left*area/ps->area_left + .5);
305     ps->samp_left -= ns;
306     ps->area_left -= area;
307     }
308     if (ns <= 0) /* XXX should be error? */
309     return(1);
310     /* buffer sample rays */
311     samp = (float *)malloc(sizeof(float)*6*ns);
312     if (samp == NULL)
313     return(0);
314 greg 2.4 for (i = ns; i--; ) { /* stratified Monte Carlo sampling */
315 greg 2.3 double sv[4];
316 greg 2.5 FVECT dv;
317 greg 2.3 multisamp(sv, 4, (i+frandom())/(double)ns);
318     sv[0] *= sv[1] = sqrt(sv[1]);
319     sv[1] = 1. - sv[1];
320     for (j = 3; j--; )
321 greg 2.5 samp[i*6 + j] = orig[j] + sv[0]*ab[j] + sv[1]*ac[j];
322     sv[2] = sqrt(sv[2]);
323     sv[3] *= 2.*PI;
324     dv[0] = tcos(sv[3]) * sv[2];
325     dv[1] = tsin(sv[3]) * sv[2];
326     dv[2] = sqrt(1. - sv[2]*sv[2]);
327     for (j = 3; j--; )
328     samp[i*6 + 3 + j] = dv[0]*ps->sdir[0][j] +
329     dv[1]*ps->sdir[1][j] +
330     dv[2]*ps->sdir[2][j] ;
331 greg 2.3 }
332     /* send to our process */
333     writebuf(ps->wd, (char *)samp, sizeof(float)*6*ns);
334     free(samp); /* that's it! */
335     return(1);
336     }
337    
338     /* Sample the given surface */
339     static int
340 greg 2.9 sample_surface(IDF_OBJECT *param, int wd)
341 greg 2.3 {
342     IDF_FIELD *fptr = get_vlist(param, NULL);
343     POLYSAMP psamp;
344     int nv;
345     Vert2_list *vlist2;
346     /* set up our polygon sampler */
347 greg 2.5 if (fptr == NULL || (nv = atoi(fptr->val)) < 3 ||
348     (vlist2 = init_poly(&psamp, fptr->next, nv)) == NULL) {
349     fprintf(stderr, "%s: bad polygon %s '%s'\n",
350 greg 2.3 progname, param->pname,
351     idf_getfield(param,NAME_FLD)->val);
352     return(0);
353     }
354 greg 2.7 psamp.samp_left = nsamps; /* assign samples & destination */
355 greg 2.3 psamp.wd = wd;
356     /* sample each subtriangle */
357     if (!polyTriangulate(vlist2, &sample_triangle))
358     return(0);
359     polyFree(vlist2); /* clean up and return */
360     return(1);
361     }
362    
363 greg 2.2 /* Compute User View Factors using open rcontrib process */
364     static int
365 greg 2.3 compute_uvfs(SUBPROC *pd, ZONE *zp)
366 greg 2.2 {
367 greg 2.9 IDF_OBJECT *pptr, *pout, *pptr1;
368 greg 2.3 float *uvfa;
369     char uvfbuf[24];
370     int n, m;
371 greg 2.9 /* create output object */
372     pout = idf_newobject(our_idf, UVF_PNAME,
373 greg 2.3 " ! computed by Radiance\n ", zp->plast);
374     if (pout == NULL) {
375     fputs(progname, stderr);
376 greg 2.9 fputs(": cannot create new IDF object\n", stderr);
377 greg 2.3 return(0);
378     }
379     if (!idf_addfield(pout, zp->zname,
380     " !- Zone Name\n ")) {
381     fputs(progname, stderr);
382     fputs(": cannot add zone name field\n", stderr);
383     return(0);
384     }
385     /* allocate read buffer */
386     uvfa = (float *)malloc(sizeof(float)*3*zp->nsurf);
387     if (uvfa == NULL)
388     return(0);
389     /* UVFs from each surface */
390     for (n = zp->nsurf, pptr = zp->pfirst; n--; pptr = pptr->dnext) {
391 greg 2.4 double vfsum = 0;
392 greg 2.3 /* send samples to rcontrib */
393     if (!sample_surface(pptr, pd->w))
394     return(0);
395     /* read results */
396     if (readbuf(pd->r, (char *)uvfa, sizeof(float)*3*zp->nsurf) !=
397     sizeof(float)*3*zp->nsurf) {
398     fputs(progname, stderr);
399     fputs(": read error from rcontrib process\n", stderr);
400     return(0);
401     }
402     /* append UVF fields */
403     for (m = 0, pptr1 = zp->pfirst;
404     m < zp->nsurf; m++, pptr1 = pptr1->dnext) {
405 greg 2.4 vfsum += uvfa[3*m + 1];
406 greg 2.3 if (pptr1 == pptr) {
407     if (uvfa[3*m + 1] > .001)
408     fprintf(stderr,
409 greg 2.5 "%s: warning - non-zero self-VF (%.1f%%) for surface '%s'\n",
410     progname, 100.*uvfa[3*m + 1],
411 greg 2.3 idf_getfield(pptr,NAME_FLD)->val);
412     continue; /* don't record self-factor */
413     }
414 greg 2.6 sprintf(uvfbuf, "%.4f", uvfa[3*m + 1]);
415 greg 2.3 if (!idf_addfield(pout,
416     idf_getfield(pptr,NAME_FLD)->val, NULL) ||
417     !idf_addfield(pout,
418     idf_getfield(pptr1,NAME_FLD)->val, NULL) ||
419     !idf_addfield(pout, uvfbuf,
420     (n || m < zp->nsurf-2) ?
421     "\n " : "\n\n")) {
422     fputs(progname, stderr);
423     fputs(": error adding UVF fields\n", stderr);
424     return(0);
425     }
426     }
427 greg 2.4 if (vfsum < 0.95)
428     fprintf(stderr,
429     "%s: warning - missing %.1f%% of energy from surface '%s'\n",
430     progname, 100.*(1.-vfsum),
431     idf_getfield(pptr,NAME_FLD)->val);
432 greg 2.3 }
433     free(uvfa); /* clean up and return */
434     return(1);
435 greg 2.2 }
436    
437 greg 2.1 /* Compute zone User View Factors */
438     static int
439     compute_zones(void)
440     {
441 greg 2.2 ZONE *zptr;
442     /* temporary octree name */
443 greg 2.5 mktemp(strcpy(temp_octree, TEMPLATE));
444 greg 2.2 /* compute each zone */
445     for (zptr = zone_list; zptr != NULL; zptr = zptr->next) {
446     SUBPROC rcproc;
447     /* start rcontrib process */
448     if (!start_rcontrib(&rcproc, zptr))
449     return(0);
450     /* compute+add view factors */
451     if (!compute_uvfs(&rcproc, zptr))
452     return(0);
453     if (close_process(&rcproc) != 0) {
454     fputs(progname, stderr);
455     fputs(": bad return status from rcontrib\n", stderr);
456     return(0);
457     }
458     }
459     unlink(temp_octree); /* remove octree file */
460     return(1);
461 greg 2.1 }
462    
463     /* Load IDF and compute User View Factors */
464     int
465     main(int argc, char *argv[])
466     {
467 greg 2.2 int incl_comments = 1;
468 greg 2.1 char *origIDF, *revIDF;
469 greg 2.9 IDF_OBJECT *pptr;
470 greg 2.1 int i;
471    
472 greg 2.7 progname = *argv++; argc--; /* get options if any */
473     while (argc > 1 && argv[0][0] == '-')
474     switch (argv[0][1]) {
475     case 'c': /* elide comments */
476     incl_comments = -1; /* header only */
477     argv++; argc--;
478     continue;
479     case 's': /* samples */
480     nsamps = 1000*atoi(*++argv);
481     argv++; argc -= 2;
482     continue;
483     default:
484     fputs(progname, stderr);
485     fputs(": unknown option '", stderr);
486     fputs(argv[0], stderr);
487     fputs("'\n", stderr);
488     goto userr;
489     }
490     if ((argc < 1) | (argc > 2))
491     goto userr;
492     origIDF = argv[0];
493     revIDF = (argc == 1) ? argv[0] : argv[1];
494 greg 2.1 /* load Input Data File */
495     our_idf = idf_load(origIDF);
496     if (our_idf == NULL) {
497 greg 2.2 fputs(progname, stderr);
498 greg 2.1 fputs(": cannot load IDF '", stderr);
499     fputs(origIDF, stderr);
500     fputs("'\n", stderr);
501     return(1);
502     }
503     /* remove existing UVFs */
504 greg 2.9 if ((pptr = idf_getobject(our_idf, UVF_PNAME)) != NULL) {
505     IDF_OBJECT *pnext;
506 greg 2.2 fputs(progname, stderr);
507 greg 2.1 fputs(": removing previous User View Factors\n", stderr);
508     do {
509     pnext = pptr->pnext;
510 greg 2.9 idf_delobject(our_idf, pptr);
511 greg 2.1 } while (pnext != NULL);
512     }
513     /* add to header */
514     if (our_idf->hrem == NULL ||
515     (i = strlen(our_idf->hrem)-strlen(ADD_HEADER)) < 0 ||
516     strcmp(our_idf->hrem+i, ADD_HEADER))
517     idf_add2hdr(our_idf, ADD_HEADER);
518     /* gather zone surfaces */
519     for (i = 0; surf_type[i].pname != NULL; i++)
520 greg 2.9 for (pptr = idf_getobject(our_idf, surf_type[i].pname);
521 greg 2.1 pptr != NULL; pptr = pptr->pnext) {
522     IDF_FIELD *fptr = idf_getfield(pptr,
523     surf_type[i].zone_fld);
524     if (fptr == NULL) {
525 greg 2.2 fputs(progname, stderr);
526     fputs(": warning - missing zone field\n", stderr);
527 greg 2.1 continue;
528     }
529 greg 2.3 if (add2zone(pptr, fptr->val) == NULL)
530 greg 2.1 return(1);
531     }
532     /* run rcontrib on each zone */
533     if (!compute_zones())
534     return(1);
535     /* write out modified IDF */
536 greg 2.2 if (!idf_write(our_idf, revIDF, incl_comments)) {
537     fputs(progname, stderr);
538 greg 2.1 fputs(": error writing IDF '", stderr);
539     fputs(revIDF, stderr);
540     fputs("'\n", stderr);
541     return(1);
542     }
543     return(0); /* finito! */
544 greg 2.7 userr:
545     fputs("Usage: ", stderr);
546     fputs(progname, stderr);
547     fputs(" [-c][-s Ksamps] Model.idf [Revised.idf]\n", stderr);
548     return(1);
549 greg 2.1 }