ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/eplus_adduvf.c
Revision: 2.5
Committed: Tue Feb 11 21:17:29 2014 UTC (10 years, 2 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.4: +31 -28 lines
Log Message:
Initial working version -- needs man page and accuracy testing

File Contents

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