readData.py 1.05 KB
Newer Older
1
2
3
4
5
6
7
# 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

8
ini = IniReader("Gebt.ini");
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
nkp = ini.getValue("Beam","nkp",typeCast=int);
nmemb = nkp-1
display = 0

x = np.linspace(0,1);
xcheck = np.linspace(0,1,nkp);

for q in [1.0,10.0,100.0]:

    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;

    sp.run(["./gebt",str(q)]);

    data = np.loadtxt("output%03d.dat"%(0));

    pts = data[0:nkp,:];
    mem = data[nkp:,:];

    err_inf = max(np.fabs(pts[:,5]+pts[:,2]-wcheck));
    print("Maximal error: %e"%(err_inf));

    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.xlabel('x');
        plt.ylabel('w(x)');
        plt.title('q = %1.1f [N/m]'%(q));
    
    #plt.clf()
if(display):
    plt.show();