import matplotlib.pyplot as plt
import numpy as np

# generate mesh of nodes and elements
# we do this here by hand
# indexing starts at zero
node  = [ [0,0],
          [0,2],
          [1,1],
          [2,0],
          [2,2]]

el = [[0,2,1],     # el is the element nodal list for each element
      [0,3,2],
      [2,3,4],
      [1,2,4]]

x = [1,2,3]
y = [1,2,3]
Aj = np.zeros(len(el))
n=0
for j in el:            # j is the cord list for element j
    for i in range(3):  # triangle elements have 3 nodes
        x[i] = node[j[i]][0]  # the x,y coords for the corners
        y[i] = node[j[i]][1]
    # calculate the element areas, using coordinates
    Aj[n] = 0.5*np.abs(x[0]*(y[1]-y[2])+x[1]*(y[2]-y[0])+x[2]*(y[0]-y[1]))
    n +=1


# now make A matrix
nodes = len(node)
A = np.zeros([nodes,nodes])

for i in range(nodes):
    for j in range(nodes):
        # look for elements that have both i and j
        n=0
        for k in el:
            if (i == j) and (i in k):  # these are the integrals on the diagonal
                p = k.index(i)
                a = node[k[p-1]][1]-node[k[p-2]][1]
                b = node[k[p-2]][0]-node[k[p-1]][0]
                A[i,i] += (a*a + b*b)/(4*Aj[n])
            else:                      # the off diagonal terms
                if (i in k) and (j in k):
                    p = k.index(i)
                    q = k.index(j)
                    for m in k:        # need to find the node that is NOT i or j
                        if (k.index(m) != p) and (k.index(m) != q):
                            s = k.index(m)
                    a = node[k[p]][1]-node[k[s]][1]
                    b = node[k[s]][1]-node[k[q]][1]
                    c = node[k[p]][0]-node[k[s]][0]
                    d = node[k[s]][0]-node[k[q]][0]
                    A[i,j] += (a*b + c*d)/(4*Aj[n])
            n +=1