VEM_poisson_2d.py

import numpy as np


def poly_geom(verts):
    x, y = verts[:, 0], verts[:, 1]

    area_components = x * np.roll(y, -1) - np.roll(x, -1) * y
    area = 0.5 * np.abs(np.dot(x, np.roll(y, 1)) - np.dot(y, np.roll(x, 1)))
    centroid = ((verts + np.roll(verts, -1, axis=0)) * np.tile(area_components, (2, 1)).T).sum(0) / (6 * area)

    diameter = 0

    for i in range(verts.shape[0] - 1):
        for j in range(i + 1, verts.shape[0]):
            diameter = max(diameter, np.linalg.norm(verts[i, :] - verts[j, :]))

    return area, diameter, centroid


def vectort(verts):
    w = np.zeros(shape=(verts.shape[0], 2))
    n_ort = np.zeros(shape=(verts.shape[0], 2))
    v = np.concatenate((verts, verts[0, :].reshape(-1, 1).T), axis=0)

    for i in range(verts.shape[0]):
        w[i, :] = v[i + 1, :] - v[i, :]
        n = np.array([w[:, 1], -w[:, 0]])
        n_ort[i, :] = n.T[i, :] / np.linalg.norm(n.T[i, :])
    return n_ort


def K(verts):
    area, diameter, _ = poly_geom(verts)
    # K = np.zeros(shape=(verts.shape[0], verts.shape[0]))
    W_ = W(verts)
    H_ = H(verts)

    P = np.matmul(H_, W_.T)
    # Ic = np.array([[0, 0, 0], [0, 1/diameter**2, 0], [0, 0, 1/diameter**2]])
    Ic = np.array([[0, 0, 0], [0, 1, 0], [0, 0, 1]])
    I = np.identity(verts.shape[0])

    # K = area * W * Ic * W.T + (I - P).T * (I - P)
    K = area * np.matmul(np.matmul(W_, Ic), W_.T) + np.matmul((I - P).T, (I - P))

    return K


def W(verts):
    area, _, _ = poly_geom(verts)

    v = np.concatenate((verts, verts[0, :].reshape(-1, 1).T), axis=0)
    n = vectort(verts)

    n_x = n[:, 0]
    n_y = n[:, 1]

    n_xx = np.roll(n_x, verts.shape[0] + 1)
    n_yy = np.roll(n_y, verts.shape[0] + 1)

    l = np.zeros(shape=(verts.shape[0],))

    for i in range(verts.shape[0]):
        l[i] = np.linalg.norm(v[i, :] - v[i + 1, :])

    l_aux = np.roll(l, verts.shape[0] - 1)

    W = np.zeros(shape=(verts.shape[0], 3))

    for i in range(verts.shape[0]):
        W[i, :] = [
            1 / verts.shape[0],
            0.5 * (l_aux[i] * n_xx[i] + l[i] * n_x[i]) / area,
            0.5 * (l_aux[i] * n_yy[i] + l[i] * n_y[i]) / area
        ]

    return W


def H(verts):
    x, y = verts[:, 0], verts[:, 1]
    area, _, centroid = poly_geom(verts)

    H = np.zeros(shape=(verts.shape[0], 3))

    for i in range(verts.shape[0]):
        H[i, :] = [1, x[i] - centroid[0], y[i] - centroid[1]]

    return H


def shapefunction(x, y, verts):
    N = np.array([
        [0.25 * (x - 1) * (y - 1), -0.25 * (x + 1) * (y - 1), 0.25 * (x + 1) * (y + 1), -0.25 * (x - 1) * (y + 1)]
    ])
    dN = np.array([
        [0.25 * (y - 1), 0.25 * (x - 1)],
        [-0.25 * (y - 1), -0.25 * (x + 1)],
        [0.25 * (y + 1), 0.25 * (x + 1)],
        [-0.25 * (y + 1), -0.25 * (x - 1)]
    ])

    vx, vy = verts[:, 0], verts[:, 1]
    dN_x, dN_y = dN[:, 0], dN[:, 1]

    xx = np.dot(dN_x, vx)
    xy = np.dot(dN_x, vy)
    yx = np.dot(dN_y, vx)
    yy = np.dot(dN_y, vy)

    j = np.array([
        [xx, xy],
        [yx, yy]
    ])

    jacobian = np.linalg.det(j)
    col = np.zeros(shape=(dN.shape[0], 2))

    for i in range(dN.shape[0]):
        col[i, 0] = np.dot(dN[i, :], np.array([yy, -xy]))
        col[i, 1] = np.dot(dN[i, :], np.array([-yx, xx]))

    return N, dN, jacobian, col


