gebt.c

#include <stdlib.h>
#include <stdio.h>
#include <math.h>

#include "readINIFile.h"
#include "gebt.h"


int check;
int i_,j_,k_; // macro iterators

void initGebt(Gebt *bm)
{
    int i,j,k;
    char *section;
    dataINI fileGebt = readINIFile("Gebt.ini");

    char*   line = NULL; // don't forget to free! 
    ssize_t read;
    int     frmt;
    size_t  len = 0;

    /* Set the mesh parameters */ 
    /*section = "Beam";
    SET_OR_DEFAULT(&(bm->analysis_flag),%d,fileGebt,section,"analysis_flag",2);
    SET_OR_DEFAULT(&(bm->nkp),%d,fileGebt,section,"nkp",2);
    SET_OR_DEFAULT(&(bm->nframe),%d,fileGebt,section,"nframe",0);// number of different cross-sectional frames
    SET_OR_DEFAULT(&(bm->ncurv),%d,fileGebt,section,"ncurv",0);// number of sets of initial curvature/twists
    SET_OR_DEFAULT(&(bm->L),%lf,fileGebt,section,"L",1.0);// number of sets of initial curvature/twists
    SET_OR_DEFAULT(&(bm->initFromFile),%d,fileGebt,section,"initFromFile",0);// number of sets of initial curvature/twists

    */
    section = "Beam";
    SET_OR_DEFAULT(&(bm->analysis_flag),%d,fileGebt,section,"analysis_flag",2);
    SET_OR_DEFAULT(&(bm->nkp),%d,fileGebt,section,"nkp",2);

    SET_OR_DEFAULT(&(bm->L),%lf,fileGebt,section,"L",1.0);// number of sets of initial curvature/twists
    SET_OR_DEFAULT(&(bm->initFromFile),%d,fileGebt,section,"initFromFile",0);// number of sets of initial curvature/twists

    SET_OR_DEFAULT(&(bm->niter),%d,fileGebt,section,"niter",100);
    SET_OR_DEFAULT(&(bm->nstep),%d,fileGebt,section,"nstep",1);

    if(bm->analysis_flag>0){
        SET(&(bm->dt),%lf,fileGebt,section,"dt");
    }

    bm->t = 0.0;
    bm->lastStep = -1;

    bm->nmemb = bm->nkp-1; // specific for beams. 
    bm->nmate = bm->nkp;   // 1 material per point
    bm->ncond_mb  = bm->nmemb; // each member has one condition
    bm->ndistrfun = bm->nmemb; // each member will have its own distributed load.
    bm->nframe    = bm->nmemb; // each member will have its own frame.
    bm->ncurv     = bm->nmemb;
    bm->ncond_pt = 2; // condition for clamped-free beam

    bm->ntimefun = 0;

    bm->initialize = 1;

    if(bm->analysis_flag==0){
    bm->ndof_el=NDOF_ND;
    for (i = 0; i < NDIM; i++) 
    {
      bm->omega_a0[i]    = 0.0;
      bm->omega_a_tf[i]  = 0;
      bm->v_root_a0[i]   = 0.0;
      bm->v_root_a_tf[i] = 0;
    }
    }
    else{
        bm->ndof_el=18;

        for (i = 0; i < NDIM; i++) 
        {
            bm->omega_a_tf[i]  = 0;
            bm->v_root_a_tf[i] = 0;
        }
        SET_OR_DEFAULT(&bm->omega_a0[0],%lf,fileGebt,section,"omegaX",0.0);
        SET_OR_DEFAULT(&bm->omega_a0[1],%lf,fileGebt,section,"omegaY",0.0);
        SET_OR_DEFAULT(&bm->omega_a0[2],%lf,fileGebt,section,"omegaZ",0.0);
        SET_OR_DEFAULT(&bm->v_root_a0[0],%lf,fileGebt,section,"vrootX",0.0);
        SET_OR_DEFAULT(&bm->v_root_a0[1],%lf,fileGebt,section,"vrootY",0.0);
        SET_OR_DEFAULT(&bm->v_root_a0[2],%lf,fileGebt,section,"vrootZ",0.0);
    }

    freeDataINI(fileGebt);
    bm->nev = 0;

    ALLOCATE2(bm->coord, double, bm->nkp, NDIM);
    ALLOCATE2(bm->member, int, bm->nmemb, MEMB_CONST);
    ALLOCATE1(bm->pt_condition, struct PrescriInf, bm->ncond_pt);
    ALLOCATE3(bm->material, double, bm->nmate, bm->ndof_el - NSTRN, NSTRN);
    ALLOCATE3(bm->frame, double, bm->nframe, NDIM, NDIM);
    ALLOCATE1(bm->mb_condition, struct PrescriInf, bm->ncond_mb);
    ALLOCATE2(bm->distr_fun, double, bm->ndistrfun, NSTRN);
    ALLOCATE2(bm->curvature, double, bm->ncurv, NDIM);
    //ALLOCATE1(bm->time_function, struct TimeFunction, bm-> ntimefun);

    for (i = 0; i < bm->nkp; i++) {
        bm->coord[i][0] = i * bm->L/((double)bm->nmemb); 
        bm->coord[i][1] = 0.0;
        bm->coord[i][2] = 0.0;
    }

    for (i=0;i<bm->nmemb;i++){ 
        bm->member[i][0] = i+1; // number of point 1
        bm->member[i][1] = i+2; // number of point 2
        bm->member[i][2] = i+1; // number of material 1
        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] = 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
    bm->pt_condition[0].id = 1;       // first point
    bm->pt_condition[1].id = bm->nkp; // last point
    for (i=0;i<6;i++){
        bm->pt_condition[0].dof[i] = i+1;
        bm->pt_condition[0].value[i] = 0.0;
        bm->pt_condition[0].time_fun_no[i] = 0;
        bm->pt_condition[0].follower[i] = 0;
        bm->pt_condition[0].value_current[i] = bm->pt_condition[0].value[i];

        bm->pt_condition[1].dof[i] = i+1+6;
        bm->pt_condition[1].value[i] = 0.0;
        bm->pt_condition[1].time_fun_no[i] = 0;
        bm->pt_condition[1].follower[i] = 0;
        bm->pt_condition[1].value_current[i] = bm->pt_condition[1].value[i];
    }

    // Initialize the materials
    gebt_setMaterial(bm);
    


    for(i=0;i<bm->nframe;i++){
        for(j=0;j<NDIM;j++){
            for(k=0;k<NDIM;k++){
                bm->frame[i][j][k]= (j==k) ? 1.0 : 0.0;
            }
        }
    }

    // Initialize the loads to zero. Use proper function later to set the distributed loads
    if(bm->ncond_mb>0){
        struct PrescriInf *tmp;
        for (i=0;i<bm->ncond_mb;i++)
        {
            tmp = &(bm->mb_condition[i]);
            tmp->id = i+1;
                
            for (j=0;j<NSTRN;j++)
            {
                tmp->dof[j] = i+1; // Set the number of the distributed function
                tmp->value[j] = 0.0; // Set the value to multiply the distributed function with
                tmp->time_fun_no[j] = 0;
                tmp->follower[j] = 0; 
                tmp->value_current[j] = tmp->value[j];
            }
        }
    }


    // Use the proper function later to set the distributed load functions
    if(bm->ndistrfun>0){
    	for(i=0;i<bm->ndistrfun;i++)
        {
            for (j=0;j<NSTRN;j++)
            {
                bm->distr_fun[i][j] = 0.0;
            }
        }
    }


    if(bm->ncurv>0)
    {    
        for(i=0;i<bm->ncurv;i++)
        {
            for(j=0;j<NDIM;j++)
            {
                bm->curvature[i][j] = 0.0;
            }
        }
    }

    // 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[1] = bm->nstep * bm->dt;

    ALLOCATE1(bm->time_function,struct TimeFunction,bm->ntimefun);
    if (bm->ntimefun>0)
    {
        printf("Unable to read timefunctions from C\n");
        exit(EXIT_FAILURE);
    }


    bm->nelem=0;
    for (i=0;i<bm->nmemb;i++)
    {
        bm->nelem += bm->member[i][5]; // total number of elements for the structure
    }

    // Initialize conditions from file
    ALLOCATE2(bm->init_cond    , double             , bm->nelem    , 12);
    for (i=0;i<bm->nelem;i++)
    {
        for (j=0;j<2*NSTRN;j++)
        {
            bm->init_cond[i][j] = 0.0;
        }
    }

    if((bm->analysis_flag==2) && (bm->initFromFile==1))
    {
        // READ INITIAL CONIDITIONS FROM THE FILE
        FILE *INIT = fopen("Gebt_initCond.ini","r"); 
        if (INIT==NULL)
        {
            printf("Cannot read %s\n","Gebt_initCond.ini");
            exit(EXIT_FAILURE);
        }

        for (i=0;i<bm->nelem;i++)
        {
            READL(line,len,INIT);
            frmt = sscanf(line,"%lf %lf %lf %lf %lf %lf",&bm->init_cond[i][0],&bm->init_cond[i][1],&bm->init_cond[i][2],&bm->init_cond[i][3],&bm->init_cond[i][4],&bm->init_cond[i][5]); CHKFRMT(frmt,6);
        }

        for (i=0;i<bm->nelem;i++)
        {
            READL(line,len,INIT);
            frmt = sscanf(line,"%lf %lf %lf %lf %lf %lf",&bm->init_cond[i][6],&bm->init_cond[i][7],&bm->init_cond[i][8],&bm->init_cond[i][9],&bm->init_cond[i][10],&bm->init_cond[i][11]); CHKFRMT(frmt,6);
        }
        fclose(INIT);
    }


    ALLOCATE3(bm->sol_pt       , double             , bm->nstep    , bm->nkp  , NDIM + NDOF_ND);    //for each point there are 12 variables
    ALLOCATE3(bm->sol_mb       , double             , bm->nstep    , bm->nelem, NDIM + bm->ndof_el);  // for each element, there are 18 variables



    free(line);

}




