src/util/eplus_adduvf.c

#ifndef lint
static const char RCSid[] = "$Id: eplus_adduvf.c,v 2.2 2014/02/09 22:19:30 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        100000                  /* 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[] =
                        "!+++ User View Factors computed by Radiance +++!\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},
                {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         3
        static char     *base_av[BASE_AC] = {
                                "rcontrib", "-ff", "-h"
                        };
        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 = 0; i < nv; i++)
                for (j = 0; j < 3; j++) {
                        if (fptr == NULL) {
                                fputs(progname, stderr);
                                fputs(": bad vertex in init_poly()\n", stderr);
                                return(NULL);
                        }
                        vl3[i][j] = atof(fptr->val);
                }
                                        /* 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]) + FTINY;
                                        /* 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*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--; ) {
                double  sv[4];
                multisamp(sv, 4, (i+frandom())/(double)ns);
                sv[0] *= sv[1] = sqrt(sv[1]);
                sv[1] = 1. - sv[1];
                for (j = 3; j--; )
                        samp[ns*6 + j] = orig[j] + sv[0]*ab[j] + sv[1]*ac[j];
                sv[3] = sqrt(sv[3]);
                sv[4] *= 2.*PI;
                samp[ns*6 + 3] = tcos(sv[4]) * sv[3];
                samp[ns*6 + 4] = tsin(sv[4]) * sv[3];
                samp[ns*6 + 5] = sqrt(1. - sv[3]*sv[3]);
        }
                                        /* 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) {
                fprintf(stderr, "%s: error in %s '%s'\n",
                                progname, param->pname,
                                idf_getfield(param,NAME_FLD)->val);
                return(0);
        }
        vlist2 = init_poly(&psamp, fptr->next, nv);
        if (vlist2 == NULL)
                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) {
                                                /* 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) {
                        if (pptr1 == pptr) {
                                if (uvfa[3*m + 1] > .001)
                                        fprintf(stderr,
                "%s: warning - non-zero self-VF (%.3e) for surface '%s'\n",
                                                progname, uvfa[3*m + 1],
                                                idf_getfield(pptr,NAME_FLD)->val);
                                continue;       /* don't record self-factor */
                        }
                        sprintf(uvfbuf, "%.6f", 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);
                        }
                }
        }
        free(uvfa);                             /* clean up and return */
        return(1);
}

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