✨ new test-cases

.gitignore

@@ -16,6 +16,8 @@ run/**/main.o
 run/**/lib
 run/**/*.deb
 
+run/SweptBeam 
+
 src/fortran90/*.mod
 
 .vscode

run/CantileverEulerBeam/Compliance.dat

+1
+0.0
 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

run/CantileverEulerBeam/Gebt.ini

 [Beam]
 analysis_flag = 0
-nkp = 11
+nkp = 31
 niter = 1000
 nstep = 1
 ncond_pt = 2

run/CantileverEulerBeam/main.c

@@ -40,7 +40,7 @@ int main(int argc, char* argv[])
     
 
     printf("Outputting the results...\n");
-    gebt_writeSolToFile(bm);
+    gebt_writeSolToFile(bm,"output");
     printf("... Done!\n");
     
     return EXIT_SUCCESS;

run/CantileverEulerBeam/readData.py

 # Launch the Euler-Bernoulli comparison for the beam
 
 import numpy as np
+import matplotlib
 import matplotlib.pyplot as plt
 from readIni import IniReader
 import subprocess as sp
@@ -8,14 +9,14 @@ import subprocess as sp
 ini = IniReader("Gebt.ini");
 nkp = ini.getValue("Beam","nkp",typeCast=int);
 nmemb = nkp-1
-display = 0
+display = 1
 
 x = np.linspace(0,1);
 xcheck = np.linspace(0,1,nkp);
 
-for q in [1.0,10.0,100.0]:
+EI = 12.5;
+for q in [5 * EI*1*1]:
 
-    EI = 12.5;
     w = q * x*x *(6-4*x+x*x)/24/EI;
     wcheck = q * xcheck*xcheck *(6-4*xcheck+xcheck*xcheck)/24/EI;
 
@@ -29,16 +30,23 @@ for q in [1.0,10.0,100.0]:
     err_inf = max(np.fabs(pts[:,5]+pts[:,2]-wcheck));
     print("Maximal error: %e"%(err_inf));
 
+    matplotlib.rcParams.update({'font.size': 16})
+
     if(display):
         #plt.plot(mem[:,0]+mem[:,3],mem[:,4]+mem[:,1]);
         plt.figure()
-        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.plot(mem[:,0]+mem[:,3],mem[:,5]+mem[:,2],'.r',label='_nolegend_');
+        plt.plot(x,w,'k-',label="E-B",lw=3);
+        plt.plot(pts[:,0]+pts[:,3],pts[:,5]+pts[:,2],'.-r',label="GEBT",lw=3,ms=12);
 
         plt.xlabel('x');
         plt.ylabel('w(x)');
-        plt.title('q = %1.1f [N/m]'%(q));
+        #plt.title('q = %1.1f [N/m]'%(q));
+        plt.grid("True");
+        plt.xlim([0,1]);
+        plt.ylim([0,w[-1]]);
+        plt.tight_layout;
+        plt.legend();
     
     #plt.clf()
 if(display):

run/CantileverMoment/readData.py

@@ -2,9 +2,12 @@
 
 import numpy as np
 import matplotlib.pyplot as plt
+import matplotlib
 from readIni import IniReader
 import subprocess as sp
 
+matplotlib.rcParams.update({'font.size': 16})
+
 ini = IniReader("Gebt.ini");
 nkp = ini.getValue("Beam","nkp",typeCast=int);
 L   = ini.getValue("Beam","L"  ,typeCast=float);
@@ -19,27 +22,31 @@ xcheck = np.linspace(0,1,nkp);
 #print(np.linalg.inv(a));
 
 EI = 86.9e3
-lam = 1.95;
-M = lam * np.pi * EI / L
-print(M);
+for lam in [0.5,1.0,1.5]:
+    M = lam * np.pi * EI / L
+    print(M);
 
-sp.run(["./gebt", "-%e"%(M)]);
+    sp.run(["./gebt", "-%e"%(M)]);
 
 
-data_pts = np.loadtxt("output000.dat",max_rows=nkp);
-data_mem = np.loadtxt("output000.dat",skiprows=nkp);
+    data_pts = np.loadtxt("output000.dat",max_rows=nkp);
+    data_mem = np.loadtxt("output000.dat",skiprows=nkp);
 
-pts = data_pts;
-mem = data_mem;
+    pts = data_pts;
+    mem = data_mem;
 
-if(display):
-    #plt.plot(mem[:,0]+mem[:,3],mem[:,4]+mem[:,1]);
-    plt.plot(pts[:,0]+pts[:,3],pts[:,5]+pts[:,2],'.-k');
-    #plt.plot(mem[:,0]+mem[:,3],mem[:,5]+mem[:,2],'.r');
+    if(display):
+        #plt.plot(mem[:,0]+mem[:,3],mem[:,4]+mem[:,1]);
+        line = plt.plot(pts[:,0]+pts[:,3],pts[:,5]+pts[:,2],'.-',label=r"$\lambda$ = %1.1f"%(lam),lw=2);
+        #plt.plot(mem[:,0]+mem[:,3],mem[:,5]+mem[:,2],'.r');
 
     plt.xlabel('x');
     plt.ylabel('w(x)');
     plt.axis('scaled');
+    plt.legend(loc="lower right")
+    plt.xlim([-4,10]);
+    plt.ylim([-2,8]);
+    plt.grid("True");
 
 
 print("u_1 = %f"%(L-(pts[-1,0]+pts[-1,3])));

run/NREL5MW_convergence/Gebt.ini

+[Beam]
+analysis_flag = 0
+nkp = 2
+niter = 1000
+nstep = 1
+ncond_pt = 2
+initFromFile = 0
+L=61.5
+dt = 0.002
\ No newline at end of file

run/NREL5MW_convergence/Makefile

+# NOTE ON THIS MAKEFILE
+# As the library is not installed in a default repo, i.e. /usr/lib or /usr/local/lib, you need to give it the path to the library at runtime. This is the rpath argument.
+DIR=$(shell dirname `pwd`)
+
+default:
+	cp -r ../../lib ./
+	gcc -g -fPIC -c -o main.o main.c -I../../include
+	gcc -g -o gebt main.o -lgebt -L./lib
+
+clean:
+	rm -f main.o
+	rm -f gebt

run/NREL5MW_convergence/Twist.dat

+49
+0.0000000E+00 1.3308000E+01  
+1.9987500E-01 1.3308000E+01  
+1.1998650E+00 1.3308000E+01  
+2.1998550E+00 1.3308000E+01  
+3.1998450E+00 1.3308000E+01  
+4.1998350E+00 1.3308000E+01  
+5.1998250E+00 1.3308000E+01  
+6.1998150E+00 1.3308000E+01  
+7.1998050E+00 1.3308000E+01  
+8.2010250E+00 1.3308000E+01  
+9.1997850E+00 1.3308000E+01  
+1.0199775E+01 1.3308000E+01  
+1.1199765E+01 1.3181000E+01  
+1.2199755E+01 1.2848000E+01  
+1.3200975E+01 1.2192000E+01  
+1.4199735E+01 1.1561000E+01  
+1.5199725E+01 1.1072000E+01  
+1.6199715E+01 1.0792000E+01  
+1.8200925E+01 1.0232000E+01  
+2.0200290E+01 9.6720000E+00  
+2.2200270E+01 9.1100000E+00  
+2.4200250E+01 8.5340000E+00  
+2.6200230E+01 7.9320000E+00  
+2.8200825E+01 7.3210000E+00  
+3.0200190E+01 6.7110000E+00  
+3.2200170E+01 6.1220000E+00  
+3.4200150E+01 5.5460000E+00  
+3.6200130E+01 4.9710000E+00  
+3.8200725E+01 4.4010000E+00  
+4.0200090E+01 3.8340000E+00  
+4.2200070E+01 3.3320000E+00  
+4.4200050E+01 2.8900000E+00  
+4.6200030E+01 2.5030000E+00  
+4.8201240E+01 2.1160000E+00  
+5.0199990E+01 1.7300000E+00  
+5.2199970E+01 1.3420000E+00  
+5.4199950E+01 9.5400000E-01  
+5.5199940E+01 7.6000000E-01  
+5.6199930E+01 5.7400000E-01  
+5.7199920E+01 4.0400000E-01  
+5.7699915E+01 3.1900000E-01  
+5.8201140E+01 2.5300000E-01  
+5.8699905E+01 2.1600000E-01  
+5.9199900E+01 1.7800000E-01  
+5.9699895E+01 1.4000000E-01  
+6.0199890E+01 1.0100000E-01  
+6.0699885E+01 6.2000000E-02  
+6.1199880E+01 2.3000000E-02  
+6.1500000E+01 0.0000000E+00
\ No newline at end of file

run/NREL5MW_convergence/main.c

+#include <stdlib.h>
+#include <stdio.h>
+#include <math.h>
+
+#include "gebt.h"
+
+
+int main(int argc, char* argv[])
+{
+    int i;
+    double load;
+
+    // Check the input argument
+    if(argc<2)
+    {
+        printf("Please specify the uniform load value as arguments\n");
+        printf("Example: ./gebt 100\n");
+        exit(EXIT_FAILURE);
+    }
+
+    // Initializing the Beam
+    printf("Initializing the beam...\n");
+    Gebt *bm = malloc(sizeof(Gebt));
+    initGebt(bm);
+    printf("...Done!\n");
+
+    gebt_setTwistFromFile(bm,"Twist.dat");
+
+    // Creating an array with the uniform load for each member
+    printf("Creating the loads...\n");
+    double *loads = malloc(sizeof(double)*bm->nmemb);
+    for(i=0;i<bm->nmemb;i++)
+    {
+        loads[i] = atof(argv[1]);
+    }
+    gebt_setLoads(bm,loads,1);
+    printf("...Done!\n");
+
+    char fname[32] = "output";
+    
+    // Performing the analysis
+    printf("Performing the analysis...\n");
+    
+    gebt_analysis(bm);
+    gebt_writeSolToFile(bm,fname);
+    printf("... Ok! \n");
+    
+    freeGebt(bm);
+    free(bm);
+    
+    return EXIT_SUCCESS;
+}
\ No newline at end of file

run/NREL5MW_convergence/prepro.py

+# This script will arange the mass and stiffness matrix from the fast input file
+# into correctly formatted input files
+
+# Also, the stiffness matrix is inverted by definition of the compliance matrix
+
+
+import numpy as np
+
+
+fid = open("ComplianceAndMass.dat","r");
+
+lines = fid.readlines();
+
+nentry = 49;
+
+stiffness = np.zeros((nentry,6,6))
+mass      = np.zeros((nentry,6,6))
+x         = np.zeros(nentry);
+
+k=0
+for i in np.arange(5,725+1,15):
+    x[k] = float(lines[i])
+    k=k+1;
+
+k=0
+for i in np.arange(6,726+1,15):
+    for j in range(6):
+        stiffness[k,j,:] = np.array(lines[i+j].split()).astype(np.float);
+    k=k+1;
+
+
+k=0
+for i in np.arange(13,734+1,15):
+    for j in range(6):
+        mass[k,j,:] = np.array(lines[i+j].split()).astype(np.float);
+    k=k+1;
+
+
+with open("Compliance.dat","w") as fid:
+    fid.write("%d\n"%(nentry));
+    for i in range(nentry):
+        comp = np.linalg.inv(stiffness[i])
+        fid.write("%1.6f\n"%(x[i]));
+        for j in range(6):
+            for k in range(6):
+                fid.write("% 1.6E "%(comp[j,k]));
+            fid.write("\n");
+        fid.write("\n");
+
+with open("Mass.dat","w") as fid:
+    fid.write("%d\n"%(nentry));
+    for i in range(nentry):
+        fid.write("%1.6f\n"%(x[i]));
+        for j in range(6):
+            for k in range(6):
+                fid.write("% 1.6E "%(mass[i,j,k]));
+            fid.write("\n");
+        fid.write("\n");
\ No newline at end of file

run/NREL5MW_convergence/readData.py

+# Launch the Euler-Bernoulli comparison for the beam
+
+import numpy as np
+import matplotlib.pyplot as plt
+from readIni import IniReader
+import subprocess as sp
+
+ini = IniReader("Gebt.ini");
+nkp = ini.getValue("Beam","nkp",typeCast=int);
+L   = ini.getValue("Beam","L"  ,typeCast=float);
+nmemb = nkp-1
+display = 1
+
+F = 1e4;
+
+u2 = [];
+u3 = [];
+nkps = [];
+
+nkp=2;
+for i in range(7):
+
+    nkp = nkp*2;
+    print("Simulation with nkp = %d"%(nkp));
+
+    with open('Gebt.ini', 'r') as file :
+        filedata = file.read()
+        filedata = filedata.replace('nkp = %d'%(nkp/2), 'nkp = %d'%(nkp))
+
+    with open('Gebt.ini', 'w') as file:
+        file.write(filedata)
+
+    out = sp.run(["./gebt", "-%e"%(F)],check=True);
+
+    data_pts = np.loadtxt("output%03d.dat"%(0),max_rows=nkp);
+    data_mem = np.loadtxt("output%03d.dat"%(0),skiprows=nkp);
+
+    pts = data_pts;
+    mem = data_mem;
+
+    u2.append(pts[-1,4])
+    u3.append(pts[-1,5])
+    nkps.append(nkp);
+
+print("Done!");
+
+ref = -9.75404289;
+with open('Gebt.ini', 'r') as file :
+        filedata = file.read()
+        filedata = filedata.replace('nkp = %d'%(nkp), 'nkp = %d'%(2))
+with open('Gebt.ini', 'w') as file:
+    file.write(filedata)
+
+plt.plot(nkps,u2,'.-');
+plt.plot(nkps,np.ones(len(nkps))*(ref),'k--');
+ax = plt.gca();
+ax.set_xscale('log', base=2);
+plt.grid(True)
+#plt.plot(nkps,u3,'.-');
+
+u2 = np.array(u2);
+nkps = np.array(nkps);
+error = np.fabs(u2-(ref));
+plt.figure();
+plt.loglog(nkps,error,'.-k');
+plt.loglog(nkps,nkps.astype(float)**(-2),'k--')
+ax = plt.gca();
+ax.set_xscale('log', base=2);
+plt.grid(True);
+
+print(u2);
+
+plt.show();
\ No newline at end of file

run/NREL5MW_dynamic/Gebt.ini

 [Beam]
 analysis_flag = 2
-nkp = 25
+nkp = 64
 niter = 1
-nstep = 10
+nstep = 5
 ncond_pt = 2
 initFromFile = 0
 L=61.5

run/NREL5MW_dynamic/main.c

@@ -29,7 +29,7 @@ int main(int argc, char* argv[])
     gebt_setLoads(bm,loads,1);   // 1 -> Fy (flapwise)
     for(i=0;i<bm->nmemb;i++)
     {
-        loads[i] = 10e3;
+        loads[i] = 1e4;
     }
     gebt_setLoads(bm,loads,2);   // 2 -> Fz (edgewise)
     printf("...Done!\n");

run/NREL5MW_dynamic/readData.py

@@ -24,7 +24,7 @@ for f in ["output"]:
     u3[f] = []
     phi1[f] = []
     for k in range(100):
-        for j in range(10):
+        for j in range(5):
             print("%s_%d_%03d.dat"%(f,k,j))
             data_pts = np.loadtxt("%s_%d_%03d.dat"%(f,k,j),max_rows=nkp);
             data_mem = np.loadtxt("%s_%d_%03d.dat"%(f,k,j),skiprows=nkp);
@@ -98,13 +98,23 @@ for f in []:
     u3[f] = np.array(u3[f]);
     t = np.linspace(0,u3[f].size*dt,u3[f].size);
 
+
+vtip = np.fromfile("/Users/ftrigaux/Documents/Beams/CBeams2/vtip.bin");
+wtip = np.fromfile("/Users/ftrigaux/Documents/Beams/CBeams2/wtip.bin");
+tv   = np.linspace(0,5,500);
 plt.figure()
 for f in ["output"]:
-    plt.plot(t,u2[f]);
-    plt.plot(t,u3[f]);
-    plt.plot(t,phi1[f]);
+    plt.plot(t,u2[f],'-k',label='Flap Gebt');
+    plt.plot(tv,vtip,'--k',label='Flap E-B');
+    plt.plot(t,u3[f],'-r',label='Edge Gebt');
+    plt.plot(tv,wtip,'--r',label='Edge E-B');
+    #plt.plot(t,phi1[f])
     plt.grid(True);
-    plt.title("GEBT");
+    plt.xlim([0,4]);
+    plt.ylim([-2,12]);
+    plt.legend();
+    plt.xlabel('Time [s]');
+    plt.ylabel('Deflection [m]');
 
 print("Done!");