void initGebtFromGebtFile(Gebt *bm)
{
  int i,j;
  char *section, dataFileName[32];
  dataINI fileGebt = readINIFile("Gebt.ini");

  int tmp_no;
  char *  line = NULL;
  ssize_t read;
  int     frmt;
  size_t  len = 0;
  double  a, b, c;

  FILE *IN = fopen("Gebt_data.dat","r");
  if (IN==NULL)
  {
      printf("Unable to read file %s for beam data\n","Gebt_data.dat");
      exit(EXIT_FAILURE);
  }

  /* Set the mesh parameters */ 
  section = "Beam";
  SET_OR_DEFAULT(&(bm->analysis_flag),%d,fileGebt,section,"analysis_flag",2);
  SET_OR_DEFAULT(&(bm->nkp),%d,fileGebt,section,"nkp",2);
  SET_OR_DEFAULT(&(bm->nmemb),%d,fileGebt,section,"nmemb",1);// number of beam members
  SET_OR_DEFAULT(&(bm->ncond_pt),%d,fileGebt,section,"ncond_pt",2);// number of point conditions
  SET_OR_DEFAULT(&(bm->nmate),%d,fileGebt,section,"nmate",1);// number of different cross-sections
  SET_OR_DEFAULT(&(bm->nframe),%d,fileGebt,section,"nframe",0);// number of different cross-sectional frames
  SET_OR_DEFAULT(&(bm->ndistrfun),%d,fileGebt,section,"ndistrfun",0);// number of distribution functions
  SET_OR_DEFAULT(&(bm->ncurv),%d,fileGebt,section,"ncurv",0);// number of sets of initial curvature/twists

  SET_OR_DEFAULT(&(bm->niter),%d,fileGebt,section,"niter",100);
  SET_OR_DEFAULT(&(bm->nstep),%d,fileGebt,section,"nstep",0);

  SET_OR_DEFAULT(&(bm->ncond_mb),%d,fileGebt,section,"ncond_mb",2);

  bm->ntimefun = 0; // impossible to pass time-functions to fortran

  if(bm->analysis_flag==0){
    bm->ndof_el=NDOF_ND;
    for (i = 0; i < NDIM; i++) 
    {
      bm->omega_a0[i]    = 0.0;
      bm->omega_a_tf[i]  = 0;
      bm->v_root_a0[i]   = 0.0;
      bm->v_root_a_tf[i] = 0;
    }
    }
    else{
        bm->ndof_el=18;

        READ3(IN,bm->omega_a0,%lf);
        READ3(IN,bm->omega_a_tf,%d);
        READ3(IN,bm->v_root_a0,%lf);
        READ3(IN,bm->v_root_a_tf,%d);

        for (i=0;i>NDIM;i++)
        {
            if ((bm->omega_a_tf[i] != 0) || (bm->v_root_a_tf[i] !=0))
            {
                printf("You can not set time functions for omega and v_root when using the C interface");
                exit(EXIT_FAILURE);
            }
        }

        if(bm->analysis_flag==3){
            READ1(IN,bm->nev,%d);
        }
        else
        {
            bm->nev = 0;
        }
    }


  ALLOCATE2(bm->coord        , double             , bm->nkp      , NDIM);
  ALLOCATE2(bm->member       , int                , bm->nmemb    , MEMB_CONST);
  ALLOCATE1(bm->pt_condition , struct PrescriInf  , bm->ncond_pt);
  ALLOCATE3(bm->material     , double             , bm->nmate   , bm->ndof_el - NSTRN, NSTRN);
  ALLOCATE3(bm->frame        , double             , bm->nframe  , NDIM               , NDIM);
  ALLOCATE1(bm->mb_condition , struct PrescriInf  , bm->ncond_mb);
  ALLOCATE2(bm->distr_fun    , double             , bm->ndistrfun, NSTRN);
  ALLOCATE2(bm->curvature    , double             , bm->ncurv    , NDIM);
  //ALLOCATE1(bm->time_function, struct TimeFunction, bm-> ntimefun);

  
  // Read the value of each variable
  for (i = 0; i < bm->nkp; i++) {
    READL(line, len, IN);
    frmt                 = sscanf(line, "%d %lf %lf %lf", &tmp_no, &a, &b, &c); CHKFRMT(frmt,4)
    bm->coord[tmp_no - 1][0] = a; 
    bm->coord[tmp_no - 1][1] = b;
    bm->coord[tmp_no - 1][2] = c;
  }

  for (i=0;i<bm->nmemb;i++){  
        READL(line,len,IN);
        frmt = sscanf(line,"%d",&tmp_no); CHKFRMT(frmt,1)
        frmt = sscanf(line,"%*d %d %d %d %d %d %d",&bm->member[tmp_no-1][0],&bm->member[tmp_no-1][1],&bm->member[tmp_no-1][2],&bm->member[tmp_no-1][3],&bm->member[tmp_no-1][4],&bm->member[tmp_no-1][5],&bm->member[tmp_no-1][6]); CHKFRMT(frmt,6);

    }

    for (i=0;i<bm->ncond_pt;i++){
        READL(line,len,IN);
        sscanf(line,"%d",&(bm->pt_condition[i].id));
        READL(line,len,IN);
        sscanf(line,"%d %d %d %d %d %d",&(bm->pt_condition[i].dof[0]),&(bm->pt_condition[i].dof[1]),&(bm->pt_condition[i].dof[2]),&(bm->pt_condition[i].dof[3]),&(bm->pt_condition[i].dof[4]),&(bm->pt_condition[i].dof[5]));
        READL(line,len,IN);
        sscanf(line,"%lf %lf %lf %lf %lf %lf",&(bm->pt_condition[i].value[0]),&(bm->pt_condition[i].value[1]),&(bm->pt_condition[i].value[2]),&(bm->pt_condition[i].value[3]),&(bm->pt_condition[i].value[4]),&(bm->pt_condition[i].value[5]));
        READL(line,len,IN);
        sscanf(line,"%d %d %d %d %d %d",&(bm->pt_condition[i].time_fun_no[0]),&(bm->pt_condition[i].time_fun_no[1]),&(bm->pt_condition[i].time_fun_no[2]),&(bm->pt_condition[i].time_fun_no[3]),&(bm->pt_condition[i].time_fun_no[4]),&(bm->pt_condition[i].time_fun_no[5]));
        READL(line,len,IN);
        sscanf(line,"%d %d %d %d %d %d",&(bm->pt_condition[i].follower[0]),&(bm->pt_condition[i].follower[1]),&(bm->pt_condition[i].follower[2]),&(bm->pt_condition[i].follower[3]),&(bm->pt_condition[i].follower[4]),&(bm->pt_condition[i].follower[5]));

        for(j=0;j<6;j++)
        {
            bm->pt_condition[i].value_current[j] = bm->pt_condition[i].value[j]; //the current value is initialized to be the original functional value, if not time varing it will remain the same. If time varying, it will be updated.
        }

    }

    for (i=0;i<bm->nmate;i++)
    {
        READ1(IN,tmp_no,%d);

        for (j=0; j<NSTRN; j++)
        {
            READL(line,len,IN);
            frmt = sscanf(line,"%lf %lf %lf %lf %lf %lf",&bm->material[tmp_no-1][j][0],&bm->material[tmp_no-1][j][1],&bm->material[tmp_no-1][j][2],&bm->material[tmp_no-1][j][3],&bm->material[tmp_no-1][j][4],&bm->material[tmp_no-1][j][5]); CHKFRMT(frmt,6);
        }
    

    if(bm->analysis_flag!=0){
            for (j=NSTRN;j<NSTRN+NSTRN;j++){
                READL(line,len,IN);
                frmt = sscanf(line,"%lf %lf %lf %lf %lf %lf",&bm->material[tmp_no-1][j][0],&bm->material[tmp_no-1][j][1],&bm->material[tmp_no-1][j][2],&bm->material[tmp_no-1][j][3],&bm->material[tmp_no-1][j][4],&bm->material[tmp_no-1][j][5]); CHKFRMT(frmt,6);  
            }   
        }
    }



    for(i=0;i<bm->nframe;i++){

        READ1(IN,tmp_no,%d);
        
        for(j=0;j<NDIM;j++){
            READL(line,len,IN);
            frmt = sscanf(line,"%lf %lf %lf",&bm->frame[tmp_no-1][j][0],&bm->frame[tmp_no-1][j][1],&bm->frame[tmp_no-1][j][2]); CHKFRMT(frmt,3);
            
        }
    }

    if(bm->ncond_mb>0){
        struct PrescriInf *tmp;
        for (i=0;i<bm->ncond_mb;i++)
        {
            tmp = &(bm->mb_condition[i]);
            READL(line,len,IN);
            frmt = sscanf(line,"%d",&tmp->id); CHKFRMT(frmt,1);
                
            READL(line,len,IN);
            frmt = sscanf(line,"%d %d %d %d %d %d",&tmp->dof[0],&tmp->dof[1],&tmp->dof[2],&tmp->dof[3],&tmp->dof[4],&tmp->dof[5]);

            READL(line,len,IN);
            frmt = sscanf(line,"%lf %lf %lf %lf %lf %lf",&tmp->value[0],&tmp->value[1],&tmp->value[2],&tmp->value[3],&tmp->value[4],&tmp->value[5]);

            READL(line,len,IN);
            frmt = sscanf(line,"%d %d %d %d %d %d",&tmp->time_fun_no[0],&tmp->time_fun_no[1],&tmp->time_fun_no[2],&tmp->time_fun_no[3],&tmp->time_fun_no[4],&tmp->time_fun_no[5]);
            
            READL(line,len,IN);
            frmt = sscanf(line,"%d %d %d %d %d %d",&tmp->follower[0],&tmp->follower[1],&tmp->follower[2],&tmp->follower[3],&tmp->follower[4],&tmp->follower[5]);

            for (j=0;j<NSTRN;j++)
            {
                tmp->value_current[j] = tmp->value[j];
            }
        }
    }

    if(bm->ndistrfun>0){
    	for(i=0;i<bm->ndistrfun;i++)
        {
            READ1(IN,tmp_no,%d);

            READL(line,len,IN);
            frmt = sscanf(line,"%lf %lf %lf %lf %lf %lf",&bm->distr_fun[tmp_no-1][0],&bm->distr_fun[tmp_no-1][1],&bm->distr_fun[tmp_no-1][2],&bm->distr_fun[tmp_no-1][3],&bm->distr_fun[tmp_no-1][4],&bm->distr_fun[tmp_no-1][5]); CHKFRMT(frmt,6);
        }

    }


    if(bm->ncurv>0){
        
        for(i=0;i<bm->ncurv;i++)
        {
            READ1(IN,tmp_no,%d);
            READ3(IN,bm->curvature[tmp_no-1],%lf);
        }
    }

    //if(bm->analysis_flag==2){
    	READL(line,len,IN);
        frmt = sscanf(line,"%lf %lf",&bm->simu_time[0],&bm->simu_time[1]); CHKFRMT(frmt,2);
    //}

    ALLOCATE1(bm->time_function,struct TimeFunction,bm->ntimefun);
    if (bm->ntimefun>0)
    {
        printf("Unable to read timefunctions from C\n");
        exit(EXIT_FAILURE);
    }


    bm->nelem=0;
    for (i=0;i<bm->nmemb;i++)
    {
        bm->nelem += bm->member[i][5]; // total number of elements for the structure
    }  

    ALLOCATE2(bm->init_cond    , double             , bm->nelem    , 12);
    if(bm->analysis_flag==2){

        // READ INITIAL CONIDITIONS FROM THE FILE
        FILE *INIT = fopen("Gebt_initCond.ini","r"); 
        if (INIT==NULL)
        {
            printf("Cannot read %s\n","Gebt_initCond.ini");
            exit(EXIT_FAILURE);
        }

        for (i=0;i<bm->nelem;i++)
        {
            READL(line,len,INIT);
            frmt = sscanf(line,"%lf %lf %lf %lf %lf %lf",&bm->init_cond[i][0],&bm->init_cond[i][1],&bm->init_cond[i][2],&bm->init_cond[i][3],&bm->init_cond[i][4],&bm->init_cond[i][5]); CHKFRMT(frmt,6);
        }

        for (i=0;i<bm->nelem;i++)
        {
            READL(line,len,INIT);
            frmt = sscanf(line,"%lf %lf %lf %lf %lf %lf",&bm->init_cond[i][6],&bm->init_cond[i][7],&bm->init_cond[i][8],&bm->init_cond[i][9],&bm->init_cond[i][10],&bm->init_cond[i][11]); CHKFRMT(frmt,6);
        }
        fclose(INIT);


    }


    ALLOCATE3(bm->sol_pt       , double             , bm->nstep    , bm->nkp  , NDIM + NDOF_ND);    //for each point there are 12 variables
    ALLOCATE3(bm->sol_mb       , double             , bm->nstep    , bm->nelem, NDIM + bm->ndof_el);  // for each element, there are 18 variables

}