def rhs(verts):
    x, y = verts[:, 0], verts[:, 1]
    f = -2 * (-1 - x) * (1 - x) - 2 * (-1 - y) * (1 - y)
    return f


def f_local(verts):
    X, Y = verts[:, 0], verts[:, 1]
    n = 20
    x = gausspoints1d(n)
    w = gaussweights1d(n)
    f_local = np.zeros(shape=(verts.shape[0], 1))

    for i in range(n):
        for j in range(n):
            N, _, jacobian, _ = shapefunction(x[i], x[j], verts)
            v = np.array([[np.dot(N, X), np.dot(N, Y)]])
            f = rhs(v)
            # print(f_local, N.T, w[i], w[j], jacobian, f)
            f_local = f_local + N.T * w[i] * w[j] * jacobian * f
    return f_local


def u_solution(verts):
    x, y = verts[:, 0], verts[:, 1]
    u = (-1 - x) * (-1 - y) * (1 - x) * (1 - y)
    return u


def K_e(verts):
    X, Y = verts[:, 0], verts[:, 1]
    n = 20
    x = gausspoints1d(n)
    w = gaussweights1d(n)
    k_elemental = np.zeros(shape=(verts.shape[0], verts.shape[0]))

    for i in range(n):
        for j in range(n):
            N, _, jacobian, col = shapefunction(x[i], x[j], verts)
            k_elemental = k_elemental + np.dot(col, col.T) * (1 / jacobian) * (1 / jacobian) * jacobian * w[i] * w[j]
    return k_elemental


def K_elasticity(verts):
    x, y = verts[:, 0], verts[:, 1]

    area, _, centroid = poly_geom(verts)

    v = np.concatenate((verts, verts[0, :].reshape(-1, 1).T), axis=0)
    n = vectort(verts)

    n_x = n[:, 0]
    n_y = n[:, 1]

    n_xx = np.roll(n_x, verts.shape[0] + 1)
    n_yy = np.roll(n_y, verts.shape[0] + 1)

    l = np.zeros(shape=(verts.shape[0],))

    for i in range(verts.shape[0]):
        l[i] = np.linalg.norm(v[i, :] - v[i + 1, :])

    l_aux = np.roll(l, verts.shape[0] - 1)

    Hr = np.zeros(shape=(2 * verts.shape[0], 3))
    Hc = np.zeros(shape=(2 * verts.shape[0], 3))

    Wr = np.zeros(shape=(2 * verts.shape[0], 3))
    Wc = np.zeros(shape=(2 * verts.shape[0], 3))

    for i in range(verts.shape[0]):
        Hr[[2 * i, 2 * i + 1], :] = np.array([
            [1, 0, y[i] - centroid[1]],
            [0, 1, -(x[i] - centroid[0])]
        ])

        Hc[[2 * i, 2 * i + 1], :] = np.array([
            [(x[i] - centroid[0]), 0, y[i] - centroid[1]],
            [0, (y[i] - centroid[1]), (x[i] - centroid[0])]
        ])

        Wr[[2 * i, 2 * i + 1], :] = np.array([
            [1 / verts.shape[0], 0, 0.5 * 0.5 * (l_aux[i] * n_yy[i] + l[i] * n_y[i]) / area],
            [0, 1 / verts.shape[0], -(0.5 * 0.5 * (l_aux[i] * n_xx[i] + l[i] * n_x[i]) / area)]
        ])

        Wc[[2 * i, 2 * i + 1], :] = np.array([
            [0.5 * (l_aux[i] * n_xx[i] + l[i] * n_x[i]) / area, 0,
             0.5 * 0.5 * (l_aux[i] * n_yy[i] + l[i] * n_y[i]) / area],
            [0, 0.5 * (l_aux[i] * n_yy[i] + l[i] * n_y[i]) / area,
             0.5 * 0.5 * (l_aux[i] * n_xx[i] + l[i] * n_x[i]) / area]
        ])

    Pc = np.matmul(Hc, Wc.T)
    Pr = np.matmul(Hr, Wr.T)
    Pp = Pr + Pc

    nu = 0.3
    E = 3 * 10 ** 7
    # D = E / ((1 + nu) * (1 - 2 * nu)) * np.array([[1 - nu, nu, 0], [nu, 1 - nu, 0], [0, 0, 2 * (1 - 2 * nu)]])
    D = E / (1 - nu ** 2) * np.array([[1, nu, 0], [nu, 1, 0], [0, 0, 2 * (1 - nu)]])
    gamma = 1

    alpha = gamma * area * np.trace(D) / np.trace(np.matmul(Hc.T, Hc))
    S = alpha * np.identity(2 * verts.shape[0])
    K = area * np.matmul(np.matmul(Wc, D), Wc.T) + np.matmul(np.matmul((np.identity(2 * verts.shape[0]) - Pp).T, S),
                                                             (np.identity(2 * verts.shape[0]) - Pp))
    return K


