surface_integral_2d_N_8_100_elements.sagews

%typeset_mode True

var('xi, eta, x, y')

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$, $\displaystyle x$, $\displaystyle y$)

nodes = [[ 0.8,  0.8],
         [ 0.8,  0.7],
         [ 0.8,  0.6],
         [ 0.9,  0.6],
         [ 1. ,  0.6],
         [ 1. ,  0.7],
         [ 1. ,  0.8],
         [ 0.9,  0.8],
         [ 0.9 , 0.7]]


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.6^2)
c_x = 1.
c_y = 1.

expand(u)

e^(-(1.71193772834421e-32)*eta^2*xi^4 - (8.55968864172105e-33)*eta^3*xi^2 - (2.56790659251631e-32)*eta*xi^4 - (2.13992216043026e-33)*eta^4 - (3.08395284618099e-17)*eta^2*xi^2 - (3.08395284618099e-17)*eta*xi^3 - (1.06996108021513e-32)*xi^4 - (1.54197642309050e-17)*eta^3 - (5.08852219619863e-16)*eta*xi^2 - (3.08395284618099e-17)*xi^3 - 0.0277777777777779*eta^2 - 0.0277777777777782*xi^2 - 0.388888888888889*eta - 0.500000000000000*xi - 3.61111111111111)

F_x = c_x * u
F_xi = 10. * F_x

F_y = c_y * u
F_eta = 10. * F_y

# Surface term

integral_00 = 0

integral_00 += L_xi[0](xi = 1) * integrate(L_eta[0] * 100. * F_xi(xi = 1), eta, -1, 1)
integral_00 += -L_xi[0](xi = -1) * integrate(L_eta[0] * 100. * F_xi(xi = -1), eta, 1, -1)
integral_00 +=  L_eta[0](eta = 1) * integrate(L_xi[0] * 100. * F_eta(eta = 1), xi, 1, -1)
integral_00 += -L_eta[0](eta = -1) *integrate(L_xi[0] * 100. * F_eta(eta =-1), xi, -1, 1)
integral_00.n()

$\displaystyle 0.0000910758972167969$

# Surface term

integral_pq = []

for p in np.arange(N_LGL):
    for q in np.arange(N_LGL):
        integral_pq.append(L_xi[p](xi = 1) * integrate(L_eta[q] * 100. * F_xi(xi = 1), eta, -1, 1))
        integral_pq[-1] += -L_xi[p](xi = -1) * integrate(L_eta[q] * 100. * F_xi(xi = -1), eta, 1, -1)

for p in np.arange(N_LGL):
    for q in np.arange(N_LGL):
        index = p * N_LGL + q
#         print(index)
        print(integral_pq[indexx])

6.03618983589522 35.6118199995744 57.6542772170424 69.7110283535537 69.7110283535533 57.6542772170429 35.6118199995742 6.03618983589531 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 2.22059014363387 13.1008564405414 21.2098232837490 25.6452541541919 25.6452541541918 21.2098232837492 13.1008564405413 2.22059014363390

︠265b8f-40-2b32568-b0b5-4cf8-ab18-df3