In [ ]:

    

%matplotlib inline


    

In [ ]:

    

import numpy as np
from numpy import newaxis
import scipy.sparse as sps
from scipy.sparse.linalg import spsolve
import matplotlib.pyplot as plt


    

In [ ]:

    

from pyfem.topo import SQuad
from pyfem.basis import LagrangeBasisQuad

topo  = SQuad()
basis = LagrangeBasisQuad(topo, 1)


    

In [ ]:

    

def eval_basis(ref):
    
    r = np.zeros((4, len(ref), 2))
    r[0,:,0] = (1-ref[:,1])
    r[1,:,1] = (1+ref[:,0])
    r[2,:,0] = (1+ref[:,1])
    r[3,:,1] = (1-ref[:,0])
    
    return r/2.0


    

In [ ]:

    

xv = np.linspace(-1, 1, 10)
yv = np.linspace(-1, 1, 10)
X, Y = np.meshgrid(xv, yv)
X = X.ravel()
Y = Y.ravel()
ref = np.hstack([X[:,newaxis], Y[:,newaxis]])
r = eval_basis(ref)


    

In [ ]:

    

plt.figure(figsize=(8,8))
for i in range(len(r)):
    plt.subplot(2, 2, i+1)
    plt.quiver(X, Y, r[i,:,0], r[i,:,1], 
               color='b')
    plt.title(str(i))


    

In [ ]:

    

Kloc = np.array([[1,-1,-1,1],
                 [-1,1,1,-1],
                 [-1,1,1,-1],
                 [1,-1,-1,1]], dtype=np.double)
Mloc = np.array([[2,1,0,0],
                 [1,2,0,0],
                 [0,0,2,1],
                 [0,0,1,2]], dtype=np.double)
Mloc *= 2.0/3.0


    

In [ ]:

    

N = 2
n_dofs = 2*N*(N+1)

etd = np.zeros((N*N, 4), dtype=np.int)
ielem = 0
for iy in range(N):
    for ix in range(N):
        etd[ielem,0] = ix+(2*N+1)*iy
        etd[ielem,1] = ix+(2*N+1)*iy+N+1
        etd[ielem,2] = ix+(2*N+1)*iy+2*N+1
        etd[ielem,3] = ix+(2*N+1)*iy+N
        ielem += 1


    

In [ ]:

    

cols = etd.ravel()
rows = np.arange(len(cols))
vals = np.ones(len(cols))
Q = sps.coo_matrix((vals, (rows, cols))).tocsr()


    

In [ ]:

    

A = sps.kron(np.eye(N*N), Mloc)
A = Q.T.dot(A.dot(Q))


    

In [ ]:

    

plt.spy(A)


    

In [ ]:

    

x, w = topo.get_quadrature(2)

dphi = basis.eval_ref(np.eye(4), x, d=1)
rhs = Q.T.dot(np.sum(eval_basis(x)*dphi, axis=-1).ravel())

spsolve(A, rhs)


    

In [ ]: