trimodel.c

来自「该程序主要用于三角网」· C语言 代码 · 共 788 行 · 第 1/2 页

                ds = sqrt(dx*dx+dy*dy);
                xou[i]=xs;
                you[i]=ys;
                agou[i]=90.-atan2(-dx,dy)*180/PI;
                fprintf(xzfp,"%f %f\n",xs,ys);
        }

        /*output interface (x,z,angle)*/
        fwrite(xou,sizeof(float),nrays,xfp);
        fwrite(you,sizeof(float),nrays,zfp);
        fwrite(agou,sizeof(float),nrays,afp);

        } /* end of if */


	/* free workspace */
	free1float(uin);
	free1float(u);
	free1float((float*)xind);
	free1float((float*)yind);
	free1float((float*)sind);

	/* set output parameters before returning */
	*nout = n;  *xout = x;  *yout = y;  *sout = s;
	*txout = tx;  *tyout = ty;  *cout = c;
}

static void smoothTwoSegments (float maxangle, int i,
	int nd, float ud[], float xd[][4], float yd[][4],
	int *n, float *u[], float *x[], float *y[])
/* used by smoothLineSegments to smooth just two adjacent line segments */
{
	int no,inew,j;
	float *uo,*xo,*yo,xm,ym,xp,yp,dot,ams,aps,cosa,angle;

	/* input/output arrays describing line segments */
	no = *n;  uo = *u;  xo = *x;  yo = *y;

	/* if at endpoint, simply return */
	if (i==0 || i==no-1) return;

	/* line segments joined at vertex i */
	xm = xo[i]-xo[i-1];  ym = yo[i]-yo[i-1];
	xp = xo[i+1]-xo[i];  yp = yo[i+1]-yo[i];

	/* angle between segments */
	dot = xm*xp+ym*yp;
	ams = xm*xm+ym*ym;
	aps = xp*xp+yp*yp;
	cosa = dot/sqrt(ams*aps);
	cosa = MAX(-1.0,MIN(1.0,cosa));
	angle = acos(cosa);

	/* if angle is small enough, simply return */
	if (angle<=maxangle) return;

	/* make room for new vertex */
	no++;
	uo = erealloc1float(uo,no);
	xo = erealloc1float(xo,no);
	yo = erealloc1float(yo,no);

	/* divide longest segment */
	inew = (ams>aps)?i:i+1;
	for (j=no-1; j>inew; --j) {
		uo[j] = uo[j-1];
		xo[j] = xo[j-1];
		yo[j] = yo[j-1];
	}
	uo[inew] = 0.5*(uo[inew-1]+uo[inew+1]);
	intcub(0,nd,ud,xd,1,&uo[inew],&xo[inew]);
	intcub(0,nd,ud,yd,1,&uo[inew],&yo[inew]);

	/* smooth line segments affected by new vertex */
	smoothTwoSegments(maxangle,inew,nd,ud,xd,yd,&no,&uo,&xo,&yo);
	smoothTwoSegments(maxangle,inew-1,nd,ud,xd,yd,&no,&uo,&xo,&yo);
	smoothTwoSegments(maxangle,inew+1,nd,ud,xd,yd,&no,&uo,&xo,&yo);

	/* set output parameters before returning */
	*n = no;  *u = uo;  *x = xo;  *y = yo;
}

static Vertex *newVertex (Model *m, Vertex *vlast,
	float x, float z, float s)
/* make a new vertex */
{
	Vertex *v;
	VertexAttributes *va;

	v = addVertexToModel(m,z,x);
	if (v==NULL) v = nearestVertexInModel(m,vlast,z,x);
	if (v->va==NULL)
		v->va = va = (VertexAttributes*)
			ealloc1(1,sizeof(VertexAttributes));
	else
		va = v->va;
	va->s = s;
	return v;
}

static void setEdgeAttributes (Vertex *v1, Vertex *v2, int k)
/* set edge attributes for fixed edge between two connected vertices */
{
	VertexUse *vu;
	EdgeUse *eu;
	EdgeAttributes *ea;

	/* determine edge use from v1 to v2 */
	vu = v1->vu;
	do {
		eu = vu->eu;
		if (eu->euMate->vu->v==v2) break;
		vu = vu->vuNext;
	} while (vu!=v1->vu);

	/* if v1 and v2 not connected, just return */
	if (eu->euMate->vu->v!=v2) return;

	/* set edge attributes for edge between v1 and v2 */
	eu->e->ea = ea = (EdgeAttributes*) ealloc1(1,sizeof(EdgeAttributes));
	ea->k = k;
}

