src/util/eplus_adduvf.c

#ifndef lint
static const char RCSid[] = "$Id: eplus_adduvf.c,v 2.5 2014/02/11 21:17:29 greg Exp $";
#endif
/*
 * Add User View Factors to EnergyPlus Input Data File
 *
 *      G.Ward for LBNL
 */

#include <stdlib.h>
#include "rtio.h"
#include "rtmath.h"
#include "random.h"
#include "eplus_idf.h"
#include "triangulate.h"
#include "rtprocess.h"

#ifndef NSAMPLES
#define NSAMPLES        80000                   /* number of samples to use */
#endif

char            *progname;                      /* global argv[0] */

char            temp_octree[128];               /* temporary octree */

const char      UVF_PNAME[] =
                        "ZoneProperty:UserViewFactor:bySurfaceName";

const char      ADD_HEADER[] =
                        "\n!+++ User View Factors computed by Radiance +++!\n\n";

#define NAME_FLD        1                       /* name field always first? */

typedef struct {
        const char      *pname;                 /* parameter type name */
        short           zone_fld;               /* zone field index */
        short           vert_fld;               /* vertex field index */
} SURF_PTYPE;           /* surface type we're interested in */

const SURF_PTYPE        surf_type[] = {
                {"BuildingSurface:Detailed", 4, 10},
                {"Floor:Detailed", 3, 9},
                {"RoofCeiling:Detailed", 3, 9},
                {"Wall:Detailed", 3, 9},
                {NULL}
        };

typedef struct s_zone {
        const char      *zname;                 /* zone name */
        struct s_zone   *next;                  /* next zone in list */
        int             nsurf;                  /* surface count */
        IDF_PARAMETER   *pfirst;                /* first matching parameter */
        IDF_PARAMETER   *plast;                 /* last matching parameter */
} ZONE;                 /* a list of collected zone surfaces */

ZONE            *zone_list = NULL;      /* our list of zones */

IDF_LOADED      *our_idf = NULL;        /* loaded/modified IDF */

/* Create a new zone and push to top of our list */
static ZONE *
new_zone(const char *zname, IDF_PARAMETER *param)
{
        ZONE    *znew = (ZONE *)malloc(sizeof(ZONE));

        if (znew == NULL)
                return(NULL);
        znew->zname = zname;                    /* assumes static */
        znew->next = zone_list;
        znew->pfirst = znew->plast = param;
        znew->nsurf = 1;
        return(zone_list = znew);
}

/* Add the detailed surface (polygon) to the named zone */
static ZONE *
add2zone(IDF_PARAMETER *param, const char *zname)
{
        ZONE    *zptr;

        for (zptr = zone_list; zptr != NULL; zptr = zptr->next)
                if (!strcmp(zptr->zname, zname))
                        break;
        if (zptr == NULL)
                return(new_zone(zname, param));
                                                /* keep surfaces together */
        if (!idf_movparam(our_idf, param, zptr->plast))
                return(NULL);
        zptr->plast = param;
        zptr->nsurf++;
        return(zptr);
}

/* Return field for vertices in the given parameter */
static IDF_FIELD *
get_vlist(IDF_PARAMETER *param, const char *zname)
{
        int             i = 0;
        IDF_FIELD       *fptr;
                                                /* check for surface type */
        while (strcmp(surf_type[i].pname, param->pname))
                if (surf_type[++i].pname == NULL)
                        return(NULL);

        if (zname != NULL) {                    /* matches specified zone? */
                fptr = idf_getfield(param, surf_type[i].zone_fld);
                if (fptr == NULL || strcmp(fptr->val, zname))
                        return(NULL);
        }
                                                /* return field for #verts */
        return(idf_getfield(param, surf_type[i].vert_fld));
}

