elamodel.c

su 的源代码库

	int n,i;
	float dx,dy,ddx,ddy,ds,
		*uin,*u,*x,*y,*tx,*ty,*c,
		(*xind)[4],(*yind)[4];
	
	/* if only one (x,z) input, do not smooth */
	if (nin==1) {
		*nout = nin;
		*xout = xin;
		*yout = yin;
		*txout = NULL;
		*tyout = NULL;
		*cout = NULL;
		return;
	}
	
	/* initially, copy input segments to output segments */
	n = nin;
	uin = alloc1float(nin);
	u = alloc1float(n);
	x = alloc1float(n);
	y = alloc1float(n);
	for (i=0; i<n; ++i) {
		x[i] = xin[i];
		y[i] = yin[i];
	}
	
	/* spline parameterization is chord length */
	uin[0] = u[0] = 0.0;
	for (i=1; i<n; ++i)
		uin[i] = u[i] = u[i-1]+
			sqrt(pow(x[i]-x[i-1],2.0)+pow(y[i]-y[i-1],2.0));
	
	/* compute cubic interpolation coefficients */
	xind = (float(*)[4])alloc1float(4*nin);
	yind = (float(*)[4])alloc1float(4*nin);
	csplin(nin,uin,xin,xind);
	csplin(nin,uin,yin,yind);
	
	/* loop over interior vertices */
	for (i=1; i<n-1; ++i)
		smoothTwoSegments(maxangle,i,nin,uin,xind,yind,&n,&u,&x,&y);
		
	/* tangent vectors, and curvatures */
	tx = alloc1float(n);
	ty = alloc1float(n);
	c = alloc1float(n);
	for (i=0; i<n; ++i) {
		intcub(1,nin,uin,xind,1,&u[i],&dx);
		intcub(1,nin,uin,yind,1,&u[i],&dy);
		intcub(2,nin,uin,xind,1,&u[i],&ddx);
		intcub(2,nin,uin,yind,1,&u[i],&ddy);
		ds = sqrt(dx*dx+dy*dy);
		tx[i] = dx/ds;
		ty[i] = dy/ds;
		c[i] = (dx*ddy-dy*ddx)/(ds*ds*ds);
	}
		
	/* free workspace */
	free1float(uin);
	free1float(u);
	free1float((float*)xind);
	free1float((float*)yind);
	
	/* set output parameters before returning */
	*nout = n;  *xout = x;  *yout = y;  ;
	*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 = realloc1float(uo,no);
	xo = realloc1float(xo,no);
	yo = realloc1float(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)
/* make a new vertex */
{
	Vertex *v;
	
	v = addVertexToModel(m,z,x);
	if (v==NULL) v = nearestVertexInModel(m,vlast,z,x);
	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*)malloc(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*)malloc(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*)malloc(sizeof(EdgeUseAttributes));
	eua->tx = tx2;
	eua->tz = tz2;
	eua->c = c2;
}

static void fillcoeff (Model *m, int nr, float **sr)
/* fill regions bounded by fixed edges with sloths */
{
	int ir,mindex;
	float x,z,a1111,a3333,a1133,a1313,a1113,a3313;
	float a1212,a2323,a1223,rho;
	

	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];
		a1111 = sr[ir][2];  a3333 = sr[ir][3];
		a1133 = sr[ir][4];  a1313 = sr[ir][5]; 
		a1113 = sr[ir][6];  a3313 = sr[ir][7];
		a1212 = sr[ir][8];  a2323 = sr[ir][9]; 
		a1223 = sr[ir][10]; rho = sr[ir][11];
		mindex = (int) sr[ir][12];

		

		/* determine triangle containing point (x,z) */
		t = insideTriInModel(m,NULL,z,x);
		
		/* flood triangles in region */
		setcoeff(t,a1111,a3333,a1133,a1313,a1113,a3313,a1212,a2323,
			a1223,rho,mindex);

	}
}

static void setcoeff (Tri *t, float a1111, float a3333, float a1133, float
 	a1313, float a1113, float a3313, float a1212, float a2323,
	float a1223, float rho, int mindex)