// 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 gebt_setTwist(Gebt *bm, double *twist)
{
    int i;
    double k1;
    for(i=0;i<bm->nframe;i++){
        bm->frame[i][0][0] = 1.0;
        bm->frame[i][0][1] = 0.0;
        bm->frame[i][0][2] = 0.0;
        bm->frame[i][1][0] = 0.0;
        bm->frame[i][1][1] = cos(twist[i]);
        bm->frame[i][1][2] = -sin(twist[i]);
        bm->frame[i][2][0] = 0.0;
        bm->frame[i][2][1] = sin(twist[i]);
        bm->frame[i][2][2] = cos(twist[i]);

        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;
        }
    }
}

// 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 gebt_setTwistFromFile(Gebt *bm, char *fname)
{
    int i,j,k;
    int idx,np;
    double coeff;
    double x,*xp,*twistp,*twist;

    char*   line = NULL; //don't forget to free!
    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 twist file to set the twist for the blade\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->nkp);
    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
    }

    gebt_setTwist(bm,twist); // NOTE: 0.0 is hardcoded

    DEALLOCATE1(twist);
    DEALLOCATE1(twistp);
    DEALLOCATE1(xp);

    free(line);

    fclose(fid);
}

void gebt_getNodeVelocities(double **vel, Gebt *bm){
    int i,j;
    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;
            }
        }

    }
    else
    {
        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] = (bm->sol_pt[step][i][j+NDIM] - bm->sol_pt[step-1][i][j+NDIM]) / bm->dt;
            }
        }
    }
}

