volume_integral_2d_N_8_100_elements_cy_1.sagews

%typeset_mode True

import numpy as np

var('xi, eta')
N_LGL = 8

xi_LGL = [-1.0,           \
          -0.87174014851, \
          -0.591700181433,\
          -0.209299217902,\
          0.209299217902, \
          0.591700181433, \
          0.87174014851,  \
          1.0]

eta_LGL = [-1.0,           \
           -0.87174014851, \
           -0.591700181433,\
           -0.209299217902,\
           0.209299217902, \
           0.591700181433, \
           0.87174014851,  \
           1.0]

($\displaystyle \xi$, $\displaystyle \eta$)

element_nodes = np.load('files/element_nodes.npy')

nodes = element_nodes[0]

N_0 = (-1.0 / 4.0) * (1 - xi)  * (1 + eta) * (1 + xi - eta)
N_1 = (1.0 / 2.0)  * (1 - xi)  * (1 - eta**2)
N_2 = (-1.0 / 4.0) * (1 - xi)  * (1 - eta) * (1 + xi + eta)
N_3 = (1.0 / 2.0)  * (1 - eta) * (1 - xi**2)
N_4 = (-1.0 / 4.0) * (1 + xi)  * (1 - eta) * (1 - xi + eta)
N_5 = (1.0 / 2.0)  * (1 + xi)  * (1 - eta**2)
N_6 = (-1.0 / 4.0) * (1 + xi)  * (1 + eta) * (1 - xi - eta)
N_7 = (1.0 / 2.0)  * (1 + eta) * (1 - xi**2)

# x =   N_0 * nodes[0][0] \
#     + N_1 * nodes[1][0] \
#     + N_2 * nodes[2][0] \
#     + N_3 * nodes[3][0] \
#     + N_4 * nodes[4][0] \
#     + N_5 * nodes[5][0] \
#     + N_6 * nodes[6][0] \
#     + N_7 * nodes[7][0]

# y =   N_0 * nodes[0][1] \
#     + N_1 * nodes[1][1] \
#     + N_2 * nodes[2][1] \
#     + N_3 * nodes[3][1] \
#     + N_4 * nodes[4][1] \
#     + N_5 * nodes[5][1] \
#     + N_6 * nodes[6][1] \
#     + N_7 * nodes[7][1]

# L_xi and L_eta polynomials given by this code are tested
import numpy as np

L_xi  = []
L_eta = []

for i in np.arange(N_LGL):
    L_xi.append(1.)
    L_eta.append(1.)
    for j in np.delete(np.arange(N_LGL), i):
        L_xi[-1]  *= ((xi - xi_LGL[j]) / (xi_LGL[i] - xi_LGL[j]))
        L_eta[-1] *= ((eta - eta_LGL[j]) / (eta_LGL[i] - eta_LGL[j]))

u = e^(-(x^2 + y^2) / 0.4^2)
c_x = 1.
c_y = 1.
F_x = c_x * u
F_xi = 10. * F_x

F_y = c_y * u
F_eta = 10. * F_y

Error in lines 1-1 Traceback (most recent call last): File "/cocalc/lib/python2.7/site-packages/smc_sagews/sage_server.py", line 996, in execute exec compile(block+'\n', '', 'single') in namespace, locals File "", line 1, in <module> NameError: name 'y' is not defined

gauss_nodes = [-0.96816024, -0.83603111, -0.61337143, -0.32425342, 0., 0.32425342, 0.61337143, 0.83603111, 0.96816024]
gauss_weights = [0.08127438836156,
                 0.18064816069485,
                 0.26061069640294,
                 0.31234707704000,
                 0.33023935500126,
                 0.31234707704000,
                 0.26061069640294,
                 0.18064816069486,
                 0.08127438836157]
integral = 0
p = 0
q = 3
for weight, node  in zip(gauss_weights, gauss_nodes):
    integral += weight * numerical_integral(100. * diff(L_xi[p], xi) * F_xi * L_eta[q], -1, 1, params = [node])[0]
print integral

-367.619734412963

# Volume integral for all the elements for all combinations of p and q

gauss_nodes = [-0.96816024, -0.83603111, -0.61337143, -0.32425342, 0., 0.32425342, 0.61337143, 0.83603111, 0.96816024]
gauss_weights = [0.08127438836156,
                 0.18064816069485,
                 0.26061069640294,
                 0.31234707704000,
                 0.33023935500126,
                 0.31234707704000,
                 0.26061069640294,
                 0.18064816069486,
                 0.08127438836157]

volume_integral = np.zeros([N_LGL**2, element_nodes.shape[0]])

for element_tag in np.arange(element_nodes.shape[0]):
#     element_tag = 0
    x =   N_0 * element_nodes[element_tag][0][0] \
        + N_1 * element_nodes[element_tag][1][0] \
        + N_2 * element_nodes[element_tag][2][0] \
        + N_3 * element_nodes[element_tag][3][0] \
        + N_4 * element_nodes[element_tag][4][0] \
        + N_5 * element_nodes[element_tag][5][0] \
        + N_6 * element_nodes[element_tag][6][0] \
        + N_7 * element_nodes[element_tag][7][0]

    y =   N_0 * element_nodes[element_tag][0][1] \
        + N_1 * element_nodes[element_tag][1][1] \
        + N_2 * element_nodes[element_tag][2][1] \
        + N_3 * element_nodes[element_tag][3][1] \
        + N_4 * element_nodes[element_tag][4][1] \
        + N_5 * element_nodes[element_tag][5][1] \
        + N_6 * element_nodes[element_tag][6][1] \
        + N_7 * element_nodes[element_tag][7][1]

    u = e^(-(x^2 + y^2) / 0.4^2)
    c_x = 1.
    c_y = 1.
    F_x = c_x * u
    F_xi = 10. * F_x

    F_y = c_y * u
    F_eta = 10. * F_y
    print(element_tag)
    p = 0
    q = 0
#     for p in np.arange(N_LGL):
#         for q in np.arange(N_LGL):
    index = p * N_LGL + q
    for weight, node  in zip(gauss_weights, gauss_nodes):
        volume_integral[index][element_tag] += weight * numerical_integral(100. * diff(L_xi[p], xi) * F_xi * L_eta[q] + 100. * diff(L_eta[q], eta) * F_eta * L_xi[p],
                                                                           -1, 1, params = [node])[0]
#             print(index)

0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25

np.save('results/volume_integral_N_8_elementtag_0.npy', volume_integral)

Error in lines 1-1 Traceback (most recent call last): File "/cocalc/lib/python2.7/site-packages/smc_sagews/sage_server.py", line 996, in execute exec compile(block+'\n', '', 'single') in namespace, locals File "", line 1, in <module> NameError: name 'np' is not defined