Skip to content

GitLab

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

🚧 interpolation of materials - to be debugged

parent 983ae8ae
Expand all Hide whitespace changes
Inline Side-by-side
src/interface/beam.c
View file @ 196b1efe
......@@ -128,56 +128,8 @@ void initBeam(Beam *bm)
}
// Initialize the materials
// TODO - This need to be reworked to interpolate the mass matrices at the nkp positions
FILE *COMP = fopen("Compliance.dat","r");
if (COMP==NULL)
{printf("Unable to read compliance matrix file\n"); exit(EXIT_FAILURE);}
FILE *MASS = fopen("Mass.dat","r");
if ((MASS==NULL) && (bm->analysis_flag==2))
{printf("Unable to read mass matrix file for dynamic analysis\n"); exit(EXIT_FAILURE);}
beam_setMaterial(bm);
for (i=0;i<bm->nmate;i++)
{
if (i==0)
{
for (j=0; j<NSTRN; j++)
{
READL(line,len,COMP);
frmt = sscanf(line,"%lf %lf %lf %lf %lf %lf",&bm->material[i][j][0],&bm->material[i][j][1],&bm->material[i][j][2],&bm->material[i][j][3],&bm->material[i][j][4],&bm->material[i][j][5]); CHKFRMT(frmt,6);
}
if(bm->analysis_flag!=0){
for (j=NSTRN;j<NSTRN+NSTRN;j++){
READL(line,len,MASS);
frmt = sscanf(line,"%lf %lf %lf %lf %lf %lf",&bm->material[i][j][0],&bm->material[i][j][1],&bm->material[i][j][2],&bm->material[i][j][3],&bm->material[i][j][4],&bm->material[i][j][5]); CHKFRMT(frmt,6);
}
}
}
else
{
// COPY compliance matrix
for (j=0;j<NSTRN;j++)
{
for(k=0;k<NSTRN;k++)
{
bm->material[i][j][k] = bm->material[0][j][k];
}
}
// Copy mass matrix
if(bm->analysis_flag!=0){
for (j=NSTRN;j<NSTRN+NSTRN;j++){
for(k=0;k<NSTRN;k++)
{
bm->material[i][j][k] = bm->material[0][j][k];
}
}
}
}
}
fclose(COMP);
fclose(MASS);
for(i=0;i<bm->nframe;i++){
......@@ -883,3 +835,286 @@ void freeBeam(Beam *bm)
DEALLOCATE1(bm->pt_condition); //freePrescriInf
printf("Done!\n");
}
// Library utils
static void beam_interpolateMatrix(double **out, double x, double *xp, double ***matrix, int np, int nx, int ny, int ixStart, int iyStart)
{
int i,j;
int idx = 0;
double coeff;
for(i=0;i<np;i++)
{
if(x>xp[i])
{
idx++;
}
}
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]);}
for(i=0;i<nx;i++)
{
for(j=0;j<ny;j++)
{
out[i+ixStart][j+iyStart] = matrix[idx][i][j];
}
}
}
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]);
for(i=0;i<nx;i++)
{
for(j=0;j<ny;j++)
{
out[i+ixStart][j+iyStart] = matrix[idx-1][i][j];
}
}
}
else
{
coeff = (x-xp[idx-1])/(xp[idx]-xp[idx-1]);
coeff = (coeff>1.0 ? 1.0 : coeff);
for(i=0;i<nx;i++)
{
for(j=0;j<ny;j++)
{
out[i+ixStart][j+iyStart] = matrix[idx][i][j] * coeff + matrix[idx-1][i][j] * (1.0-coeff);
}
}
}
}
/*
static void beam_setMaterial(Beam *bm)
{
int i,j,k;
int ncomp, nmass;
double ***matComp, ***matMass;
double *xcomp, *xmass;
double x;
char* line = NULL;
ssize_t read;
int frmt;
size_t len = 0;
// Open the file containing the compliance matrices
FILE *COMP = fopen("Compliance.dat","r");
if (COMP==NULL)
{printf("Unable to read compliance matrix file\n"); exit(EXIT_FAILURE);}
// Open the file containing the mass matrices if dynamic analysis
FILE *MASS = fopen("Mass.dat","r");
if ((MASS==NULL) && (bm->analysis_flag!=0))
{printf("Unable to read mass matrix file for dynamic analysis\n"); exit(EXIT_FAILURE);}
// Read all the compliance matrices
READ1(COMP,ncomp,%d);
ALLOCATE3(matComp,double,ncomp,NSTRN,NSTRN);
ALLOCATE1(xcomp,double,ncomp);
for (i=0;i<ncomp;i++)
{
READ1(COMP,xcomp[i],%lf);
for (j=0; j<NSTRN; j++)
{
READL(line,len,COMP);
frmt = sscanf(line,"%lf %lf %lf %lf %lf %lf",&matComp[i][j][0],&matComp[i][j][1],&matComp[i][j][2],&matComp[i][j][3],&matComp[i][j][4],&matComp[i][j][5]); CHKFRMT(frmt,6);
}
}
// Read all the mass matrices
if ((bm->analysis_flag!=0))
{
READ1(MASS,nmass,%d);
ALLOCATE3(matMass,double,nmass,NSTRN,NSTRN);
ALLOCATE1(xmass,double,nmass);
for (i=0;i<nmass;i++)
{
READ1(MASS,xmass[i],%lf);
for (j=0; j<NSTRN; j++)
{
READL(line,len,MASS);
frmt = sscanf(line,"%lf %lf %lf %lf %lf %lf",&matMass[i][j][0],&matMass[i][j][1],&matMass[i][j][2],&matMass[i][j][3],&matMass[i][j][4],&matMass[i][j][5]); CHKFRMT(frmt,6);
}
}
}
// Interpolate the matrices at the point position
for (i=0;i<bm->nmate;i++)
{
// Warning: this requires that bm->nmate = bm->nkp!!
// Position are given in arc length / L
//x = sqrt(bm->coord[i][0]*bm->coord[i][0] + bm->coord[i][1]*bm->coord[i][1] + bm->coord[i][2]*bm->coord[i][2]) / bm->L;
// Interpolate the material at the correct point position
//beam_interpolateMatrix(bm->material[i],x,xcomp,matComp,ncomp,NSTRN,NSTRN,0 ,0);
for (j=0;j<NSTRN;j++){
for(k=0;k<NSTRN;k++)
{
bm->material[i][j][k] = matComp[0][j][k];
}
}
if (bm->analysis_flag!=0)
{
//beam_interpolateMatrix(bm->material[i],x,xmass,matMass,nmass,NSTRN,NSTRN,NSTRN,0);
// Copy mass matrix
for (j=NSTRN;j<NSTRN+NSTRN;j++){
for(k=0;k<NSTRN;k++)
{
bm->material[i][j][k] = matMass[0][j-NSTRN][k];
}
}
}
//for(j=0;j<bm->ndof_el-NSTRN;j++)
//{
// for(k=0;k<NSTRN;k++)
// {
// printf("%e ",bm->material[i][j][k]);
// }
// printf("\n");
//}
//printf("\n");
}
DEALLOCATE1(xcomp);
DEALLOCATE3(matComp,ncomp,NSTRN);
if (bm->analysis_flag!=0)
{
DEALLOCATE1(xmass);
DEALLOCATE3(matMass,ncomp,NSTRN);
}
fclose(COMP);
fclose(MASS);
}*/
static void beam_setMaterial(Beam *bm)
{
int i,j,k;
char* line = NULL;
ssize_t read;
int frmt;
size_t len = 0;
int ncomp, nmass;
double x;
FILE *MASS, *COMP;
double ***matComp, ***matMass;
double *xcomp, *xmass;
// TODO - This need to be reworked to interpolate the mass matrices at the nkp positions
COMP = fopen("Compliance.dat","r");
if (COMP==NULL)
{printf("Unable to read compliance matrix file\n"); exit(EXIT_FAILURE);}
MASS = fopen("Mass.dat","r");
if ((MASS==NULL) && (bm->analysis_flag==2))
{printf("Unable to read mass matrix file for dynamic analysis\n"); exit(EXIT_FAILURE);}
// Read all the compliance matrices
// READ1(COMP,ncomp,%d);
// ALLOCATE3(matComp,double,2,2,2);
// ALLOCATE1(xcomp,double,1);
//for (i=0;i<ncomp;i++)
//{
// READ1(COMP,xcomp[i],%lf);
// for (j=0; j<NSTRN; j++)
// {
// READL(line,len,COMP);
// frmt = sscanf(line,"%lf %lf %lf %lf %lf %lf",&matComp[i][j][0],&matComp[i][j][1],&matComp[i][j][2],&matComp[i][j][3],&matComp[i][j][4],&matComp[i][j][5]); CHKFRMT(frmt,6);
// }
//}
// Read all the mass matrices
//if ((bm->analysis_flag!=0))
//{
// READ1(MASS,nmass,%d);
// ALLOCATE3(matMass,double,2,2,2);
// ALLOCATE1(xmass,double,1);
// for (i=0;i<nmass;i++)
// {
// READ1(MASS,xmass[i],%lf);
// for (j=0; j<NSTRN; j++)
// {
// READL(line,len,MASS);
// frmt = sscanf(line,"%lf %lf %lf %lf %lf %lf",&matMass[i][j][0],&matMass[i][j][1],&matMass[i][j][2],&matMass[i][j][3],&matMass[i][j][4],&matMass[i][j][5]); CHKFRMT(frmt,6);
// }
// }
//}
//fclose(COMP);
//fclose(MASS);
COMP = fopen("Compliance.dat","r");
if (COMP==NULL)
{printf("Unable to read compliance matrix file\n"); exit(EXIT_FAILURE);}
MASS = fopen("Mass.dat","r");
if ((MASS==NULL) && (bm->analysis_flag==2))
{printf("Unable to read mass matrix file for dynamic analysis\n"); exit(EXIT_FAILURE);}
READ1(COMP,ncomp,%d);
READ1(COMP,x,%lf);
READ1(MASS,nmass,%d);
READ1(MASS,x,%lf);
for (i=0;i<bm->nmate;i++)
{
if (i==0)
{
for (j=0; j<NSTRN; j++)
{
READL(line,len,COMP);
frmt = sscanf(line,"%lf %lf %lf %lf %lf %lf",&bm->material[i][j][0],&bm->material[i][j][1],&bm->material[i][j][2],&bm->material[i][j][3],&bm->material[i][j][4],&bm->material[i][j][5]); CHKFRMT(frmt,6);
}
if(bm->analysis_flag!=0){
for (j=NSTRN;j<NSTRN+NSTRN;j++){
READL(line,len,MASS);
frmt = sscanf(line,"%lf %lf %lf %lf %lf %lf",&bm->material[i][j][0],&bm->material[i][j][1],&bm->material[i][j][2],&bm->material[i][j][3],&bm->material[i][j][4],&bm->material[i][j][5]); CHKFRMT(frmt,6);
}
}
}
else
{
// COPY compliance matrix
for (j=0;j<NSTRN;j++)
{
for(k=0;k<NSTRN;k++)
{
bm->material[i][j][k] = bm->material[0][j][k];
}
}
// Copy mass matrix
if(bm->analysis_flag!=0){
for (j=NSTRN;j<NSTRN+NSTRN;j++){
for(k=0;k<NSTRN;k++)
{
bm->material[i][j][k] = bm->material[0][j][k];
}
}
}
}
}
//DEALLOCATE1(xcomp);
//DEALLOCATE3(matComp,1,1);
//if (bm->analysis_flag!=0)
//{
// DEALLOCATE1(xmass);
// DEALLOCATE3(matMass,1,1);
//}
fclose(COMP);
fclose(MASS);
}
\ No newline at end of file
src/interface/beam.h
View file @ 196b1efe
......@@ -16,7 +16,7 @@
#define ALLOCATE2(ptr,type,n1,n2) ptr = (type **)malloc((n1) * sizeof(type *)); \
for (i_=0; i_<(n1); i_++){\
ptr[i_] = (type *)malloc((n2) * sizeof(type));}
#define ALLOCATE3(ptr,type,n1,n2,n3) ptr = (type ***)malloc(((n1)) * sizeof(type **)); \
#define ALLOCATE3(ptr,type,n1,n2,n3) ptr = (type ***)malloc((n1) * sizeof(type **)); \
for (i_=0; i_<(n1); i_++){\
ptr[i_] = (type **)malloc((n2) * sizeof(type *));\
for (j_=0; j_<(n2); j_++){\
......@@ -152,7 +152,8 @@ typedef struct
void initBeam(Beam *bm);
void beam_setTwist(Beam *bm, double *twist, double twist_per_L);
void beam_setTwist(Beam *bm, double *twist);
void beam_setTwistFromFile(Beam *bm, char *fname);
void beam_setLoads(Beam *bm, double *loads, int load_no);
void beam_setEndLoad(Beam *bm, double loads, int load_no);
......@@ -164,6 +165,10 @@ void beam_writeSolToFile(Beam *bm, char *fileName);
void freeBeam(Beam *bm);
// Private functions
static void beam_interpolateMatrix(double **out, double x, double *xp, double ***matrix, int np, int nx, int ny, int ixStart, int iyStart);
static void beam_setMaterial(Beam *bm);
void fortran_analysis(int nkp,int nelem,int ndof_el,int nmemb,int ncond_pt,int nmate,int nframe,int ndistrfun,int ncurv, double **coord,
int **member,struct PrescriInf **pt_condition, double **material, int niter, int nstep, double **sol_pt, double **sol_mb,char *error, int ncond_mb,
......
src/interface/beam_save.c 0 → 100644
View file @ 196b1efe
This diff is collapsed.
Supports Markdown
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