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.18 by schorsch, Sun Mar 6 01:13:18 2016 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 "platform.h"
12 > #include "rtio.h"
13 > #include "rtmath.h"
14 > #include "random.h"
15   #include "eplus_idf.h"
16   #include "triangulate.h"
17 + #include "rtprocess.h"
18  
19   #ifndef NSAMPLES
20 < #define NSAMPLES        100000                  /* number of samples to use */
20 > #define NSAMPLES        80000                   /* default number of samples */
21   #endif
22  
23 < #define UVF_PNAME       "ZoneProperty:UserViewFactor:bySurfaceName"
23 > #define SURF_EPS        0.0005                  /* surface testing epsilon */
24  
25 + char            *progname;                      /* global argv[0] */
26 +
27 + int             nsamps = NSAMPLES;              /* number of samples to use */
28 +
29 + char            temp_octree[128];               /* temporary octree */
30 +
31 + const char      UVF_PNAME[] =
32 +                        "ZoneProperty:UserViewFactors:bySurfaceName";
33 +
34 + const char      SUBSURF_PNAME[] =
35 +                        "FenestrationSurface:Detailed";
36 +
37   const char      ADD_HEADER[] =
38 <                "!+++ User View Factors computed by Radiance +++!\n";
38 >                        "\n!+++ User View Factors computed by Radiance +++!\n\n";
39  
40   #define NAME_FLD        1                       /* name field always first? */
41  
42 + #define SS_BASE_FLD     4                       /* subsurface base surface */
43 + #define SS_VERT_FLD     10                      /* subsurface vertex count */
44 +
45   typedef struct {
46 <        const char      *pname;                 /* parameter type name */
46 >        const char      *pname;                 /* object type name */
47          short           zone_fld;               /* zone field index */
48          short           vert_fld;               /* vertex field index */
49   } SURF_PTYPE;           /* surface type we're interested in */
50  
51   const SURF_PTYPE        surf_type[] = {
52                  {"BuildingSurface:Detailed", 4, 10},
53 +                {"Floor:Detailed", 3, 9},
54 +                {"RoofCeiling:Detailed", 3, 9},
55 +                {"Wall:Detailed", 3, 9},
56                  {NULL}
57 <        };
57 >        };                              /* IDF surface types */
58  
59 < struct s_zone {
59 > typedef struct s_zone {
60          const char      *zname;                 /* zone name */
61          struct s_zone   *next;                  /* next zone in list */
62 <        IDF_PARAMETER   *pfirst;                /* first matching parameter */
63 <        IDF_PARAMETER   *plast;                 /* last matching parameter */
64 < } *zone_list = NULL;
62 >        int             nsurf;                  /* surface count */
63 >        int             ntotal;                 /* surfaces+subsurfaces */
64 >        IDF_OBJECT      *pfirst;                /* first matching object */
65 >        IDF_OBJECT      *plast;                 /* last before subsurfaces */
66 >        float           *area_redu;             /* subsurface area per surf. */
67 > } ZONE;                 /* a list of collected zone surfaces */
68  
69 + ZONE            *zone_list = NULL;      /* our list of zones */
70 +
71 + LUTAB           zonesurf_lut =          /* zone surface lookup table */
72 +                        LU_SINIT(NULL,NULL);
73 +
74   IDF_LOADED      *our_idf = NULL;        /* loaded/modified IDF */
75  
76 + typedef struct {
77 +        FVECT           norm;           /* surface normal */
78 +        double          area;           /* surface area */
79 +        int             nv;             /* number of vertices */
80 +        FVECT           vl[3];          /* vertex list (extends struct) */
81 + } SURFACE;              /* a polygonal surface */
82 +
83 + typedef struct {
84 +        FVECT           sdir[3];        /* UVW unit sampling vectors */
85 +        double          poff;           /* W-offset for plane of polygon */
86 +        double          area_left;      /* area left to sample */
87 +        int             samp_left;      /* remaining samples */
88 +        int             wd;             /* output file descriptor */
89 + } POLYSAMP;             /* structure for polygon sampling */
90 +
91   /* Create a new zone and push to top of our list */
92 < static struct s_zone *
93 < new_zone(const char *zname, IDF_PARAMETER *param)
92 > static ZONE *
93 > new_zone(const char *zname, IDF_OBJECT *param)
94   {
95 <        struct s_zone   *znew = (struct s_zone *)malloc(sizeof(struct s_zone));
95 >        ZONE    *znew = (ZONE *)malloc(sizeof(ZONE));
96  
97          if (znew == NULL)
98                  return(NULL);
99          znew->zname = zname;                    /* assumes static */
100          znew->next = zone_list;
101          znew->pfirst = znew->plast = param;
102 +        znew->ntotal = znew->nsurf = 1;
103 +        znew->area_redu = NULL;
104          return(zone_list = znew);
105   }
106  
107   /* Add the detailed surface (polygon) to the named zone */
108 < static struct s_zone *
109 < add2zone(IDF_PARAMETER *param, const char *zname)
108 > static ZONE *
109 > add2zone(IDF_OBJECT *param, const char *zname)
110   {
111 <        struct s_zone   *zptr;
111 >        IDF_FIELD       *nfp = idf_getfield(param, NAME_FLD);
112 >        ZONE            *zptr;
113 >        LUENT           *lep;
114  
115 +        if (nfp == NULL) {
116 +                fputs(progname, stderr);
117 +                fputs(": surface missing name field!\n", stderr);
118 +                return(NULL);
119 +        }
120          for (zptr = zone_list; zptr != NULL; zptr = zptr->next)
121                  if (!strcmp(zptr->zname, zname))
122                          break;
123 <        if (zptr == NULL)
124 <                return(new_zone(zname, param));
125 <                                                /* keep surfaces together */
126 <        if (!idf_movparam(our_idf, param, zptr->plast))
123 >        if (zptr == NULL) {
124 >                zptr = new_zone(zname, param);
125 >        } else {                                /* keep surfaces together */
126 >                if (!idf_movobject(our_idf, param, zptr->plast))
127 >                        return(NULL);
128 >                zptr->plast = param;
129 >                zptr->nsurf++;
130 >                zptr->ntotal++;
131 >        }
132 >                                                /* add to lookup table */
133 >        lep = lu_find(&zonesurf_lut, nfp->val);
134 >        if (lep == NULL)
135                  return(NULL);
136 <        zptr->plast = param;
136 >        if (lep->data != NULL) {
137 >                fputs(progname, stderr);
138 >                fputs(": duplicate surface name '", stderr);
139 >                fputs(nfp->val, stderr);
140 >                fputs("'\n", stderr);
141 >                return(NULL);
142 >        }
143 >        lep->key = nfp->val;
144 >        lep->data = (char *)zptr;
145          return(zptr);
146   }
147  
148 + /* Add a subsurface by finding its base surface and the corresponding zone */
149 + static ZONE *
150 + add_subsurf(IDF_OBJECT *param)
151 + {
152 +        IDF_FIELD       *bfp = idf_getfield(param, SS_BASE_FLD);
153 +        ZONE            *zptr;
154 +        LUENT           *lep;
155 +
156 +        if (bfp == NULL) {
157 +                fputs(progname, stderr);
158 +                fputs(": missing base field name in subsurface!\n", stderr);
159 +                return(NULL);
160 +        }
161 +        lep = lu_find(&zonesurf_lut, bfp->val); /* find base zone */
162 +        if (lep == NULL || lep->data == NULL) {
163 +                fputs(progname, stderr);
164 +                fputs(": cannot find referenced base surface '", stderr);
165 +                fputs(bfp->val, stderr);
166 +                fputs("'\n", stderr);
167 +                return(NULL);
168 +        }
169 +        zptr = (ZONE *)lep->data;               /* add this subsurface */
170 +        if (!idf_movobject(our_idf, param, zptr->plast))
171 +                return(NULL);
172 +        zptr->ntotal++;
173 +        return(zptr);
174 + }
175 +
176 + /* Return field for vertices in the given object */
177 + static IDF_FIELD *
178 + get_vlist(IDF_OBJECT *param, const char *zname)
179 + {
180 + #define itm_len         (sizeof(IDF_FIELD)+6)
181 +        static char             fld_buf[4*itm_len];
182 +        static char             *next_fbp = fld_buf;
183 +        int                     i;
184 +        IDF_FIELD               *res;
185 +                                                /* check if subsurface */
186 +        if (!strcmp(param->pname, SUBSURF_PNAME)) {
187 +                if (zname != NULL) {
188 +                        LUENT   *lep = lu_find(&zonesurf_lut,
189 +                                        idf_getfield(param,SS_BASE_FLD)->val);
190 +                        if (strcmp((*(ZONE *)lep->data).zname, zname))
191 +                                return(NULL);
192 +                }
193 +                res = idf_getfield(param, SS_VERT_FLD);
194 +        } else {
195 +                i = 0;                          /* check for surface type */
196 +                while (strcmp(surf_type[i].pname, param->pname))
197 +                        if (surf_type[++i].pname == NULL)
198 +                                return(NULL);
199 +
200 +                if (zname != NULL) {            /* matches specified zone? */
201 +                        IDF_FIELD       *fptr = idf_getfield(param, surf_type[i].zone_fld);
202 +                        if (fptr == NULL || strcmp(fptr->val, zname))
203 +                                return(NULL);
204 +                }
205 +                res = idf_getfield(param, surf_type[i].vert_fld);
206 +        }
207 +        if (!res->val[0]) {                     /* hack for missing #vert */
208 +                IDF_FIELD       *fptr;
209 +                if (next_fbp >= fld_buf+sizeof(fld_buf))
210 +                        next_fbp = fld_buf;
211 +                i = 0;                          /* count vertices */
212 +                for (fptr = res->next; fptr != NULL; fptr = fptr->next)
213 +                        ++i;
214 +                if (i % 3)
215 +                        return(NULL);
216 +                fptr = res->next;
217 +                res = (IDF_FIELD *)next_fbp; next_fbp += itm_len;
218 +                res->next = fptr;
219 +                res->rem = "";
220 +                sprintf(res->val, "%d", i/3);
221 +        }
222 +        return(res);
223 + #undef itm_len
224 + }
225 +
226 + /* Get/allocate surface polygon */
227 + static SURFACE *
228 + get_surface(IDF_FIELD *fptr)
229 + {
230 +        SURFACE *surf;
231 +        int     nv;
232 +        FVECT   e1, e2, vc;
233 +        int     i, j;
234 +                                        /* get number of vertices */
235 +        if (fptr == NULL || (nv = atoi(fptr->val)) < 3) {
236 +                fputs(progname, stderr);
237 +                fputs(": bad vertex count for surface!\n", stderr);
238 +                return(NULL);
239 +        }
240 +        surf = (SURFACE *)malloc(sizeof(SURFACE)+sizeof(FVECT)*(nv-3));
241 +        if (surf == NULL)
242 +                return(NULL);
243 +        surf->nv = nv;
244 +        for (i = nv; i--; )             /* reverse vertex order/orientation */
245 +                for (j = 0; j < 3; j++) {
246 +                        fptr = fptr->next;
247 +                        if (fptr == NULL) {
248 +                                fputs(progname, stderr);
249 +                                fputs(": missing vertex in get_surface()\n", stderr);
250 +                                free(surf);
251 +                                return(NULL);
252 +                        }
253 +                        if ((surf->vl[i][j] = atof(fptr->val)) == 0 &&
254 +                                                        !isflt(fptr->val)) {
255 +                                fputs(progname, stderr);
256 +                                fputs(": bad vertex in get_surface()\n", stderr);
257 +                                free(surf);
258 +                                return(NULL);
259 +                        }
260 +                }
261 +                                        /* compute area and normal */
262 +        surf->norm[0] = surf->norm[1] = surf->norm[2] = 0;
263 +        VSUB(e1, surf->vl[1], surf->vl[0]);
264 +        for (i = 2; i < nv; i++) {
265 +                VSUB(e2, surf->vl[i], surf->vl[0]);
266 +                fcross(vc, e1, e2);
267 +                surf->norm[0] += vc[0];
268 +                surf->norm[1] += vc[1];
269 +                surf->norm[2] += vc[2];
270 +                VCOPY(e1, e2);
271 +        }
272 +        surf->area = .5 * normalize(surf->norm);
273 +        if (surf->area == 0) {
274 +                fputs(progname, stderr);
275 +                fputs(": degenerate polygon in get_surface()\n", stderr);
276 +                free(surf);
277 +                return(NULL);
278 +        }
279 +        return(surf);
280 + }
281 +
282 + /* Convert surface to Radiance with modifier based on unique name */
283 + static int
284 + rad_surface(IDF_OBJECT *param, FILE *ofp)
285 + {
286 +        const char      *sname = idf_getfield(param,NAME_FLD)->val;
287 +        IDF_FIELD       *fptr = get_vlist(param, NULL);
288 +        const char      **fargs;
289 +        int             nvert, i, j;
290 +
291 +        if (fptr == NULL || (nvert = atoi(fptr->val)) < 3) {
292 +                fprintf(stderr, "%s: bad surface '%s'\n", progname, sname);
293 +                return(0);
294 +        }
295 +        fargs = (const char **)malloc(sizeof(const char *)*3*nvert);
296 +        if (fargs == NULL)
297 +                return(0);
298 +        fprintf(ofp, "\nvoid glow '%s'\n0\n0\n4 1 1 1 0\n", sname);
299 +        fprintf(ofp, "\n'%s' polygon 's_%s'\n0\n0\n%d\n", sname, sname, 3*nvert);
300 +        for (j = nvert; j--; ) {                /* get arguments in reverse */
301 +                for (i = 0; i < 3; i++) {
302 +                        fptr = fptr->next;
303 +                        if (fptr == NULL || !isflt(fptr->val)) {
304 +                                fprintf(stderr,
305 +                                        "%s: missing/bad vertex for %s '%s'\n",
306 +                                                progname, param->pname, sname);
307 +                                return(0);
308 +                        }
309 +                        fargs[3*j + i] = fptr->val;
310 +                }
311 +        }
312 +        for (j = 0; j < nvert; j++) {           /* output reversed verts */
313 +                for (i = 0; i < 3; i++) {
314 +                        fputc('\t', ofp);
315 +                        fputs(fargs[3*j + i], ofp);
316 +                }
317 +                fputc('\n', ofp);
318 +        }
319 +        free(fargs);
320 +        return(!ferror(ofp));
321 + }
322 +
323 + /* Convert subsurface to Radiance with modifier based on unique name */
324 + static double
325 + rad_subsurface(IDF_OBJECT *param, FILE *ofp)
326 + {
327 +        const char      *sname = idf_getfield(param,NAME_FLD)->val;
328 +        SURFACE         *surf = get_surface(get_vlist(param, NULL));
329 +        double          area;
330 +        int             i;
331 +
332 +        if (surf == NULL) {
333 +                fprintf(stderr, "%s: bad subsurface '%s'\n", progname, sname);
334 +                return(-1.);
335 +        }
336 +        fprintf(ofp, "\nvoid glow '%s'\n0\n0\n4 1 1 1 0\n", sname);
337 +        fprintf(ofp, "\n'%s' polygon 'ss_%s'\n0\n0\n%d\n",
338 +                                                sname, sname, 3*surf->nv);
339 +        for (i = 0; i < surf->nv; i++) {        /* offset vertices */
340 +                FVECT   vert;
341 +                VSUM(vert, surf->vl[i], surf->norm, 2.*SURF_EPS);
342 +                fprintf(ofp, "\t%.12f %.12f %.12f\n", vert[0], vert[1], vert[2]);
343 +        }
344 +        area = surf->area;
345 +        free(surf);
346 +        if (ferror(ofp))
347 +                return(-1.);
348 +        return(area);
349 + }
350 +
351 + /* Start rcontrib process */
352 + static int
353 + start_rcontrib(SUBPROC *pd, ZONE *zp)
354 + {
355 + #define BASE_AC         6
356 +        static char     *base_av[BASE_AC] = {
357 +                                "rcontrib", "-V+", "-fff", "-h", "-x", "1"
358 +                        };
359 +        char            cbuf[300];
360 +        char            **av;
361 +        FILE            *ofp;
362 +        IDF_OBJECT      *pptr;
363 +        IDF_FIELD       *fptr;
364 +        int             i, j, n;
365 +                                                /* start oconv command */
366 +        sprintf(cbuf, "oconv - > \"%s\"", temp_octree);
367 +        if ((ofp = popen(cbuf, "w")) == NULL) {
368 +                fputs(progname, stderr);
369 +                fputs(": cannot open oconv process\n", stderr);
370 +                return(0);
371 +        }
372 +                                                /* allocate argument list */
373 +        av = (char **)malloc(sizeof(char *)*(BASE_AC+4+2*(zp->ntotal)));
374 +        if (av == NULL)
375 +                return(0);
376 +        for (i = 0; i < BASE_AC; i++)
377 +                av[i] = base_av[i];
378 +        sprintf(cbuf, "%d", nsamps);
379 +        av[i++] = "-c";
380 +        av[i++] = cbuf;                         /* add modifiers & surfaces */
381 +        for (n = 0, pptr = zp->pfirst; n < zp->nsurf; n++, pptr = pptr->dnext) {
382 +                fptr = idf_getfield(pptr,NAME_FLD);
383 +                if (fptr == NULL || !fptr->val[0]) {
384 +                        fputs(progname, stderr);
385 +                        fputs(": missing name for surface object\n", stderr);
386 +                        return(0);
387 +                }
388 +                if (!rad_surface(pptr, ofp))    /* add surface to octree */
389 +                        return(0);
390 +                av[i++] = "-m";
391 +                av[i++] = fptr->val;
392 +        }
393 +                                                /* now subsurfaces */
394 +        if (zp->ntotal > zp->nsurf) {
395 +                if (zp->area_redu != NULL)
396 +                        memset(zp->area_redu, 0, sizeof(float)*zp->ntotal);
397 +                else if ((zp->area_redu = (float *)calloc(zp->ntotal,
398 +                                                sizeof(float))) == NULL)
399 +                        return(0);
400 +        }
401 +        for ( ; n < zp->ntotal; n++, pptr = pptr->dnext) {
402 +                double          ssarea;
403 +                const char      *bname;
404 +                IDF_OBJECT      *pptr1;
405 +                fptr = idf_getfield(pptr,NAME_FLD);
406 +                if (fptr == NULL || !fptr->val[0]) {
407 +                        fputs(progname, stderr);
408 +                        fputs(": missing name for subsurface object\n", stderr);
409 +                        return(0);
410 +                }
411 +                                                /* add subsurface to octree */
412 +                if ((ssarea = rad_subsurface(pptr, ofp)) < 0)
413 +                        return(0);
414 +                                                /* mark area for subtraction */
415 +                bname = idf_getfield(pptr,SS_BASE_FLD)->val;
416 +                for (j = 0, pptr1 = zp->pfirst;
417 +                                j < zp->nsurf; j++, pptr1 = pptr1->dnext)
418 +                        if (!strcmp(idf_getfield(pptr1,NAME_FLD)->val, bname)) {
419 +                                zp->area_redu[j] += ssarea;
420 +                                break;
421 +                        }
422 +                av[i++] = "-m";
423 +                av[i++] = fptr->val;
424 +        }
425 +        if (pclose(ofp) != 0) {                 /* finish oconv */
426 +                fputs(progname, stderr);
427 +                fputs(": error running oconv process\n", stderr);
428 +                return(0);
429 +        }
430 +        av[i++] = temp_octree;                  /* add final octree argument */
431 +        av[i] = NULL;
432 +        if (open_process(pd, av) <= 0) {        /* start process */
433 +                fputs(progname, stderr);
434 +                fputs(": cannot start rcontrib process\n", stderr);
435 +                return(0);
436 +        }
437 +        free(av);                               /* all done -- clean up */
438 +        return(1);
439 + #undef BASE_AC
440 + }
441 +
442 + /* Initialize polygon sampling */
443 + static Vert2_list *
444 + init_poly(POLYSAMP *ps, IDF_FIELD *f0)
445 + {
446 +        SURFACE         *surf;
447 +        Vert2_list      *vl2;
448 +        int             i;
449 +                                        /* get 3-D polygon vertices */
450 +        if ((surf = get_surface(f0)) == NULL)
451 +                return(NULL);
452 +        vl2 = polyAlloc(surf->nv);
453 +        if (vl2 == NULL)
454 +                return(NULL);
455 +                                        /* create X & Y axes */
456 +        VCOPY(ps->sdir[2], surf->norm);
457 +        VSUB(ps->sdir[0], surf->vl[1], surf->vl[0]);
458 +        if (normalize(ps->sdir[0]) == 0)
459 +                return(NULL);
460 +        fcross(ps->sdir[1], ps->sdir[2], ps->sdir[0]);
461 +                                        /* compute plane offset */
462 +        ps->poff = DOT(surf->vl[0], ps->sdir[2]);
463 +                                        /* assign 2-D vertices */
464 +        for (i = 0; i < surf->nv; i++) {
465 +                vl2->v[i].mX = DOT(surf->vl[i], ps->sdir[0]);
466 +                vl2->v[i].mY = DOT(surf->vl[i], ps->sdir[1]);
467 +        }
468 +        ps->area_left = surf->area;
469 +        free(surf);                     /* it's ready! */
470 +        vl2->p = (void *)ps;
471 +        return(vl2);
472 + }
473 +
474 + /* Generate samples on 2-D triangle */
475 + static int
476 + sample_triangle(const Vert2_list *vl2, int a, int b, int c)
477 + {
478 +        POLYSAMP        *ps = (POLYSAMP *)vl2->p;
479 +        float           *samp;
480 +        FVECT           orig;
481 +        FVECT           ab, ac;
482 +        double          area;
483 +        int             i, j, ns;
484 +                                        /* compute sampling axes */
485 +        for (i = 3; i--; ) {
486 +                orig[i] = vl2->v[a].mX*ps->sdir[0][i] +
487 +                                vl2->v[a].mY*ps->sdir[1][i] +
488 +                                (ps->poff+SURF_EPS)*ps->sdir[2][i];
489 +                ab[i] = (vl2->v[b].mX - vl2->v[a].mX)*ps->sdir[0][i] +
490 +                                (vl2->v[b].mY - vl2->v[a].mY)*ps->sdir[1][i];
491 +                ac[i] = (vl2->v[c].mX - vl2->v[a].mX)*ps->sdir[0][i] +
492 +                                (vl2->v[c].mY - vl2->v[a].mY)*ps->sdir[1][i];
493 +        }
494 +                                        /* compute number of samples to take */
495 +        area = .5*(vl2->v[a].mX*vl2->v[b].mY - vl2->v[b].mX*vl2->v[a].mY +
496 +                        vl2->v[b].mX*vl2->v[c].mY - vl2->v[c].mX*vl2->v[b].mY +
497 +                        vl2->v[c].mX*vl2->v[a].mY - vl2->v[a].mX*vl2->v[c].mY);
498 +        if (area < .0) {
499 +                fputs(progname, stderr);
500 +                fputs(": negative triangle area in sample_triangle()\n", stderr);
501 +                return(0);
502 +        }
503 +        if (area >= ps->area_left) {
504 +                ns = ps->samp_left;
505 +                ps->area_left = 0;
506 +        } else {
507 +                ns = (ps->samp_left*area/ps->area_left + .5);
508 +                ps->samp_left -= ns;
509 +                ps->area_left -= area;
510 +        }
511 +        if (ns <= 0)                    /* XXX should be error? */
512 +                return(1);
513 +                                        /* buffer sample rays */
514 +        samp = (float *)malloc(sizeof(float)*6*ns);
515 +        if (samp == NULL)
516 +                return(0);
517 +        for (i = ns; i--; ) {           /* stratified Monte Carlo sampling */
518 +                double  sv[4];
519 +                FVECT   dv;
520 +                multisamp(sv, 4, (i+frandom())/(double)ns);
521 +                sv[0] *= sv[1] = sqrt(sv[1]);
522 +                sv[1] = 1. - sv[1];
523 +                for (j = 3; j--; )
524 +                        samp[i*6 + j] = orig[j] + sv[0]*ab[j] + sv[1]*ac[j];
525 +                sv[2] = sqrt(sv[2]);
526 +                sv[3] *= 2.*PI;
527 +                dv[0] = tcos(sv[3]) * sv[2];
528 +                dv[1] = tsin(sv[3]) * sv[2];
529 +                dv[2] = sqrt(1. - sv[2]*sv[2]);
530 +                for (j = 3; j--; )
531 +                        samp[i*6 + 3 + j] = dv[0]*ps->sdir[0][j] +
532 +                                                dv[1]*ps->sdir[1][j] +
533 +                                                dv[2]*ps->sdir[2][j] ;
534 +        }
535 +                                        /* send to our process */
536 +        writebuf(ps->wd, (char *)samp, sizeof(float)*6*ns);
537 +        free(samp);                     /* that's it! */
538 +        return(1);
539 + }
540 +
541 + /* Sample the given surface */
542 + static double
543 + sample_surface(IDF_OBJECT *param, int wd)
544 + {
545 +        POLYSAMP        psamp;
546 +        double          area;
547 +        int             nv;
548 +        Vert2_list      *vlist2;
549 +                                        /* set up our polygon sampler */
550 +        if ((vlist2 = init_poly(&psamp, get_vlist(param, NULL))) == NULL) {
551 +                fprintf(stderr, "%s: bad polygon %s '%s'\n",
552 +                                progname, param->pname,
553 +                                idf_getfield(param,NAME_FLD)->val);
554 +                return(-1.);
555 +        }
556 +        psamp.samp_left = nsamps;       /* assign samples & destination */
557 +        psamp.wd = wd;
558 +                                        /* hack for subsurface sampling */
559 +        psamp.poff += 2.*SURF_EPS * !strcmp(param->pname, SUBSURF_PNAME);
560 +
561 +        area = psamp.area_left;         /* remember starting surface area */
562 +                                        /* sample each subtriangle */
563 +        if (!polyTriangulate(vlist2, &sample_triangle))
564 +                return(-1.);
565 +        polyFree(vlist2);               /* clean up and return */
566 +        return(area);
567 + }
568 +
569 + /* Compute User View Factors using open rcontrib process */
570 + static int
571 + compute_uvfs(SUBPROC *pd, ZONE *zp)
572 + {
573 +        IDF_OBJECT      *pptr, *pout, *pptr1;
574 +        float           *uvfa;
575 +        char            uvfbuf[24];
576 +        int             n, m;
577 +                                                /* create output object */
578 +        pout = idf_newobject(our_idf, UVF_PNAME,
579 +                        "    ! computed by Radiance\n        ", our_idf->plast);
580 +        if (pout == NULL) {
581 +                fputs(progname, stderr);
582 +                fputs(": cannot create new IDF object\n", stderr);
583 +                return(0);
584 +        }
585 +        if (!idf_addfield(pout, zp->zname,
586 +                        "    !- Zone Name\n        ")) {
587 +                fputs(progname, stderr);
588 +                fputs(": cannot add zone name field\n", stderr);
589 +                return(0);
590 +        }
591 +                                                /* allocate read buffer */
592 +        uvfa = (float *)malloc(sizeof(float)*3*zp->ntotal);
593 +        if (uvfa == NULL)
594 +                return(0);
595 +                                                /* UVFs from each surface */
596 +        for (n = 0, pptr = zp->pfirst; n < zp->ntotal; n++, pptr = pptr->dnext) {
597 +                double  vfsum = 0;
598 +                double  adj_factor;
599 +                                                /* send samples to rcontrib */
600 +                if ((adj_factor = sample_surface(pptr, pd->w)) < 0)
601 +                        return(0);
602 +                if (zp->area_redu == NULL)
603 +                        adj_factor = 1.;
604 +                else                            /* comp. for subsurface area */
605 +                        adj_factor /= adj_factor - zp->area_redu[n];
606 +                                                /* read results */
607 +                if (readbuf(pd->r, (char *)uvfa, sizeof(float)*3*zp->ntotal) !=
608 +                                sizeof(float)*3*zp->ntotal) {
609 +                        fputs(progname, stderr);
610 +                        fputs(": read error from rcontrib process\n", stderr);
611 +                        return(0);
612 +                }
613 +                                                /* append UVF fields */
614 +                for (m = 0, pptr1 = zp->pfirst;
615 +                                m < zp->ntotal; m++, pptr1 = pptr1->dnext) {
616 +                        const double    uvf = uvfa[3*m + 1] * adj_factor;
617 +                        vfsum += uvf;
618 +                        if (pptr1 == pptr && uvf > .001)
619 +                                fprintf(stderr,
620 +                "%s: warning - non-zero self-VF (%.1f%%) for surface '%s'\n",
621 +                                                progname, 100.*uvf,
622 +                                                idf_getfield(pptr,NAME_FLD)->val);
623 +                        sprintf(uvfbuf, "%.4f", uvf);
624 +                        if (!idf_addfield(pout,
625 +                                        idf_getfield(pptr,NAME_FLD)->val, NULL) ||
626 +                                !idf_addfield(pout,
627 +                                        idf_getfield(pptr1,NAME_FLD)->val, NULL) ||
628 +                                !idf_addfield(pout, uvfbuf,
629 +                                                (n+m < 2*zp->ntotal-2) ?
630 +                                                "\n        " : "\n\n")) {
631 +                                fputs(progname, stderr);
632 +                                fputs(": error adding UVF fields\n", stderr);
633 +                                return(0);
634 +                        }
635 +                }
636 +                if (vfsum < 0.95)
637 +                        fprintf(stderr,
638 +                "%s: warning - missing %.1f%% of energy from surface '%s'\n",
639 +                                        progname, 100.*(1.-vfsum),
640 +                                        idf_getfield(pptr,NAME_FLD)->val);
641 +        }
642 +        free(uvfa);                             /* clean up and return */
643 +        return(1);
644 + }
645 +
646   /* Compute zone User View Factors */
647   static int
648   compute_zones(void)
649   {
650 <        fputs("compute_zones() unimplemented\n", stderr);
651 <        return(0);
650 >        ZONE    *zptr;
651 >                                                /* temporary octree name */
652 >        mktemp(strcpy(temp_octree, TEMPLATE));
653 >                                                /* compute each zone */
654 >        for (zptr = zone_list; zptr != NULL; zptr = zptr->next) {
655 >                SUBPROC rcproc;
656 >                                                /* start rcontrib process */
657 >                if (!start_rcontrib(&rcproc, zptr))
658 >                        return(0);
659 >                                                /* compute+add view factors */
660 >                if (!compute_uvfs(&rcproc, zptr))
661 >                        return(0);
662 >                if (close_process(&rcproc) != 0) {
663 >                        fputs(progname, stderr);
664 >                        fputs(": bad return status from rcontrib\n", stderr);
665 >                        return(0);
666 >                }
667 >        }
668 >        unlink(temp_octree);                    /* remove octree file */
669 >        return(1);
670   }
671  
672   /* Load IDF and compute User View Factors */
673   int
674   main(int argc, char *argv[])
675   {
676 +        int             incl_comments = 1;
677          char            *origIDF, *revIDF;
678 <        IDF_PARAMETER   *pptr;
678 >        IDF_OBJECT      *pptr;
679          int             i;
680  
681 <        if ((argc < 2) | (argc > 3)) {
682 <                fputs("Usage: ", stderr);
683 <                fputs(argv[0], stderr);
684 <                fputs(" Model.idf [Revised.idf]\n", stderr);
685 <                return(1);
686 <        }
687 <        origIDF = argv[1];
688 <        revIDF = (argc == 2) ? argv[1] : argv[2];
681 >        progname = *argv++; argc--;             /* get options if any */
682 >        while (argc > 1 && argv[0][0] == '-')
683 >                switch (argv[0][1]) {
684 >                case 'c':                       /* elide comments */
685 >                        incl_comments = -1;             /* header only */
686 >                        argv++; argc--;
687 >                        continue;
688 >                case 's':                       /* samples */
689 >                        nsamps = 1000*atoi(*++argv);
690 >                        argv++; argc -= 2;
691 >                        continue;
692 >                default:
693 >                        fputs(progname, stderr);
694 >                        fputs(": unknown option '", stderr);
695 >                        fputs(argv[0], stderr);
696 >                        fputs("'\n", stderr);
697 >                        goto userr;
698 >                }
699 >        if ((argc < 1) | (argc > 2))
700 >                goto userr;
701 >        origIDF = argv[0];
702 >        revIDF = (argc == 1) ? argv[0] : argv[1];
703                                                  /* load Input Data File */
704          our_idf = idf_load(origIDF);
705          if (our_idf == NULL) {
706 <                fputs(argv[0], stderr);
706 >                fputs(progname, stderr);
707                  fputs(": cannot load IDF '", stderr);
708                  fputs(origIDF, stderr);
709                  fputs("'\n", stderr);
710                  return(1);
711          }
712                                                  /* remove existing UVFs */
713 <        if ((pptr = idf_getparam(our_idf, UVF_PNAME)) != NULL) {
714 <                IDF_PARAMETER   *pnext;
715 <                fputs(argv[0], stderr);
713 >        if ((pptr = idf_getobject(our_idf, UVF_PNAME)) != NULL) {
714 >                IDF_OBJECT      *pnext;
715 >                fputs(progname, stderr);
716                  fputs(": removing previous User View Factors\n", stderr);
717                  do {
718                          pnext = pptr->pnext;
719 <                        idf_delparam(our_idf, pptr);
719 >                        idf_delobject(our_idf, pptr);
720                  } while (pnext != NULL);
721          }
722                                                  /* add to header */
# Line 126 | Line 726 | main(int argc, char *argv[])
726                  idf_add2hdr(our_idf, ADD_HEADER);
727                                                  /* gather zone surfaces */
728          for (i = 0; surf_type[i].pname != NULL; i++)
729 <                for (pptr = idf_getparam(our_idf, surf_type[i].pname);
729 >                for (pptr = idf_getobject(our_idf, surf_type[i].pname);
730                                  pptr != NULL; pptr = pptr->pnext) {
731                          IDF_FIELD       *fptr = idf_getfield(pptr,
732                                                          surf_type[i].zone_fld);
733                          if (fptr == NULL) {
734 <                                fputs("Warning: missing zone field!\n", stderr);
734 >                                fputs(progname, stderr);
735 >                                fputs(": warning - missing zone field\n", stderr);
736                                  continue;
737                          }
738 <                        if (add2zone(pptr, fptr->arg) == NULL)
738 >                        if (add2zone(pptr, fptr->val) == NULL)
739                                  return(1);
740                  }
741 +                                                /* add subsurfaces */
742 +        for (pptr = idf_getobject(our_idf, SUBSURF_PNAME);
743 +                        pptr != NULL; pptr = pptr->pnext)
744 +                if (add_subsurf(pptr) == NULL)
745 +                        return(1);
746                                                  /* run rcontrib on each zone */
747          if (!compute_zones())
748                  return(1);
749                                                  /* write out modified IDF */
750 <        if (!idf_write(our_idf, revIDF, 1)) {
751 <                fputs(argv[0], stderr);
750 >        if (!idf_write(our_idf, revIDF, incl_comments)) {
751 >                fputs(progname, stderr);
752                  fputs(": error writing IDF '", stderr);
753                  fputs(revIDF, stderr);
754                  fputs("'\n", stderr);
755                  return(1);
756          }
757          return(0);                              /* finito! */
758 + userr:
759 +        fputs("Usage: ", stderr);
760 +        fputs(progname, stderr);
761 +        fputs(" [-c][-s Ksamps] Model.idf [Revised.idf]\n", stderr);
762 +        return(1);
763   }

Diff Legend

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