/* recursively set elastic coefficients in triangles */
{
	EdgeUse *eu;
	FaceAttributes *fa;
	
	/* if sloth already set, then return */
	if ((fa=t->fa)!=NULL)
		if (fa->a1111==a1111 && fa->a3333==a3333 && fa->a1133==a1133
		&& fa->a1313==a1313 && fa->a3313==a3313 && fa->a1113==a1113
		&& fa->a1212==a1212 && fa->a2323==a2323 && fa->a1223==a1223
		&& fa->rho==rho)
			return;
	
	/* if necessary, allocate space for attributes */
	if (fa==NULL)
		t->fa = fa = (FaceAttributes*)malloc(sizeof(FaceAttributes));
	
	/* set attributes */
	fa->a1111 = a1111;
	fa->a3333 = a3333;
	fa->a1133 = a1133;
	fa->a1313 = a1313;
	fa->a1113 = a1113;
	fa->a3313 = a3313;
	fa->a1212 = a1212;
	fa->a2323 = a2323;
	fa->a1223 = a1223;
	fa->rho   = rho;
	fa->mindex = mindex;

	/* for each edge not fixed, set attributes in opposite triangle */
	eu = t->eu;
	do {
		if (!eu->e->fixed)
		  	setcoeff(eu->euMate->f,a1111,a3333,a1133,a1313,a1113,
				a3313,a1212,a2323,a1223,rho,mindex);
		eu = eu->euCW;
	} while (eu!=t->eu);
}

static void checkfill (Model *m, int nr, float **sra)
/* check if all regions are filled */
{
	int ir, icheck;
	float fac;
	Face *f;
	FaceAttributes *fa;
	fa = (FaceAttributes*)malloc(sizeof(FaceAttributes)); 
 
        /* loop over faces */
   
	f = m->f;
	do {
		if(f->fa != NULL) {
			fa = f->fa;
		} else {
		  err("\n ERROR: You didn't fill all regions \n");
		}

		fac = fa->a1111;

		icheck = 0 ;
		for (ir=0; ir<nr; ++ir) {
			if(fac == sra[ir][2]){
				icheck = 1;
				break;
			}
		}
		if(!icheck)
		  err("\n ERROR: You didn't fill all regions \n");

		f = f->fNext;

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

}

int rottensors (float *aijkl, float angle)
/*****************************************************************************
Rotate 2-D elastic stiffness tensor
******************************************************************************
Input:
aijkl		input/output stiffness elements

Author:  Andreas Rueger, Colorado School of Mines, 03/14/94
******************************************************************************/
{
		float si,co,ss,cc,sc,a,c,f,l,m,n;
		float o,p,q;
		si=sin(angle); co=cos(angle);
		ss=si*si; cc=co*co;
		sc=si*co;
		a=aijkl[2];
		c=aijkl[3];
		f=aijkl[4];
		l=aijkl[5];
		m=aijkl[6];
		n=aijkl[7];
		o=aijkl[8];
		p=aijkl[9];
		q=aijkl[10];

		aijkl[2]=a*cc*cc+c*ss*ss+2*f*cc*ss+4*l*cc*ss-4*m*cc*sc-
			4*n*ss*sc;
		aijkl[3]=a*ss*ss+c*cc*cc+2*f*cc*ss+4*l*cc*ss+4*m*ss*sc+
			4*n*cc*sc;
		aijkl[4]=a*ss*cc+c*ss*cc+f*(ss*ss+cc*cc)-4*l*ss*cc-2*m*
			(ss*sc-cc*sc)-2*n*(cc*sc- ss*sc);
		aijkl[5]=a*ss*cc+c*ss*cc-2*f*ss*cc+l*(cc*cc+ss*ss-2*ss*cc)
			-2*m*(ss*sc-cc*sc)-2*n*(cc*sc- ss*sc);
		aijkl[6]=a*cc*sc-c*ss*sc+f*(ss*sc-cc*sc)+2*l*(ss*sc-cc*sc)+
			m*(cc*cc -3*cc*ss) +(3*cc*ss - ss*ss )*n;;
		aijkl[7]=a*ss*sc-c*cc*sc+f*(cc*sc-ss*sc)+2*l*(cc*sc-ss*sc)+
			 m*( 3*cc*ss - ss*ss ) + n*( cc*cc - 3*ss*cc);
		aijkl[8]=cc*o-2*sc*q+ss*p;
		aijkl[9]=ss*o+2*sc*q+cc*p;
		aijkl[10]=sc*o+(cc-ss)*q-sc*p;
		
		return 1;

}