src/util/eplus_adduvf.c

#ifndef lint
static const char RCSid[] = "$Id: eplus_adduvf.c,v 2.8 2014/02/12 00:12:38 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                   /* default number of samples */
#endif

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

int             nsamps = NSAMPLES;              /* number of samples to use */

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;                 /* object 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}
        };                              /* IDF surface types */

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

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

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

typedef struct {
        FVECT           sdir[3];        /* UVW unit sampling vectors */
        double          poff;           /* W-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 */

/* Create a new zone and push to top of our list */
static ZONE *
new_zone(const char *zname, IDF_OBJECT *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_OBJECT *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_movobject(our_idf, param, zptr->plast))
                return(NULL);
        zptr->plast = param;
        zptr->nsurf++;
        return(zptr);
}

/* Return field for vertices in the given object */
static IDF_FIELD *
get_vlist(IDF_OBJECT *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_OBJECT *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", "-fff", "-h", "-x", "1"
                        };
        char            cbuf[300];
        char            **av;
        FILE            *ofp;
        IDF_OBJECT      *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", nsamps);
        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 object\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) <= 0) {        /* 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
}

/* 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_OBJECT *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 = nsamps;       /* 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_OBJECT      *pptr, *pout, *pptr1;
        float           *uvfa;
        char            uvfbuf[24];
        int             n, m;
                                                /* create output object */
        pout = idf_newobject(our_idf, UVF_PNAME,
                        "    ! computed by Radiance\n        ", zp->plast);
        if (pout == NULL) {
                fputs(progname, stderr);
                fputs(": cannot create new IDF object\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_OBJECT      *pptr;
        int             i;

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