def f_elasticity(verts):
    x, y = verts[:, 0], verts[:, 1]
    area, _, _ = poly_geom(verts)
    N = np.zeros(shape=(2 * verts.shape[0], 2))

    for i in range(verts.shape[0]):
        N[[2 * i, 2 * i + 1], :] = np.array([[1 / verts.shape[0], 0], [0, 1 / verts.shape[0]]])
    f = area * np.matmul(N, np.array([[0], [1]]))
    return f


def gausspoints1d(n):
    x = np.ones(shape=(n,))

    if (n == 1):
        x[0] = 0
    elif (n == 2):
        x[0] = -0.577350269189625764509148780502
        x[1] = 0.577350269189625764509148780502
    elif (n == 3):
        x[0] = -0.774596669241483377035853079956
        x[1] = 0.000000000000000
        x[2] = 0.774596669241483377035853079956
    elif (n == 4):
        x[0] = -0.861136311594052575223946488893
        x[1] = -0.339981043584856264802665759103
        x[2] = 0.339981043584856264802665759103
        x[3] = 0.861136311594052575223946488893
    elif (n == 5):
        x[0] = -0.906179845938663992797626878299
        x[1] = -0.538469310105683091036314420700
        x[2] = 0.0
        x[3] = 0.538469310105683091036314420700
        x[4] = 0.906179845938663992797626878299
    elif (n == 6):
        x[0] = -0.932469514203152027812301554494
        x[1] = -0.661209386466264513661399595020
        x[2] = -0.238619186083196908630501721681
        x[3] = 0.238619186083196908630501721681
        x[4] = 0.661209386466264513661399595020
        x[5] = 0.932469514203152027812301554494
    elif (n == 7):
        x[0] = -0.949107912342758524526189684048
        x[1] = -0.741531185599394439863864773281
        x[2] = -0.405845151377397166906606412077
        x[3] = 0.0
        x[4] = 0.405845151377397166906606412077
        x[5] = 0.741531185599394439863864773281
        x[6] = 0.949107912342758524526189684048
    elif (n == 8):
        x[0] = -0.960289856497536231683560868569
        x[1] = -0.796666477413626739591553936476
        x[2] = -0.525532409916328985817739049189
        x[3] = -0.183434642495649804939476142360
        x[4] = 0.183434642495649804939476142360
        x[5] = 0.525532409916328985817739049189
        x[6] = 0.796666477413626739591553936476
        x[7] = 0.960289856497536231683560868569
    elif (n == 9):
        x[0] = -0.968160239507626089835576202904
        x[1] = -0.836031107326635794299429788070
        x[2] = -0.613371432700590397308702039341
        x[3] = -0.324253423403808929038538014643
        x[4] = 0.0
        x[5] = 0.324253423403808929038538014643
        x[6] = 0.613371432700590397308702039341
        x[7] = 0.836031107326635794299429788070
        x[8] = 0.968160239507626089835576202904
    elif (n == 10):
        x[0] = -0.973906528517171720077964012084
        x[1] = -0.865063366688984510732096688423
        x[2] = -0.679409568299024406234327365115
        x[3] = -0.433395394129247290799265943166
        x[4] = -0.148874338981631210884826001130
        x[5] = 0.148874338981631210884826001130
        x[6] = 0.433395394129247290799265943166
        x[7] = 0.679409568299024406234327365115
        x[8] = 0.865063366688984510732096688423
        x[9] = 0.973906528517171720077964012084
    elif (n == 11):
        x[0] = -0.978228658146056992803938001123
        x[1] = -0.887062599768095299075157769304
        x[2] = -0.730152005574049324093416252031
        x[3] = -0.519096129206811815925725669459
        x[4] = -0.269543155952344972331531985401
        x[5] = 0
        x[6] = 0.269543155952344972331531985401
        x[7] = 0.519096129206811815925725669459
        x[8] = 0.730152005574049324093416252031
        x[9] = 0.887062599768095299075157769304
        x[10] = 0.978228658146056992803938001123
    elif (n == 12):
        x[0] = -0.981560634246719250690549090149
        x[1] = -0.904117256370474856678465866119
        x[2] = -0.769902674194304687036893833213
        x[3] = -0.587317954286617447296702418941
        x[4] = -0.367831498998180193752691536644
        x[5] = -0.125233408511468915472441369464
        x[6] = 0.125233408511468915472441369464
        x[7] = 0.367831498998180193752691536644
        x[8] = 0.587317954286617447296702418941
        x[9] = 0.769902674194304687036893833213
        x[10] = 0.904117256370474856678465866119
        x[11] = 0.981560634246719250690549090149
    elif (n == 13):
        x[0] = -0.984183054718588149472829448807
        x[1] = -0.917598399222977965206547836501
        x[2] = -0.801578090733309912794206489583
        x[3] = -0.642349339440340220643984606996
        x[4] = -0.448492751036446852877912852128
        x[5] = -0.230458315955134794065528121098
        x[6] = 0
        x[7] = 0.230458315955134794065528121098
        x[8] = 0.448492751036446852877912852128
        x[9] = 0.642349339440340220643984606996
        x[10] = 0.801578090733309912794206489583
        x[11] = 0.917598399222977965206547836501
        x[12] = 0.984183054718588149472829448807
    elif (n == 14):
        x[0] = -0.986283808696812338841597266704
        x[1] = -0.928434883663573517336391139378
        x[2] = -0.827201315069764993189794742650
        x[3] = -0.687292904811685470148019803019
        x[4] = -0.515248636358154091965290718551
        x[5] = -0.319112368927889760435671824168
        x[6] = -0.108054948707343662066244650220
        x[7] = 0.108054948707343662066244650220
        x[8] = 0.319112368927889760435671824168
        x[9] = 0.515248636358154091965290718551
        x[10] = 0.687292904811685470148019803019
        x[11] = 0.827201315069764993189794742650
        x[12] = 0.928434883663573517336391139378
        x[13] = 0.986283808696812338841597266704
    elif (n == 15):
        x[0] = -0.987992518020485428489565718587
        x[1] = -0.937273392400705904307758947710
        x[2] = -0.848206583410427216200648320774
        x[3] = -0.724417731360170047416186054614
        x[4] = -0.570972172608538847537226737254
        x[5] = -0.394151347077563369897207370981
        x[6] = -0.201194093997434522300628303395
        x[7] = 0
        x[8] = 0.201194093997434522300628303395
        x[9] = 0.394151347077563369897207370981
        x[10] = 0.570972172608538847537226737254
        x[11] = 0.724417731360170047416186054614
        x[12] = 0.848206583410427216200648320774
        x[13] = 0.937273392400705904307758947710
        x[14] = 0.987992518020485428489565718587
    elif (n == 16):
        x[0] = -0.989400934991649932596154173450
        x[1] = -0.944575023073232576077988415535
        x[2] = -0.865631202387831743880467897712
        x[3] = -0.755404408355003033895101194847
        x[4] = -0.617876244402643748446671764049
        x[5] = -0.458016777657227386342419442984
        x[6] = -0.281603550779258913230460501460
        x[7] = -0.0950125098376374401853193354250
        x[8] = 0.0950125098376374401853193354250
        x[9] = 0.281603550779258913230460501460
        x[10] = 0.458016777657227386342419442984
        x[11] = 0.617876244402643748446671764049
        x[12] = 0.755404408355003033895101194847
        x[13] = 0.865631202387831743880467897712
        x[14] = 0.944575023073232576077988415535
        x[15] = 0.989400934991649932596154173450
    elif (n == 17):
        x[0] = -0.990575475314417335675434019941
        x[1] = -0.950675521768767761222716957896
        x[2] = -0.880239153726985902122955694488
        x[3] = -0.781514003896801406925230055520
        x[4] = -0.657671159216690765850302216643
        x[5] = -0.512690537086476967886246568630
        x[6] = -0.351231763453876315297185517095
        x[7] = -0.178484181495847855850677493654
        x[8] = 0
        x[9] = 0.178484181495847855850677493654
        x[10] = 0.351231763453876315297185517095
        x[11] = 0.512690537086476967886246568630
        x[12] = 0.657671159216690765850302216643
        x[13] = 0.781514003896801406925230055520
        x[14] = 0.880239153726985902122955694488
        x[15] = 0.950675521768767761222716957896
        x[16] = 0.990575475314417335675434019941
    elif (n == 18):
        x[0] = -0.991565168420930946730016004706
        x[1] = -0.955823949571397755181195892930
        x[2] = -0.892602466497555739206060591127
        x[3] = -0.803704958972523115682417455015
        x[4] = -0.691687043060353207874891081289
        x[5] = -0.559770831073947534607871548525
        x[6] = -0.411751161462842646035931793833
        x[7] = -0.251886225691505509588972854878
        x[8] = -0.0847750130417353012422618529358
        x[9] = 0.0847750130417353012422618529358
        x[10] = 0.251886225691505509588972854878
        x[11] = 0.411751161462842646035931793833
        x[12] = 0.559770831073947534607871548525
        x[13] = 0.691687043060353207874891081289
        x[14] = 0.803704958972523115682417455015
        x[15] = 0.892602466497555739206060591127
        x[16] = 0.955823949571397755181195892930
        x[17] = 0.991565168420930946730016004706
    elif (n == 19):
        x[0] = -0.992406843843584403189017670253
        x[1] = -0.960208152134830030852778840688
        x[2] = -0.903155903614817901642660928532
        x[3] = -0.822714656537142824978922486713
        x[4] = -0.720966177335229378617095860824
        x[5] = -0.600545304661681023469638164946
        x[6] = -0.464570741375960945717267148104
        x[7] = -0.316564099963629831990117328850
        x[8] = -0.160358645640225375868096115741
        x[9] = 0
        x[10] = 0.160358645640225375868096115741
        x[11] = 0.316564099963629831990117328850
        x[12] = 0.464570741375960945717267148104
        x[13] = 0.600545304661681023469638164946
        x[14] = 0.720966177335229378617095860824
        x[15] = 0.822714656537142824978922486713
        x[16] = 0.903155903614817901642660928532
        x[17] = 0.960208152134830030852778840688
        x[18] = 0.992406843843584403189017670253
    elif (n == 20):
        x[0] = -0.993128599185094924786122388471
        x[1] = -0.963971927277913791267666131197
        x[2] = -0.912234428251325905867752441203
        x[3] = -0.839116971822218823394529061702
        x[4] = -0.746331906460150792614305070356
        x[5] = -0.636053680726515025452836696226
        x[6] = -0.510867001950827098004364050955
        x[7] = -0.373706088715419560672548177025
        x[8] = -0.227785851141645078080496195369
        x[9] = -0.0765265211334973337546404093988
        x[10] = 0.0765265211334973337546404093988
        x[11] = 0.227785851141645078080496195369
        x[12] = 0.373706088715419560672548177025
        x[13] = 0.510867001950827098004364050955
        x[14] = 0.636053680726515025452836696226
        x[15] = 0.746331906460150792614305070356
        x[16] = 0.839116971822218823394529061702
        x[17] = 0.912234428251325905867752441203
        x[18] = 0.963971927277913791267666131197
        x[19] = 0.993128599185094924786122388471
    else:
        print('GAUSS_POINTS - Fatal error! Illegal value of n.')
    return x


