ViewVC Help

View File | Revision Log | Show Annotations | Download File | Root Listing

root/radiance/src/util/eplus_adduvf.c

Comparing src/util/eplus_adduvf.c (file contents):
Revision 2.2 by greg, Sun Feb 9 22:19:30 2014 UTC vs.
Revision 2.4 by greg, Mon Feb 10 17:50:02 2014 UTC

#	Line 7 | Line 7 | static const char RCSid[] = "$Id$";
7		*      G.Ward for LBNL
8		*/
9
10	–	#include <stdio.h>
10		#include <stdlib.h>
11	<	#include <string.h>
11	>	#include "rtio.h"
12	>	#include "rtmath.h"
13	>	#include "random.h"
14		#include "eplus_idf.h"
15		#include "triangulate.h"
16		#include "rtprocess.h"
#	Line 87 | Line 88 | add2zone(IDF_PARAMETER *param, const char *zname)
88		return(zptr);
89		}
90
91	<	/* Determine if a parameter is a surface in the indicated zone */
92	<	static int
93	<	in_zone(IDF_PARAMETER *param, const char *zname)
91	>	/* Return field for vertices in the given parameter */
92	>	static IDF_FIELD *
93	>	get_vlist(IDF_PARAMETER *param, const char *zname)
94		{
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	<	return(0);
101	<	/* get zone field */
102	<	fptr = idf_getfield(param, surf_type[i].zone_fld);
103	<	if (fptr == NULL)
104	<	return(0);
105	<	/* check for match */
106	<	if (strcmp(fptr->arg, zname))
106	<	return(0);
100	>	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		/* return field for #verts */
108	<	return(surf_type[i].vert_fld);
108	>	return(idf_getfield(param, surf_type[i].vert_fld));
109		}
110
111		/* Convert surface to Radiance with modifier based on unique name */
112		static int
113	<	rad_surface(const char *zname, IDF_PARAMETER *param, FILE *ofp)
113	>	rad_surface(IDF_PARAMETER *param, FILE *ofp)
114		{
115	<	const char      *sname = idf_getfield(param, NAME_FLD)->arg;
116	<	IDF_FIELD       *fptr = idf_getfield(param, in_zone(param, zname));
115	>	const char      *sname = idf_getfield(param,NAME_FLD)->val;
116	>	IDF_FIELD       *fptr = get_vlist(param, NULL);
117		int             nvert, i;
118
119	<	if (fptr == NULL || (nvert = atoi(fptr->arg)) < 3) {
119	>	if (fptr == NULL || (nvert = atoi(fptr->val)) < 3) {
120		fprintf(stderr, "%s: bad surface '%s'\n", progname, sname);
121		return(0);
122		}
#	Line 125 | Line 125 | rad_surface(const char *zname, IDF_PARAMETER *param, F
125		while (nvert--) {
126		for (i = 3; i--; ) {
127		fptr = fptr->next;
128	<	if (fptr == NULL) {
128	>	if (fptr == NULL || !isflt(fptr->val)) {
129		fprintf(stderr,
130	<	"%s: missing vertex fields in surface '%s'\n",
131	<	progname, sname);
130	>	"%s: missing/bad vertex for %s '%s'\n",
131	>	progname, param->pname, sname);
132		return(0);
133		}
134		fputc('\t', ofp);
135	<	fputs(fptr->arg, ofp);
135	>	fputs(fptr->val, ofp);
136		}
137		fputc('\n', ofp);
138		}
139		return(!ferror(ofp));
140		}
141
142	<	/* Strat rcontrib process */
142	>	/* Start rcontrib process */
143		static int
144		start_rcontrib(SUBPROC *pd, ZONE *zp)
145		{
#	Line 169 | Line 169 | start_rcontrib(SUBPROC *pd, ZONE *zp)
169		av[i++] = "-c";
170		av[i++] = cbuf;                         /* add modifier arguments */
171		for (n = zp->nsurf, pptr = zp->pfirst; n--; pptr = pptr->dnext) {
172	<	IDF_FIELD       *fptr = idf_getfield(pptr, NAME_FLD);
173	<	if (fptr == NULL) {
172	>	IDF_FIELD       *fptr = idf_getfield(pptr,NAME_FLD);
173	>	if (fptr == NULL || !fptr->val[0]) {
174		fputs(progname, stderr);
175		fputs(": missing name for surface parameter\n", stderr);
176		return(0);
177		}
178	<	av[i++] = "-m";
179	<	av[i++] = fptr->arg;
180	<	if (!rad_surface(zp->zname, pptr, ofp))
178	>	if (!rad_surface(pptr, ofp))    /* add surface to octree */
179		return(0);
180	+	av[i++] = "-m";
181	+	av[i++] = fptr->val;
182		}
183		if (pclose(ofp) != 0) {                 /* finish oconv */
184		fputs(progname, stderr);
#	Line 197 | Line 197 | start_rcontrib(SUBPROC *pd, ZONE *zp)
197		#undef BASE_AC
198		}
199
200	+	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].mX*vl2->v[b].mY - vl2->v[b].mX*vl2->v[a].mY +
287	+	vl2->v[b].mX*vl2->v[c].mY - vl2->v[c].mX*vl2->v[b].mY +
288	+	vl2->v[c].mX*vl2->v[a].mY - vl2->v[a].mX*vl2->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--; ) {           /* stratified Monte Carlo sampling */
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[ns*6 + j] = orig[j] + sv[0]*ab[j] + sv[1]*ac[j];
315	+	sv[3] = sqrt(sv[3]);
316	+	sv[4] *= 2.*PI;
317	+	samp[ns*6 + 3] = tcos(sv[4]) * sv[3];
318	+	samp[ns*6 + 4] = tsin(sv[4]) * sv[3];
319	+	samp[ns*6 + 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		/* Compute User View Factors using open rcontrib process */
355		static int
356	<	compute_uvfs(SUBPROC *pd, ZONE *sp)
356	>	compute_uvfs(SUBPROC *pd, ZONE *zp)
357		{
358	<
358	>	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	>	double  vfsum = 0;
383	>	/* send samples to rcontrib */
384	>	if (!sample_surface(pptr, pd->w))
385	>	return(0);
386	>	/* read results */
387	>	if (readbuf(pd->r, (char *)uvfa, sizeof(float)*3*zp->nsurf) !=
388	>	sizeof(float)*3*zp->nsurf) {
389	>	fputs(progname, stderr);
390	>	fputs(": read error from rcontrib process\n", stderr);
391	>	return(0);
392	>	}
393	>	/* append UVF fields */
394	>	for (m = 0, pptr1 = zp->pfirst;
395	>	m < zp->nsurf; m++, pptr1 = pptr1->dnext) {
396	>	vfsum += uvfa[3*m + 1];
397	>	if (pptr1 == pptr) {
398	>	if (uvfa[3*m + 1] > .001)
399	>	fprintf(stderr,
400	>	"%s: warning - non-zero self-VF (%.3e) for surface '%s'\n",
401	>	progname, uvfa[3*m + 1],
402	>	idf_getfield(pptr,NAME_FLD)->val);
403	>	continue;       /* don't record self-factor */
404	>	}
405	>	sprintf(uvfbuf, "%.6f", uvfa[3*m + 1]);
406	>	if (!idf_addfield(pout,
407	>	idf_getfield(pptr,NAME_FLD)->val, NULL) ||
408	>	!idf_addfield(pout,
409	>	idf_getfield(pptr1,NAME_FLD)->val, NULL) ||
410	>	!idf_addfield(pout, uvfbuf,
411	>	(n || m < zp->nsurf-2) ?
412	>	"\n        " : "\n\n")) {
413	>	fputs(progname, stderr);
414	>	fputs(": error adding UVF fields\n", stderr);
415	>	return(0);
416	>	}
417	>	}
418	>	if (vfsum < 0.95)
419	>	fprintf(stderr,
420	>	"%s: warning - missing %.1f%% of energy from surface '%s'\n",
421	>	progname, 100.*(1.-vfsum),
422	>	idf_getfield(pptr,NAME_FLD)->val);
423	>	}
424	>	free(uvfa);                             /* clean up and return */
425	>	return(1);
426		}
427
428		/* Compute zone User View Factors */
#	Line 291 | Line 512 | main(int argc, char *argv[])
512		fputs(": warning - missing zone field\n", stderr);
513		continue;
514		}
515	<	if (add2zone(pptr, fptr->arg) == NULL)
515	>	if (add2zone(pptr, fptr->val) == NULL)
516		return(1);
517		}
518		/* run rcontrib on each zone */

Diff Legend

–	Removed lines
+	Added lines
<	Changed lines (old)
>	Changed lines (new)