static void setEdgeUseAttributes (Vertex *v1, Vertex *v2,
	float tx1, float tz1, float c1, float tx2, float tz2, float c2)
/* set edge-use attributes for fixed edge between two connected vertices */
{
	VertexUse *vu;
	EdgeUse *eu;
	EdgeUseAttributes *eua;

	/* determine edge use from v1 to v2 */
	vu = v1->vu;
	do {
		eu = vu->eu;
		if (eu->euMate->vu->v==v2) break;
		vu = vu->vuNext;
	} while (vu!=v1->vu);

	/* if v1 and v2 not connected, just return */
	if (eu->euMate->vu->v!=v2) return;

	/* set edge-use attributes for edge-use from v1 to v2 */
	eu->eua = eua = (EdgeUseAttributes*)
		ealloc1(1,sizeof(EdgeUseAttributes));
	eua->tx = tx1;
	eua->tz = tz1;
	eua->c = c1;

	/* set edge-use attributes for edge-use from v2 to v1 */
	eu = eu->euMate;
	eu->eua = eua = (EdgeUseAttributes*)
		ealloc1(1,sizeof(EdgeUseAttributes));
	eua->tx = tx2;
	eua->tz = tz2;
	eua->c = c2;
}

static void makeSlothForTri (Tri *t)
/* make sloth function for triangle */
{
	FaceAttributes *fa;
	VertexAttributes *va;
	double x1,z1,x2,z2,x3,z3,s1,s2,s3,
		x2mx1,z2mz1,s2ms1,x3mx1,z3mz1,s3ms1,
		a,b,det,s00,dsdx,dsdz;
	
	x1 = t->eu->vu->v->y;
	z1 = t->eu->vu->v->x;
	va = t->eu->vu->v->va;
	s1 = va->s;
	x2 = t->eu->euCW->vu->v->y;
	z2 = t->eu->euCW->vu->v->x;
	va = t->eu->euCW->vu->v->va;
	s2 = va->s;
	x3 = t->eu->euCCW->vu->v->y;
	z3 = t->eu->euCCW->vu->v->x;
	va = t->eu->euCCW->vu->v->va;
	s3 = va->s;
	x2mx1 = x2-x1;  z2mz1 = z2-z1;  s2ms1 = s2-s1;
	x3mx1 = x3-x1;  z3mz1 = z3-z1;  s3ms1 = s3-s1;
	det = z3mz1*x2mx1-z2mz1*x3mx1;
	a = z3mz1*s2ms1-z2mz1*s3ms1;
	b = s3ms1*x2mx1-s2ms1*x3mx1;
	dsdx = a/det;
	dsdz = b/det;
	s00 = s2-dsdx*x2-dsdz*z2;
	t->fa = fa = (FaceAttributes*) ealloc1(1,sizeof(FaceAttributes));
	fa->s00 = s00;  fa->dsdx = dsdx;  fa->dsdz = dsdz;
        fa->dens = FLT_MAX;  fa->qfac = FLT_MAX;
}

static void fillsloth (Model *m, int nr, float **sr)
/* fill regions bounded by fixed edges with sloths */
{
	int ir;
	float x,z,x0,z0,s00,dsdx,dsdz;
	Tri *t;

	/* loop over regions for which sloth function is specified */
	for (ir=0; ir<nr; ++ir) {

		/* determine parameters of sloth function */
		x = sr[ir][0];  z = sr[ir][1];
		x0 = sr[ir][2];  z0 = sr[ir][3];
		s00 = sr[ir][4];  dsdx = sr[ir][5];  dsdz = sr[ir][6];

		/* adjust v0 for x0 and z0 */
		s00 -= x0*dsdx+z0*dsdz;

		/* determine triangle containing point (x,z) */
		t = insideTriInModel(m,NULL,z,x);

		/* flood triangles in region */
		setsloth(t,s00,dsdx,dsdz);
	}
}

static void setsloth (Tri *t, float s00, float dsdx, float dsdz)
/* recursively set sloth functions in triangles */
{
	EdgeUse *eu;
	FaceAttributes *fa;

	/* if sloth already set, then return */
	if ((fa=t->fa)!=NULL)
		if (fa->s00==s00 && fa->dsdx==dsdx && fa->dsdz==dsdz)
			return;

	/* if necessary, allocate space for attributes */
	if (fa==NULL)
		t->fa = fa = (FaceAttributes*)
			ealloc1(1,sizeof(FaceAttributes));

	/* set attributes */
	fa->s00 = s00;
	fa->dsdx = dsdx;
	fa->dsdz = dsdz;

	/* for each edge not fixed, set attributes in opposite triangle */
	eu = t->eu;
	do {
		if (!eu->e->fixed) setsloth(eu->euMate->f,s00,dsdx,dsdz);
		eu = eu->euCW;
	} while (eu!=t->eu);
}

static void filldens (Model *m, int nd, float **dptr)
/* fill regions bounded by fixed edges with density */
{
	int id, filled;
	float x,z,dens;
	Tri *t;
	Face *f;
	FaceAttributes *fa;

	/* loop over regions for which a density is specified */
	for (id=0; id<nd; ++id) {

		/* determine parameters of density function */
		x = dptr[id][0];  z = dptr[id][1];
		dens = dptr[id][2];  

		/* determine triangle containing point (x,z) */
		t = insideTriInModel(m,NULL,z,x);

		/* flood triangles in region */
		setdens(t,dens);
	}

	/* loop over faces to check that density is specified everywhere */
	f = m->f;
	do {
		fa = f->fa;
		dens = fa->dens;

		for (id=0,filled=0; id<nd; ++id) {
			if (dens==dptr[id][2]) {
				filled = 1;
				break;
			}
		}
		if(!filled)
			err("WARNING!!\n"
				"CHECK THAT ALL REGIONS ARE "
				"FILLED WITH A DENSITY VALUE\n");

		/* next face */
		f = f->fNext;

	} while (f!=m->f);
}

static void setdens (Tri *t, float dens)
/* recursively set density values in triangles */
{
	EdgeUse *eu;
	FaceAttributes *fa;

	/* if dens already set, then return */
	if ((fa=t->fa)!=NULL)
		if (fa->dens==dens)
			return;

	/* if necessary, allocate space for attributes */
	if (fa==NULL)
		t->fa = fa = (FaceAttributes*)
			ealloc1(1,sizeof(FaceAttributes));

	/* set attribute */
	fa->dens = dens;

	/* for each edge not fixed, set attributes in opposite triangle */
	eu = t->eu;
	do {
		if (!eu->e->fixed) setdens(eu->euMate->f,dens);
		eu = eu->euCW;
	} while (eu!=t->eu);
}

static void fillq (Model *m, int nq, float **qptr)
/* fill regions bounded by fixed edges with a Q-factor */
{
	int iq, filled;
	float x,z,qfac;
	Tri *t;
	Face *f;
	FaceAttributes *fa;

	/* loop over regions for which Q-factor is specified */
	for (iq=0; iq<nq; ++iq) {

		/* determine parameters of Q-function */
		x = qptr[iq][0];  z = qptr[iq][1];
		qfac = qptr[iq][2];  

		/* determine triangle containing point (x,z) */
		t = insideTriInModel(m,NULL,z,x);

		/* flood triangles in region */
		setq(t,qfac);
	}

	/* loop over faces to check that a Q-factor is specified everywhere */
	f = m->f;
	do {
		fa = f->fa;
		qfac = fa->qfac;

		for (iq=0,filled=0; iq<nq; ++iq) {
			if (qfac==qptr[iq][2]) {
				filled = 1;
				break;
			}
		}
		if(!filled)
			err("WARNING!!\n"
				"CHECK THAT ALL REGIONS ARE "
				"FILLED WITH A Q-FACTOR\n");

		/* next face */
		f = f->fNext;

	} while (f!=m->f);

}

static void setq (Tri *t, float qfac)
/* recursively set a Q-factor in triangles */
{
	EdgeUse *eu;
	FaceAttributes *fa;

	/* if Q-factor already set, then return */
	if ((fa=t->fa)!=NULL)
		if (fa->qfac==qfac)
			return;

	/* if necessary, allocate space for attributes */
	if (fa==NULL)
		t->fa = fa = (FaceAttributes*)
			ealloc1(1,sizeof(FaceAttributes));

	/* set attribute */
	fa->qfac = qfac;

	/* for each edge not fixed, set attributes in opposite triangle */
	eu = t->eu;
	do {
		if (!eu->e->fixed) setq(eu->euMate->f,qfac);
		eu = eu->euCW;
	} while (eu!=t->eu);
}