surface_term_xy_8LGL.sagews

%typeset_mode False

var('xi, eta')
# x_nodes = [-0.9, -1, 1, 1.1]%typeset_mode True

(xi, eta)

var('xi, eta')
#x_nodes = [0.6, 0.6, 0.8, 0.8]
#y_nodes = [0.6, 0.4, 0.4, 0.6]
x_nodes = [0,   0, 0.2, 0.2]
y_nodes = [0.2, 0, 0. , 0.2]

(xi, eta)

#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.1 0.1*eta + 0.1

x.simplify_full()
y.simplify_full()

0.1*xi + 0.1 0.1*eta + 0.1

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.

u(xi=-1)

e^(-6.25000000000000*(0.100000000000000*eta + 0.100000000000000)^2)

numerical_integral(L_xi[0](-1)* L_eta[1] * u(xi=-1), -1, 1)

(0.21048769390920383, 3.780673341015353e-15)

/cocalc/lib/python2.7/site-packages/smc_sagews/sage_server.py:1009: DeprecationWarning: Substitution using function-call syntax and unnamed arguments is deprecated and will be removed from a future release of Sage; you can use named arguments instead, like EXPR(x=..., y=...) See http://trac.sagemath.org/5930 for details. exec compile(block+'\n', '', 'single') in namespace, locals

numerical_integral(  L_xi[7](-1)* L_eta[6]     * u(xi=-1)   ,  -1, 1)[0]\
- numerical_integral(L_xi[7]    * L_eta[6](-1) * u(eta = -1),  -1, 1)[0]\
+ numerical_integral(L_xi[7](1) * L_eta[6]     * u(xi=1)    ,  -1, 1)[0]\
- numerical_integral(L_xi[7]    * L_eta[6](1)  * u(eta = 1) ,  -1, 1)[0]

0.0 0.0 0.1318272334786214 -0.0

0.21048769390920383 - 0.21048770042818485

-6.518981019928982e-9

1
0.00276958627370 - 0.00276955959998304

1 2.66737169600170e-8

integrate(L_xi[0](-1) * L_eta[0] * u(xi=-1), eta, -1, 1).n()

Error in lines 1-1 Traceback (most recent call last): File "/cocalc/lib/python2.7/site-packages/smc_sagews/sage_server.py", line 1009, in execute exec compile(block+'\n', '', 'single') in namespace, locals File "", line 1, in <module> File "/ext/sage/sage-8.0/local/lib/python2.7/site-packages/sage/misc/functional.py", line 721, in integral return x.integral(*args, **kwds) File "sage/symbolic/expression.pyx", line 12275, in sage.symbolic.expression.Expression.integral (/ext/sage/sage-8.0/src/build/cythonized/sage/symbolic/expression.cpp:69945) return integral(self, *args, **kwds) File "/ext/sage/sage-8.0/local/lib/python2.7/site-packages/sage/symbolic/integration/integral.py", line 799, in integrate return definite_integral(expression, v, a, b, hold=hold) File "sage/symbolic/function.pyx", line 996, in sage.symbolic.function.BuiltinFunction.__call__ (/ext/sage/sage-8.0/src/build/cythonized/sage/symbolic/function.cpp:11397) res = super(BuiltinFunction, self).__call__( File "sage/symbolic/function.pyx", line 486, in sage.symbolic.function.Function.__call__ (/ext/sage/sage-8.0/src/build/cythonized/sage/symbolic/function.cpp:6440) res = g_function_evalv(self._serial, vec, hold) File "sage/symbolic/function.pyx", line 1085, in sage.symbolic.function.BuiltinFunction._evalf_or_eval_ (/ext/sage/sage-8.0/src/build/cythonized/sage/symbolic/function.cpp:12757) return self._eval0_(*args) File "/ext/sage/sage-8.0/local/lib/python2.7/site-packages/sage/symbolic/integration/integral.py", line 181, in _eval_ return integrator(*args) File "/ext/sage/sage-8.0/local/lib/python2.7/site-packages/sage/symbolic/integration/external.py", line 34, in maxima_integrator result = maxima.sr_integral(expression, v, a, b) File "/ext/sage/sage-8.0/local/lib/python2.7/site-packages/sage/interfaces/maxima_lib.py", line 811, in sr_integral raise ValueError("Integral is divergent.") ValueError: Integral is divergent.

numerical_integral(L_xi[0] * L_eta[0](-1) * u(eta = -1), -1, 1)

(0.03571428835994642, 1.1856233645472606e-15)

0.00138478939334163 + 0.0013847931808617447

0.00276958257420337

###
integrate(integrate(L_xi[4] * L_eta[5] * L_xi[6] * L_eta[3], eta, -1, 1), xi, -1, 1).n()

0.000491462040977236

0.000491462040977236 - 0.0004914620419140329

-9.36796972460241e-13