/* Convert surface to Radiance with modifier based on unique name */
static int
rad_surface(IDF_PARAMETER *param, FILE *ofp)
{
        const char      *sname = idf_getfield(param,NAME_FLD)->val;
        IDF_FIELD       *fptr = get_vlist(param, NULL);
        int             nvert, i;

        if (fptr == NULL || (nvert = atoi(fptr->val)) < 3) {
                fprintf(stderr, "%s: bad surface '%s'\n", progname, sname);
                return(0);
        }
        fprintf(ofp, "\nvoid glow '%s'\n0\n0\n4 1 1 1 0\n", sname);
        fprintf(ofp, "\n'%s' polygon 's_%s'\n0\n0\n%d\n", sname, sname, 3*nvert);
        while (nvert--) {
                for (i = 3; i--; ) {
                        fptr = fptr->next;
                        if (fptr == NULL || !isflt(fptr->val)) {
                                fprintf(stderr,
                                        "%s: missing/bad vertex for %s '%s'\n",
                                                progname, param->pname, sname);
                                return(0);
                        }
                        fputc('\t', ofp);
                        fputs(fptr->val, ofp);
                }
                fputc('\n', ofp);
        }
        return(!ferror(ofp));
}

/* Start rcontrib process */
static int
start_rcontrib(SUBPROC *pd, ZONE *zp)
{
#define BASE_AC         5
        static char     *base_av[BASE_AC] = {
                                "rcontrib", "-ff", "-h", "-x", "1"
                        };
        char            cbuf[300];
        char            **av;
        FILE            *ofp;
        IDF_PARAMETER   *pptr;
        int             i, n;
                                                /* start oconv command */
        sprintf(cbuf, "oconv - > '%s'", temp_octree);
        if ((ofp = popen(cbuf, "w")) == NULL) {
                fputs(progname, stderr);
                fputs(": cannot open oconv process\n", stderr);
                return(0);
        }
                                                /* allocate argument list */
        av = (char **)malloc(sizeof(char *)*(BASE_AC+4+2*zp->nsurf));
        if (av == NULL)
                return(0);
        for (i = 0; i < BASE_AC; i++)
                av[i] = base_av[i];
        sprintf(cbuf, "%d", NSAMPLES);
        av[i++] = "-c";
        av[i++] = cbuf;                         /* add modifier arguments */
        for (n = zp->nsurf, pptr = zp->pfirst; n--; pptr = pptr->dnext) {
                IDF_FIELD       *fptr = idf_getfield(pptr,NAME_FLD);
                if (fptr == NULL || !fptr->val[0]) {
                        fputs(progname, stderr);
                        fputs(": missing name for surface parameter\n", stderr);
                        return(0);
                }
                if (!rad_surface(pptr, ofp))    /* add surface to octree */
                        return(0);
                av[i++] = "-m";
                av[i++] = fptr->val;
        }
        if (pclose(ofp) != 0) {                 /* finish oconv */
                fputs(progname, stderr);
                fputs(": error running oconv process\n", stderr);
                return(0);
        }
        av[i++] = temp_octree;                  /* add final octree argument */
        av[i] = NULL;
        if (!open_process(pd, av)) {            /* start process */
                fputs(progname, stderr);
                fputs(": cannot start rcontrib process\n", stderr);
                return(0);
        }
        free(av);                               /* all done -- clean up */
        return(1);
#undef BASE_AC
}

typedef struct {
        FVECT           sdir[3];        /* XYZ unit sampling vectors */
        double          poff;           /* Z-offset for plane of polygon */
        double          area_left;      /* area left to sample */
        int             samp_left;      /* remaining samples */
        int             wd;             /* output file descriptor */
} POLYSAMP;             /* structure for polygon sampling */

