volume_integral_2d_N_8_100_elements.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]

print(nodes)

[[-1. 1. ] [-1. 0.9] [-1. 0.8] [-0.9 0.8] [-0.8 0.8] [-0.8 0.9] [-0.8 1. ] [-0.9 1. ] [-0.9 0.9]]

# 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 / 0.6^2)
c_x = 1.
F_x = c_x * u
F_xi = 10. * F_x

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 / 0.6^2)
    c_x = 1.
    F_x = c_x * u
    F_xi = 10. * F_x
    print(element_tag)
    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],
                                                                                   -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 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99

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