Contact
CoCalc Logo Icon
StoreFeaturesDocsShareSupport News AboutSign UpSign In
| Download
Views: 1021
%typeset_mode False
var('xi, eta') # x_nodes = [-0.9, -1, 1, 1.1]%typeset_mode True
(xi, eta)
# Element 1 x_nodes = [-1, -1, -0.8, -0.8] y_nodes = [1, 0.8, 0.8, 1] # Element 28 #x_nodes = [0.6, 0.6, 0.8, 0.8] #y_nodes = [0.6, 0.4, 0.4, 0.6] # Element 45 #x_nodes = [0, 0, 0.2, 0.2] #y_nodes = [0.2, 0, 0, 0.2]
#N_0 = 0.25 * (1 - xi) * (1 - eta) #Aman #N_1 = 0.25 * (1 + xi) * (1 - eta) #N_2 = 0.25 * (1 + xi) * (1 + eta) #N_3 = 0.25 * (1 - xi) * (1 + eta) N_0 = +(1/4) * (1-xi) * (1+eta) # Balavarun N_1 = +(1/4) * (1-eta) * (1-xi) N_2 = +(1/4) * (1+xi) * (1-eta) N_3 = +(1/4) * (1+xi) * (1+eta) x = N_0 * x_nodes[0] + N_1 * x_nodes[1] + N_2 * x_nodes[2] + N_3 * x_nodes[3] y = N_0 * y_nodes[0] + N_1 * y_nodes[1] + N_2 * y_nodes[2] + N_3 * y_nodes[3]
x.simplify_full() y.simplify_full() ## y_nodes = [1.1, -1.0, -1.2, 1.4]
0.1*xi - 0.8999999999999999 0.1*eta + 0.9
x.simplify_full() y.simplify_full()
0.1*xi - 0.8999999999999999 0.1*eta + 0.9
import numpy as np 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] 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.
print(u)
e^(-6.25000000000000*(1/4*(eta + 1)*(xi + 1) - 0.200000000000000*(eta - 1)*(xi + 1) - 1/4*(eta + 1)*(xi - 1) + 0.200000000000000*(eta - 1)*(xi - 1))^2 - 6.25000000000000*(-0.200000000000000*(eta + 1)*(xi + 1) + 0.200000000000000*(eta - 1)*(xi + 1) + 1/4*(eta + 1)*(xi - 1) - 1/4*(eta - 1)*(xi - 1))^2)
integrate(integrate(u * diff(L_xi[0]) * L_eta[0],xi, -1, 1), eta, -1, 1).n()
-0.267795562744141
numerical_integral(integrate(u * diff(L_xi[0]) * L_eta[0], xi, -1, 1), -1, 1)[0]
-1.3191396607310232e-06
numerical_integral(integrate(u * diff(L_eta[0]) * L_xi[0], eta, -1, 1), -1, 1)[0]
-1.2176642079540523e-06
expand(diff(L_xi[0]) * L_eta[0])
78.6307983402193*eta^7*xi^6 - 67.3978271487594*eta^7*xi^5 - 78.6307983402193*eta^6*xi^6 - 64.8056030276711*eta^7*xi^4 + 67.3978271487594*eta^6*xi^5 - 90.7278442387395*eta^5*xi^6 + 51.8444824221369*eta^7*xi^3 + 64.8056030276711*eta^6*xi^4 + 77.7667236332053*eta^5*xi^5 + 90.7278442387395*eta^4*xi^6 + 10.6045532227039*eta^7*xi^2 - 51.8444824221369*eta^6*xi^3 + 74.7756958011795*eta^5*xi^4 - 77.7667236332053*eta^4*xi^5 + 24.7439575196426*eta^3*xi^6 - 7.06970214846930*eta^7*xi - 10.6045532227039*eta^6*xi^2 - 59.8205566409436*eta^5*xi^3 - 74.7756958011795*eta^4*xi^4 - 21.2091064454079*eta^3*xi^5 - 24.7439575196425*eta^2*xi^6 - 0.130920410156332*eta^7 + 7.06970214846930*eta^6*xi - 12.2360229492771*eta^5*xi^2 + 59.8205566409436*eta^4*xi^3 - 20.3933715821286*eta^3*xi^4 + 21.2091064454079*eta^2*xi^5 - 0.916442871094318*eta*xi^6 + 0.130920410156332*eta^6 + 8.15734863285144*eta^5*xi + 12.2360229492771*eta^4*xi^2 + 16.3146972657029*eta^3*xi^3 + 20.3933715821286*eta^2*xi^4 + 0.785522460937988*eta*xi^5 + 0.916442871094320*xi^6 + 0.151062011718886*eta^5 - 8.15734863285143*eta^4*xi + 3.33709716798283*eta^3*xi^2 - 16.3146972657029*eta^2*xi^3 + 0.755310058594427*eta*xi^4 - 0.785522460937989*xi^5 - 0.151062011718886*eta^4 - 2.22473144532188*eta^3*xi - 3.33709716798283*eta^2*xi^2 - 0.604248046875541*eta*xi^3 - 0.755310058594427*xi^4 - 0.0411987304687642*eta^3 + 2.22473144532189*eta^2*xi - 0.123596191406292*eta*xi^2 + 0.604248046875542*xi^3 + 0.0411987304687642*eta^2 + 0.0823974609375284*eta*xi + 0.123596191406292*xi^2 + 0.00152587890624461*eta - 0.0823974609375283*xi - 0.00152587890624461
expand(diff(L_xi[0]))
-23.4609375000560*xi^6 + 20.1093750000480*xi^5 + 19.3359375000515*xi^4 - 15.4687500000412*xi^3 - 3.16406250000667*xi^2 + 2.10937500000445*xi + 0.0390624999999311
for p in range(8): for q in range(8): print(numerical_integral(integrate(L_eta[q] * diff(L_xi[p]) * u, xi, -1, 1), -1, 1)[0] + numerical_integral(integrate(L_xi[p] * diff(L_eta[q]) * u, eta, -1, 1), -1, 1)[0])
-2.53680386869e-06 -6.60113341397e-06 -7.99132993168e-06 -6.39487512641e-06 -3.98570369608e-06 -2.09965123199e-06 -9.21345582622e-07 -2.40374284084e-14 -8.8094851747e-06 -3.52479801972e-07 -3.58609672579e-07 -2.12322200988e-07 -8.12271416177e-08 -1.8113011919e-08 -5.29395592034e-23 9.21345510377e-07 -2.0026039577e-05 -6.73539082638e-07 -6.60193413324e-07 -3.58198610872e-07 -1.06785708533e-07 -1.7663670819e-11 1.81055032678e-08 2.09964958812e-06 -3.7885823472e-05 -9.42836800787e-07 -8.46801363541e-07 -3.54718127537e-07 0.0 1.06754020668e-07 8.12271416177e-08 3.98570343104e-06 -6.05602583461e-05 -9.25247572284e-07 -6.47439549619e-07 -1.69406589451e-21 3.54718127537e-07 3.58146281385e-07 2.12322200988e-07 6.39487417064e-06 -7.54220521821e-05 -4.87914884422e-07 6.29249592292e-13 6.47439504418e-07 8.46801265745e-07 6.60131312381e-07 3.58609865075e-07 7.99132872838e-06 -6.2146522585e-05 0.0 4.87916111999e-07 9.25247572284e-07 9.42836800787e-07 6.73487733733e-07 3.52479801972e-07 6.60113333447e-06 1.69406589451e-21 6.2146522585e-05 7.54220515598e-05 6.05602583461e-05 3.7885823472e-05 2.00260297248e-05 8.8094851747e-06 2.53680338748e-06
-2.0026042988132806e-05 + 2.0026039577e-05 # Element no 1
-3.41113280606957e-12
-0.01017662293314848 + 0.0101766210904
-1.84274847984944e-9
(1.21766420795×1006\displaystyle -1.21766420795 \times 10^{-06}, 2.44857854363×1017\displaystyle 2.44857854363 \times 10^{-17})
︠5d9ec428-8707-4251-94ec-e86d66bbe516︠ 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 = 0 for weight, node in zip(gauss_weights, gauss_nodes): integral += weight * numerical_integral(diff(L_xi[p], xi) * F_xi * L_eta[q], -1, 1, params = [node])[0] print integral