Commit 260f1400 authored by François Trigaux's avatar François Trigaux
Browse files

added setTwistFromFile

parent 0c0abb51
......@@ -107,7 +107,7 @@ void initBeam(Beam *bm)
bm->member[i][3] = i+2; // number of material 2
bm->member[i][4] = i+1; // frame number
bm->member[i][5] = 1; // number of division per member
bm->member[i][6] = i+1; // initial curvature number
bm->member[i][6] = 0; // initial curvature number -- NOTE: Default curvature must be changed to non-zero parameter only if twist is explicitely added!!
}
// Set the clamped-free conditions
......@@ -555,10 +555,12 @@ void initBeamFromGebtFile(Beam *bm)
// Assign the twist to each of the member by changing its frame
// twist_per_L [rad/arc length] is only required if there are subdivisions to the members
void beam_setTwist(Beam *bm, double *twist, double twist_per_L)
void beam_setTwist(Beam *bm, double *twist)
{
int i;
double k1;
for(i=0;i<bm->nframe;i++){
printf("%f\n",twist[i]);
bm->frame[i][0][0] = 1.0;
bm->frame[i][0][1] = 0.0;
bm->frame[i][0][2] = 0.0;
......@@ -569,8 +571,88 @@ void beam_setTwist(Beam *bm, double *twist, double twist_per_L)
bm->frame[i][2][1] = sin(twist[i]);
bm->frame[i][2][2] = cos(twist[i]);
bm->curvature[i][0] = twist_per_L;
k1 = (twist[i+1] - twist[i])/(bm->L / ((double)bm->nmemb));
if (fabs(k1)>1e-8)
{
bm->member[i][6] = i+1; // Associate a non-zero curvature to the member -- TODO: find why it does not work when curvature is zero
bm->curvature[i][0] = k1;
}
printf("%f\n",k1);
}
}
// Assign the twist to each of the member by changing its frame
// twist_per_L [rad/arc length] is only required if there are subdivisions to the members
void beam_setTwistFromFile(Beam *bm, char *fname)
{
int i,j,k;
int idx,np;
double coeff;
double x,*xp,*twistp,*twist;
char* line = NULL;
ssize_t read;
int frmt;
size_t len = 0;
// Open the file containing the mass matrices if dynamic analysis
FILE *fid = fopen("Twist.dat","r");
if (fid==NULL)
{printf("Unable to read mass matrix file for dynamic analysis\n"); exit(EXIT_FAILURE);}
// fill the vectors for interpolation
READ1(fid,np,%d);
ALLOCATE1(xp,double,np);
ALLOCATE1(twistp,double,np);
for(k=0;k<np;k++)
{
READL(line,len,fid);
frmt = sscanf(line,"%lf %lf",&xp[k],&twistp[k]); CHKFRMT(frmt,2);
}
ALLOCATE1(twist,double,bm->nmemb);
for (k=0;k<bm->nkp;k++)
{
//x = ( sqrt(bm->coord[k][0]*bm->coord[k][0] + bm->coord[k][1]*bm->coord[k][1] + bm->coord[k][2]*bm->coord[k][2]) + sqrt(bm->coord[k+1][0]*bm->coord[k+1][0] + bm->coord[k+1][1]*bm->coord[k+1][1] + bm->coord[k+1][2]*bm->coord[k+1][2]) ) / 2.0; //member coordinates
x = ( sqrt(bm->coord[k][0]*bm->coord[k][0] + bm->coord[k][1]*bm->coord[k][1] + bm->coord[k][2]*bm->coord[k][2]));
idx = 0;
for(i=0;i<np;i++)
{
if(x>xp[i])
idx++;
}
if (idx<=0)
{
if (np>1) printf("Warning: in Twist Interpolation : value of xp=%f is below or equal to xp[0]=%f \n",x,xp[0]);
twist[k] = twistp[idx];
}
else if (idx>=np)
{
if (np>1) printf("Warning: in Twist Interpolation : value of xp=%f is higher than xp[np-1]=%f \n",x,xp[np-1]);
twist[k] = twistp[idx-1];
}
else
{
coeff = (x-xp[idx-1])/(xp[idx]-xp[idx-1]);
coeff = (coeff>1.0 ? 1.0 : coeff);
twist[k] = twistp[idx]* coeff + twistp[idx-1] * (1.0-coeff);
}
twist[k] *= M_PI/180.0; // convert to radian
}
beam_setTwist(bm,twist); // NOTE: 0.0 is hardcoded
DEALLOCATE1(twist);
DEALLOCATE1(twistp);
DEALLOCATE1(xp);
fclose(fid);
}
void beam_writeSolToFile(Beam *bm,char* fileName)
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment