ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/util/eplus_adduvf.c
(Generate patch)

Comparing ray/src/util/eplus_adduvf.c (file contents):
Revision 2.1 by greg, Sun Feb 9 05:49:21 2014 UTC vs.
Revision 2.6 by greg, Tue Feb 11 22:01:44 2014 UTC

# Line 7 | Line 7 | static const char RCSid[] = "$Id$";
7   *      G.Ward for LBNL
8   */
9  
10 – #include <stdio.h>
10   #include <stdlib.h>
11 < #include <string.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 */
19 > #define NSAMPLES        80000                   /* number of samples to use */
20   #endif
21  
22 < #define UVF_PNAME       "ZoneProperty:UserViewFactor:bySurfaceName"
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";
30 >                        "\n!+++ User View Factors computed by Radiance +++!\n\n";
31  
32   #define NAME_FLD        1                       /* name field always first? */
33  
# Line 32 | Line 39 | typedef struct {
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 < struct s_zone {
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_list = NULL;
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 struct s_zone *
61 > static ZONE *
62   new_zone(const char *zname, IDF_PARAMETER *param)
63   {
64 <        struct s_zone   *znew = (struct s_zone *)malloc(sizeof(struct s_zone));
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 struct s_zone *
76 > static ZONE *
77   add2zone(IDF_PARAMETER *param, const char *zname)
78   {
79 <        struct s_zone   *zptr;
79 >        ZONE    *zptr;
80  
81          for (zptr = zone_list; zptr != NULL; zptr = zptr->next)
82                  if (!strcmp(zptr->zname, zname))
# Line 73 | Line 87 | add2zone(IDF_PARAMETER *param, const char *zname)
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, "%.4f", 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 <        fputs("compute_zones() unimplemented\n", stderr);
440 <        return(0);
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(argv[0], stderr);
478 <                fputs(" Model.idf [Revised.idf]\n", stderr);
477 >                fputs(progname, stderr);
478 >                fputs(" [-c] Model.idf [Revised.idf]\n", stderr);
479                  return(1);
480          }
481          origIDF = argv[1];
# Line 103 | Line 483 | main(int argc, char *argv[])
483                                                  /* load Input Data File */
484          our_idf = idf_load(origIDF);
485          if (our_idf == NULL) {
486 <                fputs(argv[0], stderr);
486 >                fputs(progname, stderr);
487                  fputs(": cannot load IDF '", stderr);
488                  fputs(origIDF, stderr);
489                  fputs("'\n", stderr);
# Line 112 | Line 492 | main(int argc, char *argv[])
492                                                  /* remove existing UVFs */
493          if ((pptr = idf_getparam(our_idf, UVF_PNAME)) != NULL) {
494                  IDF_PARAMETER   *pnext;
495 <                fputs(argv[0], stderr);
495 >                fputs(progname, stderr);
496                  fputs(": removing previous User View Factors\n", stderr);
497                  do {
498                          pnext = pptr->pnext;
# Line 131 | Line 511 | main(int argc, char *argv[])
511                          IDF_FIELD       *fptr = idf_getfield(pptr,
512                                                          surf_type[i].zone_fld);
513                          if (fptr == NULL) {
514 <                                fputs("Warning: missing zone field!\n", stderr);
514 >                                fputs(progname, stderr);
515 >                                fputs(": warning - missing zone field\n", stderr);
516                                  continue;
517                          }
518 <                        if (add2zone(pptr, fptr->arg) == NULL)
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, 1)) {
526 <                fputs(argv[0], stderr);
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);

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines