ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/eplus_adduvf.c
Revision: 2.8
Committed: Wed Feb 12 00:12:38 2014 UTC (10 years, 2 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.7: +3 -3 lines
Log Message:
Fixed error-checking

File Contents

# Content
1 #ifndef lint
2 static const char RCSid[] = "$Id: eplus_adduvf.c,v 2.7 2014/02/11 23:44:11 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 80000 /* default number of samples */
20 #endif
21
22 char *progname; /* global argv[0] */
23
24 int nsamps = NSAMPLES; /* number of samples to use */
25
26 char temp_octree[128]; /* temporary octree */
27
28 const char UVF_PNAME[] =
29 "ZoneProperty:UserViewFactor:bySurfaceName";
30
31 const char ADD_HEADER[] =
32 "\n!+++ User View Factors computed by Radiance +++!\n\n";
33
34 #define NAME_FLD 1 /* name field always first? */
35
36 typedef struct {
37 const char *pname; /* parameter type name */
38 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 {"Floor:Detailed", 3, 9},
45 {"RoofCeiling:Detailed", 3, 9},
46 {"Wall:Detailed", 3, 9},
47 {NULL}
48 }; /* IDF surface types */
49
50 typedef struct s_zone {
51 const char *zname; /* zone name */
52 struct s_zone *next; /* next zone in list */
53 int nsurf; /* surface count */
54 IDF_PARAMETER *pfirst; /* first matching parameter */
55 IDF_PARAMETER *plast; /* last matching parameter */
56 } ZONE; /* a list of collected zone surfaces */
57
58 ZONE *zone_list = NULL; /* our list of zones */
59
60 IDF_LOADED *our_idf = NULL; /* loaded/modified IDF */
61
62 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 /* Create a new zone and push to top of our list */
71 static ZONE *
72 new_zone(const char *zname, IDF_PARAMETER *param)
73 {
74 ZONE *znew = (ZONE *)malloc(sizeof(ZONE));
75
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 znew->nsurf = 1;
82 return(zone_list = znew);
83 }
84
85 /* Add the detailed surface (polygon) to the named zone */
86 static ZONE *
87 add2zone(IDF_PARAMETER *param, const char *zname)
88 {
89 ZONE *zptr;
90
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 if (!idf_movparam(our_idf, param, zptr->plast))
98 return(NULL);
99 zptr->plast = param;
100 zptr->nsurf++;
101 return(zptr);
102 }
103
104 /* Return field for vertices in the given parameter */
105 static IDF_FIELD *
106 get_vlist(IDF_PARAMETER *param, const char *zname)
107 {
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 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 /* return field for #verts */
121 return(idf_getfield(param, surf_type[i].vert_fld));
122 }
123
124 /* Convert surface to Radiance with modifier based on unique name */
125 static int
126 rad_surface(IDF_PARAMETER *param, FILE *ofp)
127 {
128 const char *sname = idf_getfield(param,NAME_FLD)->val;
129 IDF_FIELD *fptr = get_vlist(param, NULL);
130 int nvert, i;
131
132 if (fptr == NULL || (nvert = atoi(fptr->val)) < 3) {
133 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 if (fptr == NULL || !isflt(fptr->val)) {
142 fprintf(stderr,
143 "%s: missing/bad vertex for %s '%s'\n",
144 progname, param->pname, sname);
145 return(0);
146 }
147 fputc('\t', ofp);
148 fputs(fptr->val, ofp);
149 }
150 fputc('\n', ofp);
151 }
152 return(!ferror(ofp));
153 }
154
155 /* Start rcontrib process */
156 static int
157 start_rcontrib(SUBPROC *pd, ZONE *zp)
158 {
159 #define BASE_AC 5
160 static char *base_av[BASE_AC] = {
161 "rcontrib", "-fff", "-h", "-x", "1"
162 };
163 char cbuf[300];
164 char **av;
165 FILE *ofp;
166 IDF_PARAMETER *pptr;
167 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 sprintf(cbuf, "%d", nsamps);
182 av[i++] = "-c";
183 av[i++] = cbuf; /* add modifier arguments */
184 for (n = zp->nsurf, pptr = zp->pfirst; n--; pptr = pptr->dnext) {
185 IDF_FIELD *fptr = idf_getfield(pptr,NAME_FLD);
186 if (fptr == NULL || !fptr->val[0]) {
187 fputs(progname, stderr);
188 fputs(": missing name for surface parameter\n", stderr);
189 return(0);
190 }
191 if (!rad_surface(pptr, ofp)) /* add surface to octree */
192 return(0);
193 av[i++] = "-m";
194 av[i++] = fptr->val;
195 }
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 if (open_process(pd, av) <= 0) { /* start process */
204 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 /* 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 for (i = nv; i--; ) /* reverse vertex ordering */
230 for (j = 0; j < 3; j++) {
231 if (fptr == NULL) {
232 fputs(progname, stderr);
233 fputs(": missing vertex in init_poly()\n", stderr);
234 return(NULL);
235 }
236 vl3[i][j] = atof(fptr->val);
237 fptr = fptr->next;
238 }
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 ps->poff = DOT(vl3[0], ps->sdir[2]);
262 /* 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 (ps->poff+.001)*ps->sdir[2][i];
286 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 for (i = ns; i--; ) { /* stratified Monte Carlo sampling */
315 double sv[4];
316 FVECT dv;
317 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 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 }
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 sample_surface(IDF_PARAMETER *param, int wd)
341 {
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 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 progname, param->pname,
351 idf_getfield(param,NAME_FLD)->val);
352 return(0);
353 }
354 psamp.samp_left = nsamps; /* assign samples & destination */
355 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 /* Compute User View Factors using open rcontrib process */
364 static int
365 compute_uvfs(SUBPROC *pd, ZONE *zp)
366 {
367 IDF_PARAMETER *pptr, *pout, *pptr1;
368 float *uvfa;
369 char uvfbuf[24];
370 int n, m;
371 /* create output parameter */
372 pout = idf_newparam(our_idf, UVF_PNAME,
373 " ! computed by Radiance\n ", zp->plast);
374 if (pout == NULL) {
375 fputs(progname, stderr);
376 fputs(": cannot create new IDF parameter\n", stderr);
377 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 double vfsum = 0;
392 /* 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 vfsum += uvfa[3*m + 1];
406 if (pptr1 == pptr) {
407 if (uvfa[3*m + 1] > .001)
408 fprintf(stderr,
409 "%s: warning - non-zero self-VF (%.1f%%) for surface '%s'\n",
410 progname, 100.*uvfa[3*m + 1],
411 idf_getfield(pptr,NAME_FLD)->val);
412 continue; /* don't record self-factor */
413 }
414 sprintf(uvfbuf, "%.4f", uvfa[3*m + 1]);
415 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 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 }
433 free(uvfa); /* clean up and return */
434 return(1);
435 }
436
437 /* Compute zone User View Factors */
438 static int
439 compute_zones(void)
440 {
441 ZONE *zptr;
442 /* temporary octree name */
443 mktemp(strcpy(temp_octree, TEMPLATE));
444 /* 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 }
462
463 /* Load IDF and compute User View Factors */
464 int
465 main(int argc, char *argv[])
466 {
467 int incl_comments = 1;
468 char *origIDF, *revIDF;
469 IDF_PARAMETER *pptr;
470 int i;
471
472 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 /* load Input Data File */
495 our_idf = idf_load(origIDF);
496 if (our_idf == NULL) {
497 fputs(progname, stderr);
498 fputs(": cannot load IDF '", stderr);
499 fputs(origIDF, stderr);
500 fputs("'\n", stderr);
501 return(1);
502 }
503 /* remove existing UVFs */
504 if ((pptr = idf_getparam(our_idf, UVF_PNAME)) != NULL) {
505 IDF_PARAMETER *pnext;
506 fputs(progname, stderr);
507 fputs(": removing previous User View Factors\n", stderr);
508 do {
509 pnext = pptr->pnext;
510 idf_delparam(our_idf, pptr);
511 } 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 for (pptr = idf_getparam(our_idf, surf_type[i].pname);
521 pptr != NULL; pptr = pptr->pnext) {
522 IDF_FIELD *fptr = idf_getfield(pptr,
523 surf_type[i].zone_fld);
524 if (fptr == NULL) {
525 fputs(progname, stderr);
526 fputs(": warning - missing zone field\n", stderr);
527 continue;
528 }
529 if (add2zone(pptr, fptr->val) == NULL)
530 return(1);
531 }
532 /* run rcontrib on each zone */
533 if (!compute_zones())
534 return(1);
535 /* write out modified IDF */
536 if (!idf_write(our_idf, revIDF, incl_comments)) {
537 fputs(progname, stderr);
538 fputs(": error writing IDF '", stderr);
539 fputs(revIDF, stderr);
540 fputs("'\n", stderr);
541 return(1);
542 }
543 return(0); /* finito! */
544 userr:
545 fputs("Usage: ", stderr);
546 fputs(progname, stderr);
547 fputs(" [-c][-s Ksamps] Model.idf [Revised.idf]\n", stderr);
548 return(1);
549 }