/* Initialize polygon sampling */
static Vert2_list *
init_poly(POLYSAMP *ps, IDF_FIELD *f0, int nv)
{
        IDF_FIELD       *fptr = f0;
        int             i, j;
        FVECT           *vl3, e1, e2, vc;
        Vert2_list      *vl2 = polyAlloc(nv);

        if (vl2 == NULL)
                return(NULL);
        vl2->p = ps;
                                        /* get 3-D vertices */
        vl3 = (FVECT *)malloc(sizeof(FVECT)*nv);
        if (vl3 == NULL)
                return(NULL);
        for (i = nv; i--; )             /* reverse vertex ordering */
                for (j = 0; j < 3; j++) {
                        if (fptr == NULL) {
                                fputs(progname, stderr);
                                fputs(": missing vertex in init_poly()\n", stderr);
                                return(NULL);
                        }
                        vl3[i][j] = atof(fptr->val);
                        fptr = fptr->next;
                }
                                        /* compute area and normal */
        ps->sdir[2][0] = ps->sdir[2][1] = ps->sdir[2][2] = 0;
        VSUB(e1, vl3[1], vl3[0]);
        for (i = 2; i < nv; i++) {
                VSUB(e2, vl3[i], vl3[0]); 
                fcross(vc, e1, e2);
                ps->sdir[2][0] += vc[0];
                ps->sdir[2][1] += vc[1];
                ps->sdir[2][2] += vc[2];
                VCOPY(e1, e2);
        }
        ps->area_left = .5 * normalize(ps->sdir[2]);
        if (ps->area_left == .0) {
                fputs(progname, stderr);
                fputs(": degenerate polygon in init_poly()\n", stderr);
                return(0);
        }
                                        /* create X & Y axes */
        VCOPY(ps->sdir[0], e1);
        normalize(ps->sdir[0]);
        fcross(ps->sdir[1], ps->sdir[2], ps->sdir[0]);
                                        /* compute plane offset */
        ps->poff = DOT(vl3[0], ps->sdir[2]);
                                        /* assign 2-D vertices */
        for (i = 0; i < nv; i++) {
                vl2->v[i].mX = DOT(vl3[i], ps->sdir[0]);
                vl2->v[i].mY = DOT(vl3[i], ps->sdir[1]);
        }
        free(vl3);                      /* it's ready! */
        return(vl2);
}

/* Generate samples on 2-D triangle */
static int
sample_triangle(const Vert2_list *vl2, int a, int b, int c)
{
        POLYSAMP        *ps = (POLYSAMP *)vl2->p;
        float           *samp;
        FVECT           orig;
        FVECT           ab, ac;
        double          area;
        int             i, j, ns;
                                        /* compute sampling axes */
        for (i = 3; i--; ) {
                orig[i] = vl2->v[a].mX*ps->sdir[0][i] +
                                vl2->v[a].mY*ps->sdir[1][i] +
                                (ps->poff+.001)*ps->sdir[2][i];
                ab[i] = (vl2->v[b].mX - vl2->v[a].mX)*ps->sdir[0][i] +
                                (vl2->v[b].mY - vl2->v[a].mY)*ps->sdir[1][i];
                ac[i] = (vl2->v[c].mX - vl2->v[a].mX)*ps->sdir[0][i] +
                                (vl2->v[c].mY - vl2->v[a].mY)*ps->sdir[1][i];
        }
                                        /* compute number of samples to take */
        area = .5*(vl2->v[a].mX*vl2->v[b].mY - vl2->v[b].mX*vl2->v[a].mY +
                        vl2->v[b].mX*vl2->v[c].mY - vl2->v[c].mX*vl2->v[b].mY +
                        vl2->v[c].mX*vl2->v[a].mY - vl2->v[a].mX*vl2->v[c].mY);
        if (area < .0) {
                fputs(progname, stderr);
                fputs(": negative triangle area in sample_triangle()\n", stderr);
                return(0);
        }
        if (area >= ps->area_left) {
                ns = ps->samp_left;
                ps->area_left = 0;
        } else {
                ns = (ps->samp_left*area/ps->area_left + .5);
                ps->samp_left -= ns;
                ps->area_left -= area;
        }
        if (ns <= 0)                    /* XXX should be error? */
                return(1);
                                        /* buffer sample rays */
        samp = (float *)malloc(sizeof(float)*6*ns);
        if (samp == NULL)
                return(0);
        for (i = ns; i--; ) {           /* stratified Monte Carlo sampling */
                double  sv[4];
                FVECT   dv;
                multisamp(sv, 4, (i+frandom())/(double)ns);
                sv[0] *= sv[1] = sqrt(sv[1]);
                sv[1] = 1. - sv[1];
                for (j = 3; j--; )
                        samp[i*6 + j] = orig[j] + sv[0]*ab[j] + sv[1]*ac[j];
                sv[2] = sqrt(sv[2]);
                sv[3] *= 2.*PI;
                dv[0] = tcos(sv[3]) * sv[2];
                dv[1] = tsin(sv[3]) * sv[2];
                dv[2] = sqrt(1. - sv[2]*sv[2]);
                for (j = 3; j--; )
                        samp[i*6 + 3 + j] = dv[0]*ps->sdir[0][j] +
                                                dv[1]*ps->sdir[1][j] +
                                                dv[2]*ps->sdir[2][j] ;
        }
                                        /* send to our process */
        writebuf(ps->wd, (char *)samp, sizeof(float)*6*ns);
        free(samp);                     /* that's it! */
        return(1);
}

/* Sample the given surface */
static int
sample_surface(IDF_PARAMETER *param, int wd)
{
        IDF_FIELD       *fptr = get_vlist(param, NULL);
        POLYSAMP        psamp;
        int             nv;
        Vert2_list      *vlist2;
                                        /* set up our polygon sampler */
        if (fptr == NULL || (nv = atoi(fptr->val)) < 3 ||
                        (vlist2 = init_poly(&psamp, fptr->next, nv)) == NULL) {
                fprintf(stderr, "%s: bad polygon %s '%s'\n",
                                progname, param->pname,
                                idf_getfield(param,NAME_FLD)->val);
                return(0);
        }
        psamp.samp_left = NSAMPLES;     /* assign samples & destination */
        psamp.wd = wd;
                                        /* sample each subtriangle */
        if (!polyTriangulate(vlist2, &sample_triangle))
                return(0);
        polyFree(vlist2);               /* clean up and return */
        return(1);
}

/* Compute User View Factors using open rcontrib process */
static int
compute_uvfs(SUBPROC *pd, ZONE *zp)
{
        IDF_PARAMETER   *pptr, *pout, *pptr1;
        float           *uvfa;
        char            uvfbuf[24];
        int             n, m;
                                                /* create output parameter */
        pout = idf_newparam(our_idf, UVF_PNAME,
                        "    ! computed by Radiance\n        ", zp->plast);
        if (pout == NULL) {
                fputs(progname, stderr);
                fputs(": cannot create new IDF parameter\n", stderr);
                return(0);
        }
        if (!idf_addfield(pout, zp->zname,
                        "    !- Zone Name\n        ")) {
                fputs(progname, stderr);
                fputs(": cannot add zone name field\n", stderr);
                return(0);
        }
                                                /* allocate read buffer */
        uvfa = (float *)malloc(sizeof(float)*3*zp->nsurf);
        if (uvfa == NULL)
                return(0);
                                                /* UVFs from each surface */
        for (n = zp->nsurf, pptr = zp->pfirst; n--; pptr = pptr->dnext) {
                double  vfsum = 0;
                                                /* send samples to rcontrib */
                if (!sample_surface(pptr, pd->w))
                        return(0);
                                                /* read results */
                if (readbuf(pd->r, (char *)uvfa, sizeof(float)*3*zp->nsurf) !=
                                sizeof(float)*3*zp->nsurf) {
                        fputs(progname, stderr);
                        fputs(": read error from rcontrib process\n", stderr);
                        return(0);
                }
                                                /* append UVF fields */
                for (m = 0, pptr1 = zp->pfirst;
                                m < zp->nsurf; m++, pptr1 = pptr1->dnext) {
                        vfsum += uvfa[3*m + 1];
                        if (pptr1 == pptr) {
                                if (uvfa[3*m + 1] > .001)
                                        fprintf(stderr,
                "%s: warning - non-zero self-VF (%.1f%%) for surface '%s'\n",
                                                progname, 100.*uvfa[3*m + 1],
                                                idf_getfield(pptr,NAME_FLD)->val);
                                continue;       /* don't record self-factor */
                        }
                        sprintf(uvfbuf, "%.4f", uvfa[3*m + 1]);
                        if (!idf_addfield(pout,
                                        idf_getfield(pptr,NAME_FLD)->val, NULL) ||
                                !idf_addfield(pout,
                                        idf_getfield(pptr1,NAME_FLD)->val, NULL) ||
                                !idf_addfield(pout, uvfbuf,
                                                (n || m < zp->nsurf-2) ?
                                                        "\n        " : "\n\n")) {
                                fputs(progname, stderr);
                                fputs(": error adding UVF fields\n", stderr);
                                return(0);
                        }
                }
                if (vfsum < 0.95)
                        fprintf(stderr,
                "%s: warning - missing %.1f%% of energy from surface '%s'\n",
                                        progname, 100.*(1.-vfsum),
                                        idf_getfield(pptr,NAME_FLD)->val);
        }
        free(uvfa);                             /* clean up and return */
        return(1);
}