def gaussweights1d(n):
    w = np.ones(shape=(n,))

    if (n == 1):
        w[0] = 2.0
    elif (n == 2):
        w[0] = 1.0
        w[1] = w[0]
    elif (n == 3):
        w[0] = 0.5555555555555555555555555555565
        w[1] = 0.8888888888888888888888888888889
        w[2] = 0.5555555555555555555555555555565
    elif (n == 4):
        w[0] = 0.347854845137453857373063949222
        w[1] = 0.652145154862546142626936050778
        w[2] = 0.652145154862546142626936050778
        w[3] = 0.347854845137453857373063949222
    elif (n == 5):
        w[0] = 0.236926885056189087514264040720
        w[1] = 0.478628670499366468041291514836
        w[2] = 0.568888888888888888888888888889
        w[3] = 0.478628670499366468041291514836
        w[4] = 0.236926885056189087514264040720
    elif (n == 6):
        w[0] = 0.171324492379170345040296142173
        w[1] = 0.360761573048138607569833513838
        w[2] = 0.467913934572691047389870343990
        w[3] = 0.467913934572691047389870343990
        w[4] = 0.360761573048138607569833513838
        w[5] = 0.171324492379170345040296142173
    elif (n == 7):
        w[0] = 0.129484966168869693270611432679
        w[1] = 0.279705391489276667901467771424
        w[2] = 0.381830050505118944950369775489
        w[3] = 0.417959183673469387755102040816
        w[4] = 0.381830050505118944950369775489
        w[5] = 0.279705391489276667901467771424
        w[6] = 0.129484966168869693270611432679
    elif (n == 8):
        w[0] = 0.101228536290376259152531354310
        w[1] = 0.222381034453374470544355994426
        w[2] = 0.313706645877887287337962201987
        w[3] = 0.362683783378361982965150449277
        w[4] = 0.362683783378361982965150449277
        w[5] = 0.313706645877887287337962201987
        w[6] = 0.222381034453374470544355994426
        w[7] = 0.101228536290376259152531354310
    elif (n == 9):
        w[0] = 0.0812743883615744119718921581105
        w[1] = 0.180648160694857404058472031243
        w[2] = 0.260610696402935462318742869419
        w[3] = 0.312347077040002840068630406584
        w[4] = 0.330239355001259763164525069287
        w[5] = 0.312347077040002840068630406584
        w[6] = 0.260610696402935462318742869419
        w[7] = 0.180648160694857404058472031243
        w[8] = 0.0812743883615744119718921581105
    elif (n == 10):
        w[0] = 0.0666713443086881375935688098933
        w[1] = 0.149451349150580593145776339658
        w[2] = 0.219086362515982043995534934228
        w[3] = 0.269266719309996355091226921569
        w[4] = 0.295524224714752870173892994651
        w[5] = 0.295524224714752870173892994651
        w[6] = 0.269266719309996355091226921569
        w[7] = 0.219086362515982043995534934228
        w[8] = 0.149451349150580593145776339658
        w[9] = 0.0666713443086881375935688098933
    elif (n == 11):
        w[0] = 0.0556685671161736664827537204425
        w[1] = 0.125580369464904624634694299224
        w[2] = 0.186290210927734251426097641432
        w[3] = 0.233193764591990479918523704843
        w[4] = 0.262804544510246662180688869891
        w[5] = 0.272925086777900630714483528336
        w[6] = 0.262804544510246662180688869891
        w[7] = 0.233193764591990479918523704843
        w[8] = 0.186290210927734251426097641432
        w[9] = 0.125580369464904624634694299224
        w[10] = 0.0556685671161736664827537204425
    elif (n == 12):
        w[0] = 0.0471753363865118271946159614850
        w[1] = 0.106939325995318430960254718194
        w[2] = 0.160078328543346226334652529543
        w[3] = 0.203167426723065921749064455810
        w[4] = 0.233492536538354808760849898925
        w[5] = 0.249147045813402785000562436043
        w[6] = 0.249147045813402785000562436043
        w[7] = 0.233492536538354808760849898925
        w[8] = 0.203167426723065921749064455810
        w[9] = 0.160078328543346226334652529543
        w[10] = 0.106939325995318430960254718194
        w[11] = 0.0471753363865118271946159614850
    elif (n == 13):
        w[0] = 0.0404840047653158795200215922010
        w[1] = 0.0921214998377284479144217759538
        w[2] = 0.138873510219787238463601776869
        w[3] = 0.178145980761945738280046691996
        w[4] = 0.207816047536888502312523219306
        w[5] = 0.226283180262897238412090186040
        w[6] = 0.232551553230873910194589515269
        w[7] = 0.226283180262897238412090186040
        w[8] = 0.207816047536888502312523219306
        w[9] = 0.178145980761945738280046691996
        w[10] = 0.138873510219787238463601776869
        w[11] = 0.0921214998377284479144217759538
        w[12] = 0.0404840047653158795200215922010
    elif (n == 14):
        w[0] = 0.0351194603317518630318328761382
        w[1] = 0.0801580871597602098056332770629
        w[2] = 0.121518570687903184689414809072
        w[3] = 0.157203167158193534569601938624
        w[4] = 0.185538397477937813741716590125
        w[5] = 0.205198463721295603965924065661
        w[6] = 0.215263853463157790195876443316
        w[7] = 0.215263853463157790195876443316
        w[8] = 0.205198463721295603965924065661
        w[9] = 0.185538397477937813741716590125
        w[10] = 0.157203167158193534569601938624
        w[11] = 0.121518570687903184689414809072
        w[12] = 0.0801580871597602098056332770629
        w[13] = 0.0351194603317518630318328761382
    elif (n == 15):
        w[0] = 0.0307532419961172683546283935772
        w[1] = 0.0703660474881081247092674164507
        w[2] = 0.107159220467171935011869546686
        w[3] = 0.139570677926154314447804794511
        w[4] = 0.166269205816993933553200860481
        w[5] = 0.186161000015562211026800561866
        w[6] = 0.198431485327111576456118326444
        w[7] = 0.202578241925561272880620199968
        w[8] = 0.198431485327111576456118326444
        w[9] = 0.186161000015562211026800561866
        w[10] = 0.166269205816993933553200860481
        w[11] = 0.139570677926154314447804794511
        w[12] = 0.107159220467171935011869546686
        w[13] = 0.0703660474881081247092674164507
        w[14] = 0.0307532419961172683546283935772
    elif (n == 16):
        w[0] = 0.0271524594117540948517805724560
        w[1] = 0.0622535239386478928628438369944
        w[2] = 0.0951585116824927848099251076022
        w[3] = 0.124628971255533872052476282192
        w[4] = 0.149595988816576732081501730547
        w[5] = 0.169156519395002538189312079030
        w[6] = 0.182603415044923588866763667969
        w[7] = 0.189450610455068496285396723208
        w[8] = 0.189450610455068496285396723208
        w[9] = 0.182603415044923588866763667969
        w[10] = 0.169156519395002538189312079030
        w[11] = 0.149595988816576732081501730547
        w[14] = 0.124628971255533872052476282192
        w[15] = 0.0951585116824927848099251076022
        w[16] = 0.0622535239386478928628438369944
        w[17] = 0.0271524594117540948517805724560
    elif (n == 17):
        w[0] = 0.0241483028685479319601100262876
        w[1] = 0.0554595293739872011294401653582
        w[2] = 0.0850361483171791808835353701911
        w[3] = 0.111883847193403971094788385626
        w[4] = 0.135136368468525473286319981702
        w[5] = 0.154045761076810288081431594802
        w[6] = 0.168004102156450044509970663788
        w[7] = 0.176562705366992646325270990113
        w[8] = 0.179446470356206525458265644262
        w[9] = 0.176562705366992646325270990113
        w[10] = 0.168004102156450044509970663788
        w[11] = 0.154045761076810288081431594802
        w[12] = 0.135136368468525473286319981702
        w[13] = 0.111883847193403971094788385626
        w[14] = 0.0850361483171791808835353701911
        w[15] = 0.0554595293739872011294401653582
        w[16] = 0.0241483028685479319601100262876
    elif (n == 18):
        w[0] = 0.0216160135264833103133427102665
        w[1] = 0.0497145488949697964533349462026
        w[2] = 0.0764257302548890565291296776166
        w[3] = 0.100942044106287165562813984925
        w[4] = 0.122555206711478460184519126800
        w[5] = 0.140642914670650651204731303752
        w[6] = 0.154684675126265244925418003836
        w[7] = 0.164276483745832722986053776466
        w[8] = 0.169142382963143591840656470135
        w[9] = 0.169142382963143591840656470135
        w[10] = 0.164276483745832722986053776466
        w[11] = 0.154684675126265244925418003836
        w[12] = 0.140642914670650651204731303752
        w[13] = 0.122555206711478460184519126800
        w[14] = 0.100942044106287165562813984925
        w[15] = 0.0764257302548890565291296776166
        w[16] = 0.0497145488949697964533349462026
        w[17] = 0.0216160135264833103133427102665
    elif (n == 19):
        w[0] = 0.0194617882297264770363120414644
        w[1] = 0.0448142267656996003328381574020
        w[2] = 0.0690445427376412265807082580060
        w[3] = 0.0914900216224499994644620941238
        w[4] = 0.111566645547333994716023901682
        w[5] = 0.128753962539336227675515784857
        w[6] = 0.142606702173606611775746109442
        w[7] = 0.152766042065859666778855400898
        w[8] = 0.158968843393954347649956439465
        w[9] = 0.161054449848783695979163625321
        w[10] = 0.158968843393954347649956439465
        w[11] = 0.152766042065859666778855400898
        w[12] = 0.142606702173606611775746109442
        w[13] = 0.128753962539336227675515784857
        w[14] = 0.111566645547333994716023901682
        w[15] = 0.0914900216224499994644620941238
        w[16] = 0.0690445427376412265807082580060
        w[17] = 0.0448142267656996003328381574020
        w[18] = 0.0194617882297264770363120414644
    elif (n == 20):
        w[0] = 0.0176140071391521183118619623519
        w[1] = 0.0406014298003869413310399522749
        w[2] = 0.0626720483341090635695065351870
        w[3] = 0.0832767415767047487247581432220
        w[4] = 0.101930119817240435036750135480
        w[5] = 0.118194531961518417312377377711
        w[6] = 0.131688638449176626898494499748
        w[7] = 0.142096109318382051329298325067
        w[8] = 0.149172986472603746787828737002
        w[9] = 0.152753387130725850698084331955
        w[10] = 0.152753387130725850698084331955
        w[11] = 0.149172986472603746787828737002
        w[12] = 0.142096109318382051329298325067
        w[13] = 0.131688638449176626898494499748
        w[14] = 0.118194531961518417312377377711
        w[15] = 0.101930119817240435036750135480
        w[16] = 0.0832767415767047487247581432220
        w[17] = 0.0626720483341090635695065351870
        w[18] = 0.0406014298003869413310399522749
        w[19] = 0.0176140071391521183118619623519
    else:
        print('GAUSS_WEIGHTS - Fatal error! Illegal value of n.')
    return w