ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/eplus_adduvf.c
Revision: 2.3
Committed: Mon Feb 10 04:51:26 2014 UTC (10 years, 2 months ago) by greg
Content type: text/plain
Branch: MAIN
Changes since 2.2: +246 -32 lines
Log Message:
Complete but untested implmentation of eplus_adduvf

File Contents

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