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

Added exchange functions

parent 1805e914
...@@ -43,6 +43,9 @@ $(OBJF90): $(OBJDIR)/%.o: $(DIRF90)/%.f90 ...@@ -43,6 +43,9 @@ $(OBJF90): $(OBJDIR)/%.o: $(DIRF90)/%.f90
$(OBJF): $(OBJDIR)/%.o: $(DIRF)/%.f $(OBJF): $(OBJDIR)/%.o: $(DIRF)/%.f
$(CF) $(FFLAGS) -fallow-argument-mismatch -c -o $@ $< -J$(OBJDIR) $(CF) $(FFLAGS) -fallow-argument-mismatch -c -o $@ $< -J$(OBJDIR)
install:
ln -f lib/libgebt.so /usr/local/lib
ln -f include/gebt.h /usr/local/include
clean: clean:
rm -vf $(OBJDIR)/*.o rm -vf $(OBJDIR)/*.o
......
...@@ -44,6 +44,9 @@ void initGebt(Gebt *bm) ...@@ -44,6 +44,9 @@ void initGebt(Gebt *bm)
SET(&(bm->dt),%lf,fileGebt,section,"dt"); SET(&(bm->dt),%lf,fileGebt,section,"dt");
} }
bm->t = 0.0;
bm->lastStep = -1;
bm->nmemb = bm->nkp-1; // specific for beams. bm->nmemb = bm->nkp-1; // specific for beams.
bm->nmate = bm->nkp; // 1 material per point bm->nmate = bm->nkp; // 1 material per point
bm->ncond_mb = bm->nmemb; // each member has one condition bm->ncond_mb = bm->nmemb; // each member has one condition
...@@ -184,7 +187,7 @@ void initGebt(Gebt *bm) ...@@ -184,7 +187,7 @@ void initGebt(Gebt *bm)
} }
} }
// This must be initialize but has no consequence as no time-functions are used // This is used to compute dt in Fortran. MUST BE UPDATED IF nstep OR dt ARE CHANGED!!!
bm->simu_time[0] = 0.0; bm->simu_time[0] = 0.0;
bm->simu_time[1] = bm->nstep * bm->dt; bm->simu_time[1] = bm->nstep * bm->dt;
...@@ -574,12 +577,12 @@ void gebt_setTwistFromFile(Gebt *bm, char *fname) ...@@ -574,12 +577,12 @@ void gebt_setTwistFromFile(Gebt *bm, char *fname)
if (idx<=0) 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]); //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]; twist[k] = twistp[idx];
} }
else if (idx>=np) 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]); //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]; twist[k] = twistp[idx-1];
} }
else else
...@@ -604,28 +607,57 @@ void gebt_setTwistFromFile(Gebt *bm, char *fname) ...@@ -604,28 +607,57 @@ void gebt_setTwistFromFile(Gebt *bm, char *fname)
void gebt_getNodeVelocities(double **vel, Gebt *bm){ void gebt_getNodeVelocities(double **vel, Gebt *bm){
int i,j; int i,j;
int step = bm->nstep-1; int step = bm->lastStep;
if(step<1)
{
for(i=0;i<bm->nkp;i++)
{
for(j=0;j<NSTRN;j++)
{
// !!! it is vel(j,i) to ease the access from interpolation methods!
vel[j][i] = 0.0;
}
}
for(i=0;i<bm->nkp;i++) }
else
{ {
for(j=0;j<NSTRN;j++) for(i=0;i<bm->nkp;i++)
{ {
// !!! it is vel(j,i) to ease the access from interpolation methods! for(j=0;j<NSTRN;j++)
vel[j][i] = (bm->sol_pt[step][i][j+NDIM] - bm->sol_pt[step-1][i][j+NDIM]) / bm->dt; {
// !!! it is vel(j,i) to ease the access from interpolation methods!
vel[j][i] = (bm->sol_pt[step][i][j+NDIM] - bm->sol_pt[step-1][i][j+NDIM]) / bm->dt;
}
} }
} }
} }
void gebt_getMemberVelocities(double **vel, Gebt *bm){ void gebt_getMemberVelocities(double **vel, Gebt *bm){
int i,j; int i,j;
int step = bm->nstep-1; int step = bm->lastStep;
for(i=0;i<bm->nelem;i++) if(step<1)
{ {
for(j=0;j<NSTRN;j++) for(i=0;i<bm->nelem;i++)
{ {
// !!! it is vel(j,i) to ease the access from interpolation methods! for(j=0;j<NSTRN;j++)
vel[j][i] = (bm->sol_mb[step][i][j+NDIM] - bm->sol_mb[step-1][i][j+NDIM]) / bm->dt; {
// !!! it is vel(j,i) to ease the access from interpolation methods!
vel[j][i] = bm->init_cond[i][j+NSTRN];
}
}
}
else
{
for(i=0;i<bm->nelem;i++)
{
for(j=0;j<NSTRN;j++)
{
// !!! it is vel(j,i) to ease the access from interpolation methods!
vel[j][i] = (bm->sol_mb[step][i][j+NDIM] - bm->sol_mb[step-1][i][j+NDIM]) / bm->dt;
}
} }
} }
} }
...@@ -635,16 +667,29 @@ void gebt_getMemberVelocities(double **vel, Gebt *bm){ ...@@ -635,16 +667,29 @@ void gebt_getMemberVelocities(double **vel, Gebt *bm){
*/ */
void gebt_getNodeDisplacement(double **disp, Gebt *bm){ void gebt_getNodeDisplacement(double **disp, Gebt *bm){
int i,j; int i,j;
int step = bm->nstep-1; int step = bm->lastStep;
for(i=0;i<bm->nkp;i++) if(step<0)
{ {
for(j=0;j<NSTRN;j++) for (i = 0; i < bm->nkp; i++)
{ {
// !!! it is disp(j,i) to ease the access from interpolation methods! for (j = 0; j < NSTRN; j++)
disp[j][i] = bm->sol_pt[step][i][j+NDIM]; {
disp[j][i] = 0.0;
}
} }
} }
else
{
for (i = 0; i < bm->nkp; i++)
{
for (j = 0; j < NSTRN; j++)
{
// !!! it is disp(j,i) to ease the access from interpolation methods!
disp[j][i] = bm->sol_pt[step][i][j + NDIM];
}
}
}
} }
/* disp = allocated double pointer of size (6,nelem) /* disp = allocated double pointer of size (6,nelem)
...@@ -652,53 +697,94 @@ void gebt_getNodeDisplacement(double **disp, Gebt *bm){ ...@@ -652,53 +697,94 @@ void gebt_getNodeDisplacement(double **disp, Gebt *bm){
*/ */
void gebt_getMemberDisplacement(double **disp, Gebt *bm){ void gebt_getMemberDisplacement(double **disp, Gebt *bm){
int i,j; int i,j;
int step = bm->nstep-1; int step = bm->lastStep;
for(i=0;i<bm->nelem;i++) if(step<0)
{
for(i=0;i<bm->nelem;i++)
{
for(j=0;j<NSTRN;j++)
{
// !!! it is vel(j,i) to ease the access from interpolation methods!
disp[j][i] = bm->init_cond[i][j];
}
}
}
else
{ {
for(j=0;j<NSTRN;j++) for(i=0;i<bm->nelem;i++)
{ {
// !!! it is disp(j,i) to ease the access from interpolation methods! for(j=0;j<NSTRN;j++)
disp[j][i] = bm->sol_mb[step][i][j+NDIM]; {
// !!! it is disp(j,i) to ease the access from interpolation methods!
disp[j][i] = bm->sol_mb[step][i][j+NDIM];
}
} }
} }
} }
/* Get the rotation matrix from the internal Wiener-Milenkovic parameters of GEBT.
*/
void gebt_getRotationMatrix(double out[3][3], double c[3])
{
double c0,fac;
c0 = 2.0 - (c[0]*c[0] + c[1]*c[1] + c[2]*c[2])/8.0;
fac = 1.0/(4.0-c0)/(4.0-c0);
out[0][0] = fac * (c0*c0 + c[0]*c[0] - c[1]*c[1] - c[2]*c[2]);
out[0][1] = fac * (2*(c[0]*c[1] + c0*c[2]));
out[0][2] = fac * (2*(c[0]*c[2]-c0*c[1]));
out[1][0] = fac * (2*(c[0]*c[1] - c0*c[2]));
out[1][1] = fac * (c0*c0 - c[0]*c[0] + c[1]*c[1] - c[2]*c[2]);
out[1][2] = fac * (2*(c[1]*c[2] + c0*c[0]));
out[2][0] = fac * (2*(c[0]*c[2] + c0*c[1]));
out[2][1] = fac * (2*(c[1]*c[2] - c0*c[0]));
out[2][2] = fac * (c0*c0 - c[0]*c[0] - c[1]*c[1] + c[2]*c[2]);
}
void gebt_writeSolToFile(Gebt *bm,char* fileName) void gebt_writeSolToFile(Gebt *bm,char* fileName)
{ {
int i,j,k; int i,j,k;
char fname[32]; char fname[32];
FILE *fid; FILE *fid;
for (k=0;k<bm->nstep;k++){ if(bm->lastStep<0)
{
printf("Unable to write the solution: No time-step have been performed with bm_analysis");
exit(EXIT_FAILURE);
}
sprintf(fname,"%s%03d.dat",fileName,k); //for (k=0;k<bm->nstep;k++){
fid = fopen(fname,"w"); k = bm->lastStep;
if (fid==NULL) sprintf(fname,"%s.dat",fileName);
{ fid = fopen(fname,"w");
printf("Unable to open file %s to write the solution of the structural analysis\n",fileName); if (fid==NULL)
exit(EXIT_FAILURE); {
} printf("Unable to open file %s to write the solution of the structural analysis\n",fileName);
exit(EXIT_FAILURE);
}
for (i=0;i<bm->nkp;i++) fprintf(fid,"Time = %f\n",bm->t);
for (i=0;i<bm->nkp;i++)
{
for(j=0;j<NDIM + NDOF_ND;j++)
{ {
for(j=0;j<NDIM + NDOF_ND;j++) fprintf(fid,"%1.8e ",bm->sol_pt[k][i][j]);
{
fprintf(fid,"%1.8e ",bm->sol_pt[k][i][j]);
}
fprintf(fid,"\n");
} }
fprintf(fid,"\n");
}
for (i=0;i<bm->nelem;i++) for (i=0;i<bm->nelem;i++)
{
for(j=0;j<NDIM + bm->ndof_el;j++)
{ {
for(j=0;j<NDIM + bm->ndof_el;j++) fprintf(fid,"%1.8e ",bm->sol_mb[k][i][j]);
{
fprintf(fid,"%1.8e ",bm->sol_mb[k][i][j]);
}
fprintf(fid,"\n");
} }
fclose(fid); fprintf(fid,"\n");
} }
fclose(fid);
//}
} }
...@@ -770,6 +856,10 @@ void gebt_analysis(Gebt *bm) ...@@ -770,6 +856,10 @@ void gebt_analysis(Gebt *bm)
material_tmp[i] = rav_material[i]; material_tmp[i] = rav_material[i];
} }
// In case dt or nstep have been changed, must recompute the simu-time because fortran uses it to compute the dt.
bm->simu_time[0] = 0.0;
bm->simu_time[1] = bm->nstep * bm->dt;
// Call the fortran function // Call the fortran function
fortran_analysis(bm->nkp,bm->nelem,bm->ndof_el,bm->nmemb,bm->ncond_pt,bm->nmate, bm->nframe,bm->ndistrfun, fortran_analysis(bm->nkp,bm->nelem,bm->ndof_el,bm->nmemb,bm->ncond_pt,bm->nmate, bm->nframe,bm->ndistrfun,
bm->ncurv,&rav_coord,&rav_member,&(bm->pt_condition),&material_tmp,bm->niter,bm->nstep, bm->ncurv,&rav_coord,&rav_member,&(bm->pt_condition),&material_tmp,bm->niter,bm->nstep,
...@@ -793,6 +883,8 @@ void gebt_analysis(Gebt *bm) ...@@ -793,6 +883,8 @@ void gebt_analysis(Gebt *bm)
// Set initialize to zero to perform dynamic analysis without redoing the initialization step // Set initialize to zero to perform dynamic analysis without redoing the initialization step
bm->initialize = 0; bm->initialize = 0;
bm->lastStep = bm->nstep-1;
bm->t += bm->nstep * bm->dt;
free(rav_eigen_val); free(rav_eigen_val);
free(rav_eigen_vec_mb); free(rav_eigen_vec_mb);
...@@ -918,7 +1010,7 @@ static void gebt_interpolateMatrix(double **out, double x, double *xp, double ** ...@@ -918,7 +1010,7 @@ static void gebt_interpolateMatrix(double **out, double x, double *xp, double **
if (idx==0) if (idx==0)
{ {
if (np>1){ printf("Warning: in Matrix Interpolation : value of xp=%f is below or equal to xp[0]=%f \n",x,xp[0]);} //if (np>1){ printf("Warning: in Matrix Interpolation : value of xp=%f is below or equal to xp[0]=%f \n",x,xp[0]);}
for(i=0;i<nx;i++) for(i=0;i<nx;i++)
{ {
for(j=0;j<ny;j++) for(j=0;j<ny;j++)
...@@ -929,7 +1021,7 @@ static void gebt_interpolateMatrix(double **out, double x, double *xp, double ** ...@@ -929,7 +1021,7 @@ static void gebt_interpolateMatrix(double **out, double x, double *xp, double **
} }
else if (idx==np) else if (idx==np)
{ {
if (np>1) printf("Warning: in Matrix Interpolation : value of xp=%f is higher than xp[np-1]=%f \n",x,xp[np-1]); //if (np>1) printf("Warning: in Matrix Interpolation : value of xp=%f is higher than xp[np-1]=%f \n",x,xp[np-1]);
for(i=0;i<nx;i++) for(i=0;i<nx;i++)
{ {
for(j=0;j<ny;j++) for(j=0;j<ny;j++)
......
...@@ -144,8 +144,9 @@ typedef struct ...@@ -144,8 +144,9 @@ typedef struct
// Personal add // Personal add
double L; // Beam Length double L; // Beam Length
int initFromFile; int initFromFile;
double dt; double dt,t;
int initialize; // flag indicating if an initialisation step is required to find the compatible state variables int initialize; // flag indicating if an initialisation step is required to find the compatible state variables
int lastStep; // number of the last solved step
} Gebt; } Gebt;
...@@ -166,6 +167,8 @@ void gebt_getMemberVelocities(double **vel, Gebt *bm); ...@@ -166,6 +167,8 @@ void gebt_getMemberVelocities(double **vel, Gebt *bm);
void gebt_getNodeDisplacement(double **disp, Gebt *bm); void gebt_getNodeDisplacement(double **disp, Gebt *bm);
void gebt_getMemberDisplacement(double **disp, Gebt *bm); void gebt_getMemberDisplacement(double **disp, Gebt *bm);
void gebt_getRotationMatrix(double out[3][3],double c[3]);
void gebt_writeGebtToFile(Gebt *bm,char *fname); void gebt_writeGebtToFile(Gebt *bm,char *fname);
void gebt_writeSolToFile(Gebt *bm, char *fileName); void gebt_writeSolToFile(Gebt *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