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.2 by greg, Sun Feb 9 22:19:30 2014 UTC vs.
Revision 2.5 by greg, Tue Feb 11 21:17:29 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        10000                   /* number of samples to use */
20   #endif
21  
22   char            *progname;                      /* global argv[0] */
23  
24 < char            temp_octree[128];                       /* temporary octree */
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 38 | 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  
# Line 87 | Line 91 | add2zone(IDF_PARAMETER *param, const char *zname)
91          return(zptr);
92   }
93  
94 < /* Determine if a parameter is a surface in the indicated zone */
95 < static int
96 < in_zone(IDF_PARAMETER *param, const char *zname)
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(0);
104 <                                                /* get zone field */
105 <        fptr = idf_getfield(param, surf_type[i].zone_fld);
106 <        if (fptr == NULL)
107 <                return(0);
108 <                                                /* check for match */
109 <        if (strcmp(fptr->arg, zname))
106 <                return(0);
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(surf_type[i].vert_fld);
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(const char *zname, IDF_PARAMETER *param, FILE *ofp)
116 > rad_surface(IDF_PARAMETER *param, FILE *ofp)
117   {
118 <        const char      *sname = idf_getfield(param, NAME_FLD)->arg;
119 <        IDF_FIELD       *fptr = idf_getfield(param, in_zone(param, zname));
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->arg)) < 3) {
122 >        if (fptr == NULL || (nvert = atoi(fptr->val)) < 3) {
123                  fprintf(stderr, "%s: bad surface '%s'\n", progname, sname);
124                  return(0);
125          }
# Line 125 | Line 128 | rad_surface(const char *zname, IDF_PARAMETER *param, F
128          while (nvert--) {
129                  for (i = 3; i--; ) {
130                          fptr = fptr->next;
131 <                        if (fptr == NULL) {
131 >                        if (fptr == NULL || !isflt(fptr->val)) {
132                                  fprintf(stderr,
133 <                                "%s: missing vertex fields in surface '%s'\n",
134 <                                                progname, sname);
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->arg, ofp);
138 >                        fputs(fptr->val, ofp);
139                  }
140                  fputc('\n', ofp);
141          }
142          return(!ferror(ofp));
143   }
144  
145 < /* Strat rcontrib process */
145 > /* Start rcontrib process */
146   static int
147   start_rcontrib(SUBPROC *pd, ZONE *zp)
148   {
149 < #define BASE_AC         3
149 > #define BASE_AC         5
150          static char     *base_av[BASE_AC] = {
151 <                                "rcontrib", "-ff", "-h"
151 >                                "rcontrib", "-ff", "-h", "-x", "1"
152                          };
153          char            cbuf[300];
154          char            **av;
# Line 169 | Line 172 | start_rcontrib(SUBPROC *pd, ZONE *zp)
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) {
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 <                av[i++] = "-m";
179 <                av[i++] = fptr->arg;
180 <                if (!rad_surface(zp->zname, pptr, ofp))
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);
# Line 197 | Line 200 | start_rcontrib(SUBPROC *pd, ZONE *zp)
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 *sp)
363 > compute_uvfs(SUBPROC *pd, ZONE *zp)
364   {
365 <        
365 >        IDF_PARAMETER   *pptr, *pout, *pptr1;
366 >        float           *uvfa;
367 >        char            uvfbuf[24];
368 >        int             n, m;
369 >                                                /* create output parameter */
370 >        pout = idf_newparam(our_idf, UVF_PNAME,
371 >                        "    ! computed by Radiance\n        ", zp->plast);
372 >        if (pout == NULL) {
373 >                fputs(progname, stderr);
374 >                fputs(": cannot create new IDF parameter\n", stderr);
375 >                return(0);
376 >        }
377 >        if (!idf_addfield(pout, zp->zname,
378 >                        "    !- Zone Name\n        ")) {
379 >                fputs(progname, stderr);
380 >                fputs(": cannot add zone name field\n", stderr);
381 >                return(0);
382 >        }
383 >                                                /* allocate read buffer */
384 >        uvfa = (float *)malloc(sizeof(float)*3*zp->nsurf);
385 >        if (uvfa == NULL)
386 >                return(0);
387 >                                                /* UVFs from each surface */
388 >        for (n = zp->nsurf, pptr = zp->pfirst; n--; pptr = pptr->dnext) {
389 >                double  vfsum = 0;
390 >                                                /* send samples to rcontrib */
391 >                if (!sample_surface(pptr, pd->w))
392 >                        return(0);
393 >                                                /* read results */
394 >                if (readbuf(pd->r, (char *)uvfa, sizeof(float)*3*zp->nsurf) !=
395 >                                sizeof(float)*3*zp->nsurf) {
396 >                        fputs(progname, stderr);
397 >                        fputs(": read error from rcontrib process\n", stderr);
398 >                        return(0);
399 >                }
400 >                                                /* append UVF fields */
401 >                for (m = 0, pptr1 = zp->pfirst;
402 >                                m < zp->nsurf; m++, pptr1 = pptr1->dnext) {
403 >                        vfsum += uvfa[3*m + 1];
404 >                        if (pptr1 == pptr) {
405 >                                if (uvfa[3*m + 1] > .001)
406 >                                        fprintf(stderr,
407 >                "%s: warning - non-zero self-VF (%.1f%%) for surface '%s'\n",
408 >                                                progname, 100.*uvfa[3*m + 1],
409 >                                                idf_getfield(pptr,NAME_FLD)->val);
410 >                                continue;       /* don't record self-factor */
411 >                        }
412 >                        sprintf(uvfbuf, "%.6f", uvfa[3*m + 1]);
413 >                        if (!idf_addfield(pout,
414 >                                        idf_getfield(pptr,NAME_FLD)->val, NULL) ||
415 >                                !idf_addfield(pout,
416 >                                        idf_getfield(pptr1,NAME_FLD)->val, NULL) ||
417 >                                !idf_addfield(pout, uvfbuf,
418 >                                                (n || m < zp->nsurf-2) ?
419 >                                                        "\n        " : "\n\n")) {
420 >                                fputs(progname, stderr);
421 >                                fputs(": error adding UVF fields\n", stderr);
422 >                                return(0);
423 >                        }
424 >                }
425 >                if (vfsum < 0.95)
426 >                        fprintf(stderr,
427 >                "%s: warning - missing %.1f%% of energy from surface '%s'\n",
428 >                                        progname, 100.*(1.-vfsum),
429 >                                        idf_getfield(pptr,NAME_FLD)->val);
430 >        }
431 >        free(uvfa);                             /* clean up and return */
432 >        return(1);
433   }
434  
435   /* Compute zone User View Factors */
# Line 210 | Line 438 | compute_zones(void)
438   {
439          ZONE    *zptr;
440                                                  /* temporary octree name */
441 <        if (temp_filename(temp_octree, sizeof(temp_octree), TEMPLATE) == NULL) {
214 <                fputs(progname, stderr);
215 <                fputs(": cannot create temporary octree\n", stderr);
216 <                return(0);
217 <        }
441 >        mktemp(strcpy(temp_octree, TEMPLATE));
442                                                  /* compute each zone */
443          for (zptr = zone_list; zptr != NULL; zptr = zptr->next) {
444                  SUBPROC rcproc;
# Line 291 | Line 515 | main(int argc, char *argv[])
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 */

Diff Legend

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