In [62]:
import os, sys
lib_path = os.path.abspath('../')
sys.path.append(lib_path)
from interfaces import *
from utilities import *
import numpy as np
#This command is for a prettier print of numpy arrays
np.set_printoptions(formatter={'float': '{: 0.3f}'.format})
from scipy.sparse import *
import time
import matplotlib.pyplot as plt
In [63]:
#Vector spaceya
vs=IteratedVectorSpace(UniformLagrangeVectorSpace(3), np.linspace(0,1,5))
vs.print_info()
In [64]:
cells=vs.n_cells
cells
Out[64]:
In [65]:
nbasis=vs.n_dofs
nbasis
Out[65]:
In [66]:
#test using sparce functions.
cum=dok_matrix((nbasis,nbasis))
(qst,wst)=np.polynomial.legendre.leggauss(5)
for lcell in range(vs.n_cells):
a,b=vs.cells[lcell],vs.cells[lcell+1]
q=((a+b)+(b-a)*qst)/2
w=(b-a)*wst/2
fnts=vs.cell_span(lcell)
nfnts=len(fnts)
fq=np.zeros((nfnts,5))
it=range(nfnts)
for i in it:
fq[i]=vs.basis(fnts[i])(q)
out=np.dot(fq*w,np.transpose(fq))
cont=-1
for rout, rcum in zip(out, fnts):
for vout, ccum in zip(rout, fnts):
cum[rcum,ccum]+=vout
#print using 3 formats
tmp=cum.todense()
print (type(tmp), "\n", tmp)
tmp=cum.toarray()
print (type(tmp), "\n", tmp)
print (type(cum), "\n", cum)
In [ ]:
#Using the default function that returns a numpy array no sparced.
t=MassMatrix(vs,5)
print (t)
In [68]:
#Using the CSC function
t2=MassMatrix(vs,5,format="CSC")
print (t2)
In [69]:
#Using the CSR option.
t2=MassMatrix(vs,5,format="CSR")
print (t2)
In [70]:
for i in range(2):
for j in range(2,i):
vs=IteratedVectorSpace(UniformLagrangeVectorSpace(j), np.linspace(0,1,i))
v=np.ones(vs.n_dofs)
for k in range(2,10):
tem=MassMatrix(vs,k)
assert abs(v.dot(tem.dot(v))-1.0)<10**(-10)
In [72]:
#test for to decide between using LIL or DOK
results=[]
print ("cels","nbasis","CSR","LIL")
for i in np.arange(2,5):
for j in np.arange(2,10):
vs=IteratedVectorSpace(UniformLagrangeVectorSpace(j), np.linspace(0,1,4**i))
v=np.ones(vs.n_dofs)
t1 = int(round(time.time() * 1000))
tdok=MassMatrix(vs,3)
t2 = int(round(time.time() * 1000))
tlil=MassMatrix(vs,3,internal="LIL")
t3 = int(round(time.time() * 1000))
lres=[4**i, j, t2-t1,t3-t2]
results.append(lres)
print (lres)
print ("done")