void gebt_getMemberVelocities(double **vel, Gebt *bm){
    int i,j;
    int step = bm->lastStep;

    if(step<1)
    {
        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->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;
            }
        }
    }
}

/* disp = allocated double pointer of size (6,nkp)
   /!\ the order of is reversed!
*/ 
void gebt_getNodeDisplacement(double **disp, Gebt *bm){
    int i,j;
    int step = bm->lastStep;

    if(step<0)
    {
        for (i = 0; i < bm->nkp; i++) 
        {
          for (j = 0; j < NSTRN; j++) 
          {
            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)
   /!\ the order of is reversed!
*/ 
void gebt_getMemberDisplacement(double **disp, Gebt *bm){
    int i,j;
    int step = bm->lastStep;

    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(i=0;i<bm->nelem;i++)
        {
            for(j=0;j<NSTRN;j++)
            {
                // !!! 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)
{
    int i,j,k;
    char fname[32];
    FILE *fid;

    if(bm->lastStep<0)
    {
        printf("Unable to write the solution: No time-step have been performed with bm_analysis");
        exit(EXIT_FAILURE);
    }

    //for (k=0;k<bm->nstep;k++){
    k = bm->lastStep;
    sprintf(fname,"%s.dat",fileName);
    fid = fopen(fname,"w");
    if (fid==NULL)
    {
        printf("Unable to open file %s to write the solution of the structural analysis\n",fileName);
        exit(EXIT_FAILURE);
    }

    fprintf(fid,"Time = %f\n",bm->t);
    for (i=0;i<bm->nkp;i++)
    {
        for(j=0;j<NDIM + NDOF_ND;j++)
        {
            fprintf(fid,"%1.8e ",bm->sol_pt[k][i][j]);
        }
        fprintf(fid,"\n");
    }

    for (i=0;i<bm->nelem;i++)
    {
        for(j=0;j<NDIM + bm->ndof_el;j++)
        {
            fprintf(fid,"%1.8e ",bm->sol_mb[k][i][j]);
        }
        fprintf(fid,"\n");
    }
    fclose(fid);
    //}
}



void gebt_setLoads(Gebt *bm, double *loads, int load_no)
{
    int i;
    for(i=0;i<bm->ncond_mb;i++)
    {
        bm->mb_condition[i].dof[load_no] = i+1;
        bm->mb_condition[i].value[load_no] = loads[i];
        bm->mb_condition[i].value_current[load_no] = loads[i];
    }

    for(i=0;i<bm->ndistrfun;i++)
    {
        bm->distr_fun[i][0] = 1.0; //Constant chebyshev polynomial
    }
}

void gebt_setEndLoad(Gebt *bm, double loads, int load_no)
{
    int i;
    for(i=0;i<bm->ncond_mb;i++)
    {
        bm->pt_condition[1].dof[load_no+NSTRN] = i+1;
        bm->pt_condition[1].value[load_no] = loads;
        bm->pt_condition[1].value_current[load_no] = loads;
    }

    for(i=0;i<bm->ndistrfun;i++)
    {
        bm->distr_fun[i][0] = 1.0; //Constant chebyshev polynomial
    }
}


void gebt_analysis(Gebt *bm)
{

    int i;
    char *error="";


    // Transform every pointer >1D into 1D pointers
    RAVEL2(bm,coord,bm->nkp,NDIM,double);
    RAVEL2(bm,member,bm->nmemb,MEMB_CONST,int);
    RAVEL3(bm,material,bm->nmate,bm->ndof_el-NSTRN,NSTRN,double);
    RAVEL3(bm,sol_pt,bm->nstep,bm->nkp,NDIM+NDOF_ND,double);
    RAVEL3(bm,sol_mb,bm->nstep,bm->nelem,NDIM+bm->ndof_el,double);
    RAVEL3(bm,frame,bm->nframe,NDIM,NDIM,double);
    RAVEL2(bm,distr_fun,bm->ndistrfun,NSTRN,double);
    RAVEL2(bm,curvature,bm->ncurv,NDIM,double);
    RAVEL2(bm,init_cond,bm->nelem,12,double);
    
    //RAVEL2(bm,eigen_val,2,bm->nev+1,double);
    //RAVEL3(bm,eigen_vec_pt,bm->nev+1,bm->nkp,NDIM+NDOF_ND,double);
    //RAVEL3(bm,eigen_vec_mb,bm->nev+1,bm->nelem,NDIM+bm->ndof_el,double);

    // Empty pointers are given for the eigen values as we don't want to perform such analysis
    double *rav_eigen_val   = malloc(sizeof(double));
    double *rav_eigen_vec_pt= malloc(sizeof(double));
    double *rav_eigen_vec_mb= malloc(sizeof(double));

    // The material tensor is modified by fortran_analysis. Thus, we will give it a copy
    double *material_tmp = (double*) malloc(sizeof(double) * bm->nmate * (bm->ndof_el - NSTRN) * NSTRN);
    for(i=0;i<bm->nmate * (bm->ndof_el - NSTRN) * NSTRN;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
    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,
               &rav_sol_pt,&rav_sol_mb,error,bm->ncond_mb,bm->ntimefun,&rav_frame,&(bm->mb_condition),&rav_distr_fun,
               &rav_curvature,bm->omega_a0,bm->omega_a_tf,bm->v_root_a0,bm->v_root_a_tf,bm->simu_time,
               &bm->time_function,bm->analysis_flag,&rav_init_cond,&bm->nev,&rav_eigen_val,&rav_eigen_vec_pt,&rav_eigen_vec_mb,bm->initialize);

    // Transform 1D pointers to 2D pointers if necessary
    UNRAVEL2(bm,coord,bm->nkp,NDIM,double);
    UNRAVEL2(bm,member,bm->nmemb,MEMB_CONST,int);
    UNRAVEL3(bm,material,bm->nmate,bm->ndof_el-NSTRN,NSTRN,int);
    UNRAVEL3(bm,sol_pt,bm->nstep,bm->nkp,NDIM+NDOF_ND,double);
    UNRAVEL3(bm,sol_mb,bm->nstep,bm->nelem,NDIM+bm->ndof_el,double);
    UNRAVEL3(bm,frame,bm->nframe,NDIM,NDIM,double);
    UNRAVEL2(bm,distr_fun,bm->ndistrfun,NSTRN,double);
    UNRAVEL2(bm,curvature,bm->ncurv,NDIM,double);
    UNRAVEL2(bm,init_cond,bm->nelem,12,double);
    //UNRAVEL2(eigen_val,2,nev+1,double);
    //UNRAVEL3(eigen_vec_pt,nev+1,nkp,NDIM+NDOF_ND,double);
    //UNRAVEL3(eigen_vec_mb,nev+1,nelem,NDIM+ndof_el,double);

    // Set initialize to zero to perform dynamic analysis without redoing the initialization step
    bm->initialize = 0;
    bm->lastStep = bm->nstep-1;
    bm->t += bm->nstep * bm->dt;

    free(rav_eigen_val);
    free(rav_eigen_vec_mb);
    free(rav_eigen_vec_pt);
    
    free(material_tmp);



}

/*
This updates the initial condition, but when the analysis does its initial step, it does not work
void gebt_updateInitialConditions(Gebt *bm)
{
    int i,j;
    for (i=0;i<bm->nelem;i++)
    {
        // copy u_dot, theta_dot from the initial conditions
        for (j=0;j<NSTRN;j++)
        {
            bm->init_cond[i][j+NSTRN] = bm->init_cond[i][j];
        }
        // copy u, theta from the solution member
        for (j=0;j<NSTRN;j++)
        {
            bm->init_cond[i][j] = bm->sol_mb[bm->nstep-1][i][j+NDIM];
        }
        
    }
}
*/

#define WRITEF3(ptr,fid,type,n1,n2,n3)  for (i_=0; i_<(n1); i_++){\
                                     for (j_=0; j_<(n2); j_++){\
                                     for (k_=0; k_<(n3); k_++){\
                                      fprintf(fid,#type " ",ptr[i_][j_][k_]);}\
                                      fprintf(fid,"\n");}\
                                      fprintf(fid,"\n");}

#define WRITEF2(ptr,fid,type,n1,n2)  for (i_=0; i_<(n1); i_++){\
                                     for (j_=0; j_<(n2); j_++){\
                                      fprintf(fid,#type " ",ptr[i_][j_]);}\
                                      fprintf(fid,"\n");}

#define WRITEF1(ptr,fid,type,n1)  for (i_=0; i_<(n1); i_++){\
                                      fprintf(fid,#type " ",ptr[i_][j_]);\
                                      fprintf(fid,"\n");}

void gebt_writeGebtToFile(Gebt *bm,char *fname)
{
    int i,j,k;
    FILE *fid = fopen(fname,"w");

    fprintf(fid,"nkp       = %d\n",bm->nkp);
    fprintf(fid,"nelem     = %d\n",bm->nelem);
    fprintf(fid,"ndof      = %d\n",bm->ndof_el);
    fprintf(fid,"nmemb     = %d\n",bm->nmemb);
    fprintf(fid,"ncond_pt  = %d\n",bm->ncond_pt);
    fprintf(fid,"nmate     = %d\n",bm->nmate);
    fprintf(fid,"nframe    = %d\n",bm->nframe);
    fprintf(fid,"ndistrfun = %d\n",bm->ndistrfun);
    fprintf(fid,"ncurv     = %d\n",bm->ncurv);
    fprintf(fid,"Coordinates: \n");
    WRITEF2(bm->coord, fid, %f, bm->nkp, NDIM);
    fprintf(fid,"Member: \n");
    WRITEF2(bm->member, fid, %d, bm->nmemb, MEMB_CONST);
    fprintf(fid,"Material: \n");
    WRITEF3(bm->material, fid, % 1.6E, bm->nmate, bm->ndof_el - NSTRN, NSTRN);
    fprintf(fid,"Frame: \n");
    WRITEF3(bm->frame, fid, %f, bm->nframe, NDIM, NDIM);
    fprintf(fid,"Distributed Functions: \n");
    WRITEF2(bm->distr_fun,fid,  %f, bm->ndistrfun, NSTRN);
    fprintf(fid,"Curvatures: \n");
    WRITEF2(bm->curvature,fid,  %f, bm->ncurv, NDIM);

    fclose(fid);
}


void freeGebt(Gebt *bm)
{
    
    printf("Freeing the memory...\n");
    int i;
    DEALLOCATE3(bm->sol_mb,bm->nstep,bm->nelem);
    DEALLOCATE3(bm->sol_pt,bm->nstep,bm->nkp);
    DEALLOCATE2(bm->init_cond,bm->nelem);
    DEALLOCATE2(bm->curvature,bm->ncurv);
    DEALLOCATE2(bm->distr_fun,bm->ndistrfun);
    DEALLOCATE3(bm->frame,bm->nframe,NDIM);
    DEALLOCATE3(bm->material,bm->nmate,bm->ndof_el-NSTRN);
    DEALLOCATE2(bm->member,bm->nmemb);
    DEALLOCATE2(bm->coord,bm->nkp);

    // Free time-functions
    //for (i=0;i<bm->ntimefun;i++)
    //    freeTimeFunction(&time_function[i]);
    DEALLOCATE1(bm->time_function); // freeTF

    DEALLOCATE1(bm->mb_condition); // freePrescriInf
    DEALLOCATE1(bm->pt_condition); //freePrescriInf
    printf("Done!\n");
}





// Library utils
static void gebt_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 gebt_setMaterial(Gebt *bm)
{
    int i,j,k;
    int ncomp, nmass;
    double ***matComp, ***matMass;
    double *xcomp, *xmass;
    double x;

    char*   line = NULL; //don't forget to free! 
    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
        gebt_interpolateMatrix(bm->material[i],x,xcomp,matComp,ncomp,NSTRN,NSTRN,0    ,0);
        if (bm->analysis_flag!=0) 
        {
            gebt_interpolateMatrix(bm->material[i],x,xmass,matMass,nmass,NSTRN,NSTRN,NSTRN,0);
        }

    }

    DEALLOCATE1(xcomp);
    DEALLOCATE3(matComp,ncomp,NSTRN);
    if (bm->analysis_flag!=0) 
    {
        DEALLOCATE1(xmass);
        DEALLOCATE3(matMass,ncomp,NSTRN);
    }

    free(line);

    fclose(COMP);
    fclose(MASS);
}