/* Compute zone User View Factors */
static int
compute_zones(void)
{
        ZONE    *zptr;
                                                /* temporary octree name */
        mktemp(strcpy(temp_octree, TEMPLATE));
                                                /* compute each zone */
        for (zptr = zone_list; zptr != NULL; zptr = zptr->next) {
                SUBPROC rcproc;
                                                /* start rcontrib process */
                if (!start_rcontrib(&rcproc, zptr))
                        return(0);
                                                /* compute+add view factors */
                if (!compute_uvfs(&rcproc, zptr))
                        return(0);
                if (close_process(&rcproc) != 0) {
                        fputs(progname, stderr);
                        fputs(": bad return status from rcontrib\n", stderr);
                        return(0);
                }
        }
        unlink(temp_octree);                    /* remove octree file */
        return(1);
}

/* Load IDF and compute User View Factors */
int
main(int argc, char *argv[])
{
        int             incl_comments = 1;
        char            *origIDF, *revIDF;
        IDF_PARAMETER   *pptr;
        int             i;

        progname = argv[0];
        if (argc > 2 && !strcmp(argv[1], "-c")) {
                incl_comments = -1;             /* output header only */
                ++argv; --argc;
        }
        if ((argc < 2) | (argc > 3)) {
                fputs("Usage: ", stderr);
                fputs(progname, stderr);
                fputs(" [-c] Model.idf [Revised.idf]\n", stderr);
                return(1);
        }
        origIDF = argv[1];
        revIDF = (argc == 2) ? argv[1] : argv[2];
                                                /* load Input Data File */
        our_idf = idf_load(origIDF);
        if (our_idf == NULL) {
                fputs(progname, stderr);
                fputs(": cannot load IDF '", stderr);
                fputs(origIDF, stderr);
                fputs("'\n", stderr);
                return(1);
        }
                                                /* remove existing UVFs */
        if ((pptr = idf_getparam(our_idf, UVF_PNAME)) != NULL) {
                IDF_PARAMETER   *pnext;
                fputs(progname, stderr);
                fputs(": removing previous User View Factors\n", stderr);
                do {
                        pnext = pptr->pnext;
                        idf_delparam(our_idf, pptr);
                } while (pnext != NULL);
        }
                                                /* add to header */
        if (our_idf->hrem == NULL ||
                (i = strlen(our_idf->hrem)-strlen(ADD_HEADER)) < 0 ||
                        strcmp(our_idf->hrem+i, ADD_HEADER))
                idf_add2hdr(our_idf, ADD_HEADER);
                                                /* gather zone surfaces */
        for (i = 0; surf_type[i].pname != NULL; i++)
                for (pptr = idf_getparam(our_idf, surf_type[i].pname);
                                pptr != NULL; pptr = pptr->pnext) {
                        IDF_FIELD       *fptr = idf_getfield(pptr,
                                                        surf_type[i].zone_fld);
                        if (fptr == NULL) {
                                fputs(progname, stderr);
                                fputs(": warning - missing zone field\n", stderr);
                                continue;
                        }
                        if (add2zone(pptr, fptr->val) == NULL)
                                return(1);
                }
                                                /* run rcontrib on each zone */
        if (!compute_zones())
                return(1);
                                                /* write out modified IDF */
        if (!idf_write(our_idf, revIDF, incl_comments)) {
                fputs(progname, stderr);
                fputs(": error writing IDF '", stderr);
                fputs(revIDF, stderr);
                fputs("'\n", stderr);
                return(1);
        }
        return(0);                              /* finito! */
}
Revision:	2.6
Committed:	Tue Feb 11 22:01:44 2014 UTC (11 years, 8 months ago) by greg
Content type:	text/plain
Branch:	MAIN
Changes since 2.5:	+3 -3 lines
Log Message:	Checked calculation accuracy to better than 1%
#	User	Rev	Content
1	greg	2.1	#ifndef lint
2	greg	2.6	static const char RCSid[] = "$Id: eplus_adduvf.c,v 2.5 2014/02/11 21:17:29 greg Exp $";
3	greg	2.1	#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	greg	2.3	#include "rtio.h"
12			#include "rtmath.h"
13			#include "random.h"
14	greg	2.1	#include "eplus_idf.h"
15			#include "triangulate.h"
16	greg	2.2	#include "rtprocess.h"
17	greg	2.1
18			#ifndef NSAMPLES
19	greg	2.6	#define NSAMPLES 80000 /* number of samples to use */
20	greg	2.1	#endif
21
22	greg	2.2	char progname; / global argv[0] */
23
24	greg	2.5	char temp_octree[128]; /* temporary octree */
25	greg	2.2
26			const char UVF_PNAME[] =
27			"ZoneProperty:UserViewFactor:bySurfaceName";
28	greg	2.1
29			const char ADD_HEADER[] =
30	greg	2.5	"\n!+++ User View Factors computed by Radiance +++!\n\n";
31	greg	2.1
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	greg	2.5	{"Floor:Detailed", 3, 9},
43			{"RoofCeiling:Detailed", 3, 9},
44			{"Wall:Detailed", 3, 9},
45	greg	2.1	{NULL}
46			};
47
48	greg	2.2	typedef struct s_zone {
49	greg	2.1	const char zname; / zone name */
50			struct s_zone next; / next zone in list */
51	greg	2.2	int nsurf; /* surface count */
52	greg	2.1	IDF_PARAMETER pfirst; / first matching parameter */
53			IDF_PARAMETER plast; / last matching parameter */
54	greg	2.2	} ZONE; /* a list of collected zone surfaces */
55
56			ZONE zone_list = NULL; / our list of zones */
57	greg	2.1
58			IDF_LOADED our_idf = NULL; / loaded/modified IDF */
59
60			/* Create a new zone and push to top of our list */
61	greg	2.2	static ZONE *
62	greg	2.1	new_zone(const char zname, IDF_PARAMETER param)
63			{
64	greg	2.2	ZONE znew = (ZONE )malloc(sizeof(ZONE));
65	greg	2.1
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	greg	2.2	znew->nsurf = 1;
72	greg	2.1	return(zone_list = znew);
73			}
74
75			/* Add the detailed surface (polygon) to the named zone */
76	greg	2.2	static ZONE *
77	greg	2.1	add2zone(IDF_PARAMETER param, const char zname)
78			{
79	greg	2.2	ZONE *zptr;
80	greg	2.1
81			for (zptr = zone_list; zptr != NULL; zptr = zptr->next)
82			if (!strcmp(zptr->zname, zname))
83			break;
84			if (zptr == NULL)
85			return(new_zone(zname, param));
86			/* keep surfaces together */
87			if (!idf_movparam(our_idf, param, zptr->plast))
88			return(NULL);
89			zptr->plast = param;
90	greg	2.2	zptr->nsurf++;
91	greg	2.1	return(zptr);
92			}
93
94	greg	2.3	/* Return field for vertices in the given parameter */
95			static IDF_FIELD *
96			get_vlist(IDF_PARAMETER param, const char zname)
97	greg	2.2	{
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	greg	2.3	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	greg	2.2	/* return field for #verts */
111	greg	2.3	return(idf_getfield(param, surf_type[i].vert_fld));
112	greg	2.2	}
113
114			/* Convert surface to Radiance with modifier based on unique name */
115			static int
116	greg	2.3	rad_surface(IDF_PARAMETER param, FILE ofp)
117	greg	2.2	{
118	greg	2.3	const char *sname = idf_getfield(param,NAME_FLD)->val;
119			IDF_FIELD *fptr = get_vlist(param, NULL);
120	greg	2.2	int nvert, i;
121
122	greg	2.3	if (fptr == NULL \|\| (nvert = atoi(fptr->val)) < 3) {
123	greg	2.2	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	greg	2.3	if (fptr == NULL \|\| !isflt(fptr->val)) {
132	greg	2.2	fprintf(stderr,
133	greg	2.3	"%s: missing/bad vertex for %s '%s'\n",
134			progname, param->pname, sname);
135	greg	2.2	return(0);
136			}
137			fputc('\t', ofp);
138	greg	2.3	fputs(fptr->val, ofp);
139	greg	2.2	}
140			fputc('\n', ofp);
141			}
142			return(!ferror(ofp));
143			}
144
145	greg	2.3	/* Start rcontrib process */
146	greg	2.2	static int
147			start_rcontrib(SUBPROC pd, ZONE zp)
148			{
149	greg	2.5	#define BASE_AC 5
150	greg	2.2	static char *base_av[BASE_AC] = {
151	greg	2.5	"rcontrib", "-ff", "-h", "-x", "1"
152	greg	2.2	};
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+2zp->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	greg	2.3	IDF_FIELD *fptr = idf_getfield(pptr,NAME_FLD);
176			if (fptr == NULL \|\| !fptr->val[0]) {
177	greg	2.2	fputs(progname, stderr);
178			fputs(": missing name for surface parameter\n", stderr);
179			return(0);
180			}
181	greg	2.3	if (!rad_surface(pptr, ofp)) /* add surface to octree */
182			return(0);
183	greg	2.2	av[i++] = "-m";
184	greg	2.3	av[i++] = fptr->val;
185	greg	2.2	}
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	greg	2.3	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	greg	2.5	for (i = nv; i--; ) /* reverse vertex ordering */
228	greg	2.3	for (j = 0; j < 3; j++) {
229			if (fptr == NULL) {
230			fputs(progname, stderr);
231	greg	2.5	fputs(": missing vertex in init_poly()\n", stderr);
232	greg	2.3	return(NULL);
233			}
234			vl3[i][j] = atof(fptr->val);
235	greg	2.5	fptr = fptr->next;
236	greg	2.3	}
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	greg	2.5	ps->poff = DOT(vl3[0], ps->sdir[2]);
260	greg	2.3	/* 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	greg	2.5	(ps->poff+.001)*ps->sdir[2][i];
284	greg	2.3	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].mXvl2->v[b].mY - vl2->v[b].mX*vl2->v[a].mY +
291			vl2->v[b].mXvl2->v[c].mY - vl2->v[c].mXvl2->v[b].mY +
292			vl2->v[c].mXvl2->v[a].mY - vl2->v[a].mXvl2->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	greg	2.4	for (i = ns; i--; ) { /* stratified Monte Carlo sampling */
313	greg	2.3	double sv[4];
314	greg	2.5	FVECT dv;
315	greg	2.3	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	greg	2.5	samp[i6 + 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[i6 + 3 + j] = dv[0]ps->sdir[0][j] +
327			dv[1]*ps->sdir[1][j] +
328			dv[2]*ps->sdir[2][j] ;
329	greg	2.3	}
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	greg	2.5	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	greg	2.3	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	greg	2.2	/* Compute User View Factors using open rcontrib process */
362			static int
363	greg	2.3	compute_uvfs(SUBPROC pd, ZONE zp)
364	greg	2.2	{
365	greg	2.3	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	greg	2.4	double vfsum = 0;
390	greg	2.3	/* 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)3zp->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	greg	2.4	vfsum += uvfa[3*m + 1];
404	greg	2.3	if (pptr1 == pptr) {
405			if (uvfa[3*m + 1] > .001)
406			fprintf(stderr,
407	greg	2.5	"%s: warning - non-zero self-VF (%.1f%%) for surface '%s'\n",
408			progname, 100.uvfa[3m + 1],
409	greg	2.3	idf_getfield(pptr,NAME_FLD)->val);
410			continue; /* don't record self-factor */
411			}
412	greg	2.6	sprintf(uvfbuf, "%.4f", uvfa[3*m + 1]);
413	greg	2.3	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	greg	2.4	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	greg	2.3	}
431			free(uvfa); /* clean up and return */
432			return(1);
433	greg	2.2	}
434
435	greg	2.1	/* Compute zone User View Factors */
436			static int
437			compute_zones(void)
438			{
439	greg	2.2	ZONE *zptr;
440			/* temporary octree name */
441	greg	2.5	mktemp(strcpy(temp_octree, TEMPLATE));
442	greg	2.2	/* 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	greg	2.1	}
460
461			/* Load IDF and compute User View Factors */
462			int
463			main(int argc, char *argv[])
464			{
465	greg	2.2	int incl_comments = 1;
466	greg	2.1	char origIDF, revIDF;
467			IDF_PARAMETER *pptr;
468			int i;
469
470	greg	2.2	progname = argv[0];
471			if (argc > 2 && !strcmp(argv[1], "-c")) {
472			incl_comments = -1; /* output header only */
473			++argv; --argc;
474			}
475	greg	2.1	if ((argc < 2) \| (argc > 3)) {
476			fputs("Usage: ", stderr);
477	greg	2.2	fputs(progname, stderr);
478			fputs(" [-c] Model.idf [Revised.idf]\n", stderr);
479	greg	2.1	return(1);
480			}
481			origIDF = argv[1];
482			revIDF = (argc == 2) ? argv[1] : argv[2];
483			/* load Input Data File */
484			our_idf = idf_load(origIDF);
485			if (our_idf == NULL) {
486	greg	2.2	fputs(progname, stderr);
487	greg	2.1	fputs(": cannot load IDF '", stderr);
488			fputs(origIDF, stderr);
489			fputs("'\n", stderr);
490			return(1);
491			}
492			/* remove existing UVFs */
493			if ((pptr = idf_getparam(our_idf, UVF_PNAME)) != NULL) {
494			IDF_PARAMETER *pnext;
495	greg	2.2	fputs(progname, stderr);
496	greg	2.1	fputs(": removing previous User View Factors\n", stderr);
497			do {
498			pnext = pptr->pnext;
499			idf_delparam(our_idf, pptr);
500			} while (pnext != NULL);
501			}
502			/* add to header */
503			if (our_idf->hrem == NULL \|\|
504			(i = strlen(our_idf->hrem)-strlen(ADD_HEADER)) < 0 \|\|
505			strcmp(our_idf->hrem+i, ADD_HEADER))
506			idf_add2hdr(our_idf, ADD_HEADER);
507			/* gather zone surfaces */
508			for (i = 0; surf_type[i].pname != NULL; i++)
509			for (pptr = idf_getparam(our_idf, surf_type[i].pname);
510			pptr != NULL; pptr = pptr->pnext) {
511			IDF_FIELD *fptr = idf_getfield(pptr,
512			surf_type[i].zone_fld);
513			if (fptr == NULL) {
514	greg	2.2	fputs(progname, stderr);
515			fputs(": warning - missing zone field\n", stderr);
516	greg	2.1	continue;
517			}
518	greg	2.3	if (add2zone(pptr, fptr->val) == NULL)
519	greg	2.1	return(1);
520			}
521			/* run rcontrib on each zone */
522			if (!compute_zones())
523			return(1);
524			/* write out modified IDF */
525	greg	2.2	if (!idf_write(our_idf, revIDF, incl_comments)) {
526			fputs(progname, stderr);
527	greg	2.1	fputs(": error writing IDF '", stderr);
528			fputs(revIDF, stderr);
529			fputs("'\n", stderr);
530			return(1);
531			}
532			return(0); /* finito! */
533			}