In [47]:

    

import numpy as np
from matplotlib import pylab as plt
%matplotlib inline

def u(x,y):
    return (1-x)*(1-np.square(y))

def f(x,y):
    term1 = (x-1)*np.square(x)*(np.square(y)-1)*np.sin(x*y)
    term2 = (x-1)*np.square(y)*(np.square(y)-1)*np.sin(x*y)
    term3 = 2*y*(np.square(y)-1)*np.cos(x*y) - 2*(x-1)*np.sin(x*y)
    term4 = 4*(x-1)*x*y*np.cos(x*y)
    return term1 + term2 - term3 - term4

def N_2_1_1d(xi):
    return -xi*(1 - xi)/2

def N_2_2_1d(xi):
    return xi*(1 + xi)/2

def N_2_3_1d(xi):
    return np.square(1-xi)


#2d lagrange
def N_11(xi, eta):
    return N_2_1_1d(xi)*N_2_1_1d(eta)

def N_12(xi, eta):
    return N_2_1_1d(xi)*N_2_2_1d(eta)

def N_13(xi, eta):
    return N_2_1_1d(xi)*N_2_3_1d(eta)

def N_21(xi, eta):
    return N_2_2_1d(xi)*N_2_1_1d(eta)

def N_22(xi, eta):
    return N_2_2_1d(xi)*N_2_2_1d(eta)

def N_23(xi, eta):
    return N_2_2_1d(xi)*N_2_3_1d(eta)

def N_31(xi, eta):
    return N_2_3_1d(xi)*N_2_1_1d(eta)

def N_32(xi, eta):
    return N_2_3_1d(xi)*N_2_2_1d(eta)

def N_33(xi, eta):
    return N_2_3_1d(xi)*N_2_3_1d(eta)

def U_can(xi,eta, c):
    return 
        c[0]*N_11(xi, eta)+
        c[1]*N_12(xi, eta)+
        c[2]*N_13(xi, eta)+
        c[3]*N_21(xi, eta)+
        c[4]*N_22(xi, eta)+
        c[5]*N_23(xi, eta)+
        c[6]*N_31(xi, eta)+
        c[7]*N_32(xi, eta)+
        c[8]*N_33(xi, eta)
        
def N_11_1d(xi):
    return (1 - xi)/2

def N_12_1d(xi):
    return (1 + xi)/2

def map_can_to_xy(xi,eta, borders):
    # borders is a list with [x_ij,y_ij] lists. Each small list is a corner of a square
    # (this func is generalised to rectangles)
    # numeration: 11 is 0, 12 is 1, 21 is 2, 22 is 3 !!!! remember
    borders = np.array(borders)
    res = borders[0]*N_11_1d(xi)*N_11_1d(eta)
    res += borders[1]*N_11_1d(xi)*N_12_1d(eta)
    res += borders[2]*N_12_1d(xi)*N_11_1d(eta)
    res += borders[3]*N_12_1d(xi)*N_12_1d(eta)
    return res


    

In [35]:

    

def grid_vect(N):
    x = np.linspace(0, 1, N, endpoint=True)
    h = x[1] - x[0]
    return np.mgrid[0:1+h:h, 0:1+h:h].reshape(2,-1).T


    

In [46]:

    

xy = grid_vect(6)
#print(xy)
fig = plt.figure(figsize=(10,5))
plt.grid(True)
plt.scatter(xy[:,0],xy[:,1], marker='x', color = 'r', s = 50)


    

    
Out[46]:





<matplotlib.collections.PathCollection at 0x7f546554b630>






    

In [ ]:

    

def K_j():
    return 1
def I_j():
    return 0


    

In [ ]:

    




    

In [ ]:

    




    

In [ ]:

    




    

In [ ]:

    




    

In [ ]:

    




    

In [ ]:

    




    

In [ ]:

    




    

In [ ]:

    




    

In [ ]:

    




    

In [ ]:

    




    

In [ ]: