--- ray/src/util/eplus_adduvf.c	2014/02/09 22:19:30	2.2
+++ ray/src/util/eplus_adduvf.c	2014/02/10 04:51:26	2.3
@@ -1,5 +1,5 @@
 #ifndef lint
-static const char RCSid[] = "$Id: eplus_adduvf.c,v 2.2 2014/02/09 22:19:30 greg Exp $";
+static const char RCSid[] = "$Id: eplus_adduvf.c,v 2.3 2014/02/10 04:51:26 greg Exp $";
 #endif
 /*
  * Add User View Factors to EnergyPlus Input Data File
@@ -7,9 +7,10 @@ static const char RCSid[] = "$Id: eplus_adduvf.c,v 2.2
  *	G.Ward for LBNL
  */
 
-#include <stdio.h>
 #include <stdlib.h>
-#include <string.h>
+#include "rtio.h"
+#include "rtmath.h"
+#include "random.h"
 #include "eplus_idf.h"
 #include "triangulate.h"
 #include "rtprocess.h"
@@ -87,36 +88,35 @@ add2zone(IDF_PARAMETER *param, const char *zname)
 	return(zptr);
 }
 
-/* Determine if a parameter is a surface in the indicated zone */
-static int
-in_zone(IDF_PARAMETER *param, const char *zname)
+/* 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(0);
-						/* get zone field */
-	fptr = idf_getfield(param, surf_type[i].zone_fld);
-	if (fptr == NULL)
-		return(0);
-						/* check for match */
-	if (strcmp(fptr->arg, zname))
-		return(0);
+			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(surf_type[i].vert_fld);
+	return(idf_getfield(param, surf_type[i].vert_fld));
 }
 
 /* Convert surface to Radiance with modifier based on unique name */
 static int
-rad_surface(const char *zname, IDF_PARAMETER *param, FILE *ofp)
+rad_surface(IDF_PARAMETER *param, FILE *ofp)
 {
-	const char	*sname = idf_getfield(param, NAME_FLD)->arg;
-	IDF_FIELD	*fptr = idf_getfield(param, in_zone(param, zname));
+	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->arg)) < 3) {
+	if (fptr == NULL || (nvert = atoi(fptr->val)) < 3) {
 		fprintf(stderr, "%s: bad surface '%s'\n", progname, sname);
 		return(0);
 	}
@@ -125,21 +125,21 @@ rad_surface(const char *zname, IDF_PARAMETER *param, F
 	while (nvert--) {
 		for (i = 3; i--; ) {
 			fptr = fptr->next;
-			if (fptr == NULL) {
+			if (fptr == NULL || !isflt(fptr->val)) {
 				fprintf(stderr,
-				"%s: missing vertex fields in surface '%s'\n",
-						progname, sname);
+					"%s: missing/bad vertex for %s '%s'\n",
+						progname, param->pname, sname);
 				return(0);
 			}
 			fputc('\t', ofp);
-			fputs(fptr->arg, ofp);
+			fputs(fptr->val, ofp);
 		}
 		fputc('\n', ofp);
 	}
 	return(!ferror(ofp));
 }
 
-/* Strat rcontrib process */
+/* Start rcontrib process */
 static int
 start_rcontrib(SUBPROC *pd, ZONE *zp)
 {
@@ -169,16 +169,16 @@ start_rcontrib(SUBPROC *pd, ZONE *zp)
 	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) {
+		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);
 		}
-		av[i++] = "-m";
-		av[i++] = fptr->arg;
-		if (!rad_surface(zp->zname, pptr, ofp))
+		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);
@@ -197,11 +197,225 @@ start_rcontrib(SUBPROC *pd, ZONE *zp)
 #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 *sp)
+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 */
@@ -291,7 +505,7 @@ main(int argc, char *argv[])
 				fputs(": warning - missing zone field\n", stderr);
 				continue;
 			}
-			if (add2zone(pptr, fptr->arg) == NULL)
+			if (add2zone(pptr, fptr->val) == NULL)
 				return(1);
 		}
 						/* run rcontrib on each zone */