🍺 Added function for distributed loads

.gitignore

 obj/*
 bin/*
+run/gebt*
+run/output*

run/Beam.ini

 [Beam]
 analysis_flag = 0
-nkp = 2
+nkp = 400
 niter = 1000
 nstep = 1
-nmemb = 1
 ncond_pt = 2
-nmate = 1
 nframe = 0
-ndistrfun = 1
 ncurv = 0
-ncond_mb = 1
\ No newline at end of file
+initFromFile = 0
+L=1.0
\ No newline at end of file

run/Compliance.dat

+0.00000000000000E+00     0.00000000000000E+00     0.00000000000000E+00     0.00000000000000E+00     0.00000000000000E+00    0.00000000000000E+00
+0.00000000000000E+00     0.00000000000000E+00     0.00000000000000E+00     0.00000000000000E+00     0.00000000000000E+00    0.00000000000000E+00
+0.00000000000000E+00     0.00000000000000E+00     0.00000000000000E+00     0.00000000000000E+00     0.00000000000000E+00    0.00000000000000E+00
+0.00000000000000E+00     0.00000000000000E+00     0.00000000000000E+00     0.00000000000000E+00     0.00000000000000E+00    0.00000000000000E+00
+0.00000000000000E+00     0.00000000000000E+00     0.00000000000000E+00     0.00000000000000E+00     1.33333333333333E-02    0.00000000000000E+00
+0.00000000000000E+00     0.00000000000000E+00     0.00000000000000E+00     0.00000000000000E+00     0.00000000000000E+00    1.33333333333333E-02
\ No newline at end of file

run/readData.py

+
+
 import numpy as np
 import matplotlib.pyplot as plt
-npts = 3
-nmemb = 10
+from readIni import IniReader
+
+ini = IniReader("Beam.ini");
+nkp = ini.getValue("Beam","nkp",typeCast=int);
+nmemb = nkp-1
 
 
+x = np.linspace(0,1);
+q  = 1;
+w_obj = 0.01;
+EI = 6 * q * 1 *(6-4+1)/24 / w_obj;
+print(EI);
+print(1/EI)
+w = 6 * q * x*x *(6-4*x+x*x)/24/EI;
+
 for k in range(1):
     data = np.loadtxt("output%03d.dat"%(k));
 
-    pts = data[0:npts,:];
-    mem = data[npts:,:];
+    pts = data[0:nkp,:];
+    mem = data[nkp:,:];
 
     #plt.plot(mem[:,0]+mem[:,3],mem[:,4]+mem[:,1]);
-    plt.plot(mem[:,0]+mem[:,3],mem[:,5]+mem[:,2]);
+    plt.plot(pts[:,0]+pts[:,3],pts[:,5]+pts[:,2],'.-b');
+    plt.plot(mem[:,0]+mem[:,3],mem[:,5]+mem[:,2],'.r');
+
+    plt.plot(x,w,'k');
     #plt.axis('scaled');
     #plt.xlim([-1,1]);
     #plt.ylim([-1,1]);

src/fortran90/Analysis.f90

@@ -94,6 +94,7 @@ IF(DEBUG) THEN
   IF(FileOpen(IOUT, deb_name, 'UNKNOWN', 'WRITE',error))	 RETURN 
 ENDIF  
 WRITE(*,*) error
+
 CALL Preprocess(nkp,nelem,ndof_el,member,material,frame,coord,curvature,&
                     & dof_con,memb_info,error)
 IF(error/='') RETURN
@@ -182,7 +183,6 @@ IF(analysis_flag==2) THEN
 
 ENDIF
 
-
 DO i=1,nstep
     ! evaluate the prescribed information for this step
 	!---------------------------------------------------
@@ -210,7 +210,6 @@ DO i=1,nstep
 		IF(DEBUG) WRITE(IOUT,*) "STEP=",i
 	ENDIF
 
-
     IF(niter==1) THEN  ! solve the linear system for ith step
       CALL LinearSolution(ndof_el,memb_info,v_root_a,omega_a,member,error,&
 			& ncond_mb,mb_condition,distr_fun,(dof_con),x,init_cond)

src/fortran90/Wrappers.f90

@@ -83,8 +83,6 @@ SUBROUTINE  c_Analysis(nkp,nelem,ndof_el,nmemb,ncond_pt,nmate, nframe,ndistrfun,
         TYPE(PrescriInf),POINTER       :: f_pt_condition(:)
         TYPE(PrescriInf),POINTER       :: f_mb_condition(:)
 
-        WRITE(*,*) simu_time
-
         
         WRITE(*,*) "Converting all pointers..."
         CALL c_f_pointer(coord,f_coord, [nkp,NDIM]);
@@ -102,7 +100,6 @@ SUBROUTINE  c_Analysis(nkp,nelem,ndof_el,nmemb,ncond_pt,nmate, nframe,ndistrfun,
         !CALL c_f_pointer(omega_a_tf, f_omega_a_tf, [NDIM]);
         !CALL c_f_pointer(v_root_a_tf, f_v_root_a_tf, [NDIM]);
         !CALL c_f_pointer(simu_time, f_simu_time, [2])
-        WRITE(*,*) ntimefun
         CALL c_f_pointer(time_function, f_time_function, [ntimefun]);
         CALL c_f_pointer(init_cond, f_init_cond, [nelem,12]);
         CALL c_f_pointer(eigen_val, f_eigen_val, [2,nev+1]);
@@ -112,6 +109,7 @@ SUBROUTINE  c_Analysis(nkp,nelem,ndof_el,nmemb,ncond_pt,nmate, nframe,ndistrfun,
         ! --------------------------------------------------
         ! Call the original fortran function
         ! --------------------------------------------------
+        WRITE(*,*)   '--- Calling GEBT ---'
         CALL Analysis(nkp,nelem,ndof_el,nmemb,ncond_pt,nmate, nframe,ndistrfun,ncurv,f_coord, &
         &  f_member,f_pt_condition,f_material,niter,nstep,f_sol_pt,f_sol_mb,f_error,ncond_mb, &
         &  ntimefun,f_frame,f_mb_condition,f_distr_fun,f_curvature,omega_a0,omega_a_tf, &

src/interface/beam.c

 #include <stdlib.h>
 #include <stdio.h>
 
-#include "beam.h"
 #include "readINIFile.h"
+#include "beam.h"
+
 
 int check;
 int i_,j_,k_; // macro iterators
 
 void initBeam(Beam *bm)
+{
+    int i,j,k;
+    char *section;
+    dataINI fileBeam = readINIFile("Beam.ini");
+
+    char*   line = NULL;
+    ssize_t read;
+    int     frmt;
+    size_t  len = 0;
+
+    /* Set the mesh parameters */ 
+    /*section = "Beam";
+    SET_OR_DEFAULT(&(bm->analysis_flag),%d,fileBeam,section,"analysis_flag",2);
+    SET_OR_DEFAULT(&(bm->nkp),%d,fileBeam,section,"nkp",2);
+    SET_OR_DEFAULT(&(bm->nframe),%d,fileBeam,section,"nframe",0);// number of different cross-sectional frames
+    SET_OR_DEFAULT(&(bm->ncurv),%d,fileBeam,section,"ncurv",0);// number of sets of initial curvature/twists
+    SET_OR_DEFAULT(&(bm->L),%lf,fileBeam,section,"L",1.0);// number of sets of initial curvature/twists
+    SET_OR_DEFAULT(&(bm->initFromFile),%d,fileBeam,section,"initFromFile",0);// number of sets of initial curvature/twists
+
+    */
+    section = "Beam";
+    SET_OR_DEFAULT(&(bm->analysis_flag),%d,fileBeam,section,"analysis_flag",2);
+    SET_OR_DEFAULT(&(bm->nkp),%d,fileBeam,section,"nkp",2);
+    SET_OR_DEFAULT(&(bm->nframe),%d,fileBeam,section,"nframe",0);// number of different cross-sectional frames
+    SET_OR_DEFAULT(&(bm->ncurv),%d,fileBeam,section,"ncurv",0);// number of sets of initial curvature/twists
+
+    SET_OR_DEFAULT(&(bm->L),%lf,fileBeam,section,"L",1.0);// number of sets of initial curvature/twists
+    SET_OR_DEFAULT(&(bm->initFromFile),%d,fileBeam,section,"initFromFile",0);// number of sets of initial curvature/twists
+
+    SET_OR_DEFAULT(&(bm->niter),%d,fileBeam,section,"niter",100);
+    SET_OR_DEFAULT(&(bm->nstep),%d,fileBeam,section,"nstep",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->ncond_pt = 2; // condition for clamped-free beam
+
+    bm->ntimefun = 0;
+
+    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,fileBeam,section,"omegaX",0.0);
+        SET_OR_DEFAULT(&bm->omega_a0[1],%lf,fileBeam,section,"omegaY",0.0);
+        SET_OR_DEFAULT(&bm->omega_a0[2],%lf,fileBeam,section,"omegaZ",0.0);
+        SET_OR_DEFAULT(&bm->v_root_a0[0],%lf,fileBeam,section,"vrootX",0.0);
+        SET_OR_DEFAULT(&bm->v_root_a0[1],%lf,fileBeam,section,"vrootY",0.0);
+        SET_OR_DEFAULT(&bm->v_root_a0[2],%lf,fileBeam,section,"vrootZ",0.0);
+    }
+
+    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->nkp); 
+        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] = 0;   // frame number
+        bm->member[i][5] = 1;   // number of division
+        bm->member[i][6] = 0;   // initial curvature number
+    }
+
+    // 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[0].value[i];
+    }
+
+    // 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);}
+    
+    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++){
+        for(j=0;j<NDIM;j++){
+            for(k=0;k<NDIM;k++){
+                bm->frame[i][j][k]= (i==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++)
+        {
+            // TODO - Allow for curvature
+            printf("Cannot set curvature\n");
+            exit(EXIT_FAILURE);
+        }
+    }
+
+    // This must be initialize but has no consequence as no time-functions are used
+    bm->simu_time[0] = 0.0;
+    bm->simu_time[1] = 1.0;
+
+    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);
+    if(bm->analysis_flag==2){
+
+        if (bm->initFromFile==1)
+        {
+            // READ INITIAL CONIDITIONS FROM THE FILE
+            FILE *INIT = fopen("Beam_initCond.ini","r"); 
+            if (INIT==NULL)
+            {
+                printf("Cannot read %s\n","Beam_init.dat");
+                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);
+        }
+        else
+        {
+             for (i=0;i<bm->nelem;i++)
+             {
+                 for (j=0;j<2*NSTRN;j++)
+                 {
+                    bm->init_cond[i][6] = 0.0;
+                 }
+             }
+
+        }
+
+
+    }
+
+
+    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
+
+}
+
+
+void initBeamFromGebtFile(Beam *bm)
 {
   int i,j;
   char *section, dataFileName[32];
@@ -279,7 +562,7 @@ void beam_writeSolToFile(Beam *bm)
 
         for (i=0;i<bm->nkp;i++)
         {
-            for(j=0;j<15;j++)
+            for(j=0;j<NDIM + NDOF_ND;j++)
             {
                 fprintf(fid,"%1.8e ",bm->sol_pt[k][i][j]);
             }
@@ -288,7 +571,7 @@ void beam_writeSolToFile(Beam *bm)
 
         for (i=0;i<bm->nelem;i++)
         {
-            for(j=0;j<15;j++)
+            for(j=0;j<NDIM + bm->ndof_el;j++)
             {
                 fprintf(fid,"%1.8e ",bm->sol_mb[k][i][j]);
             }
@@ -299,6 +582,23 @@ void beam_writeSolToFile(Beam *bm)
 }
 
 
+
+void beam_setLoads(Beam *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 freeBeam(Beam *bm)
 {
     

src/interface/beam.h

@@ -108,10 +108,15 @@ typedef struct
     double               **init_cond;
     int                  analysis_flag;
     int                  nev;             // upon return, nev is the number of converged eigen values might be the original nev given by the user in the inputs or nev+1 due to complex conjugate pair 
+
+    // Personal add
+    double              L;                // Beam Length
+    int                 initFromFile;
 } Beam;
 
 
 
 void initBeam(Beam *bm);
+void beam_setLoads(Beam *bm, double *loads, int load_no);
 void beam_writeSolToFile(Beam *bm);
 void freeBeam(Beam *bm);
\ No newline at end of file

src/interface/main.c

@@ -7,10 +7,20 @@
 
 int main(int argc, char* argv[])
 {
+    int i;
+
     printf("Reading...\n");
     Beam *bm = malloc(sizeof(Beam));
     initBeam(bm);
-    int i,j,k;
+    printf("...Done!\n");
+
+    printf("Creating the loads...\n");
+    double *loads = malloc(sizeof(double)*bm->nmemb);
+    for(i=0;i<bm->nmemb;i++)
+    {
+        loads[i] = 1.0;
+    }
+    beam_setLoads(bm,loads,2);
     printf("...Done!\n");