In [ ]:

    

%matplotlib inline
%load_ext autoreload
%autoreload 2
%load_ext line_profiler


    

In [ ]:

    

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


    

In [ ]:

    

from pyamg.classical import ruge_stuben_solver
norm  = lambda x: np.max(np.abs(x)) if len(x)>0 else 0.0
kron3 = lambda x,y,z: sps.kron(x,sps.kron(y,z))


    

In [ ]:

    

from tensormesh import HexCubePoisson
from maps import LinearIsopMap
from topology import CubicTopology
from poisson import PoissonProblem


    

In [ ]:

    

class CB(object):
    def __init__(self):
        self.n_iter = 0
    def __call__(self, x):
        self.n_iter += 1


    

In [ ]:

    

def solve_problem(f, f2, E, N, periodic):
    
    Ex, Ey, Ez = E
    
    lmap = LinearIsopMap()

    topo = CubicTopology(N, (Ex, Ey, Ez),
                        periodic=periodic)
    topo.build()

    etn = topo.elem_to_vertex
    Q, etd = topo.Q, topo.elem_to_dof
    R = topo.R
    boundary_dofs = topo.boundary_dofs

    vertex_ref = topo.get_vertex_ref()

    vertex_phys = vertex_ref.copy()

    chi, eta, zeta = vertex_ref.T
    sx = sy = sz = 0.1
    vp = vertex_phys
    vp[:,0] = chi +sx*np.sin(np.pi*chi)*np.sin(np.pi*eta)*np.sin(np.pi*zeta)
    vp[:,1] = eta +sy*np.sin(np.pi*chi)*np.sin(np.pi*eta)*np.sin(np.pi*zeta)
    vp[:,2] = zeta+sz*np.sin(np.pi*chi)*np.sin(np.pi*eta)*np.sin(np.pi*zeta)

    poisson = PoissonProblem(topo, lmap)
    poisson.build(vertex_phys)

    p = poisson
    if periodic:
        nn = p.n_dofs
    else:
        nn = (p.nz_dofs-2)*(p.ny_dofs-2)*(p.nx_dofs-2)

    linOp = LinearOperator((nn, nn), matvec=p.apply_A)

    M = HexCubePoisson(N,Ex,L=2,periodic=periodic)
    M.build_mesh()
    precond = LinearOperator((nn,nn), 
                             matvec=M.solve)

    dof_phys = p.dof_phys
    fh  = f2(dof_phys)
    fl = fh
    rhs = p.B.dot(fl)
    radj = np.zeros(p.nx_dofs*p.ny_dofs*p.nz_dofs)
    radj[boundary_dofs] = f(dof_phys)[boundary_dofs]
    rhs = R.dot(rhs-p.apply_A(radj, apply_R=False))
    exact = f(dof_phys)

    if periodic:
        rhs -= np.mean(rhs)    


    cb = CB()
    solcg, errc = sps.linalg.cg(linOp, rhs, tol=1e-14, 
                                maxiter=2000, callback=cb,
                                M=precond)

    solcg = R.T.dot(solcg)
    if periodic:
        solcg -= np.mean(solcg)
        exact -= np.mean(exact)
    else:
        solcg[boundary_dofs] = f(dof_phys[boundary_dofs])

    print cb.n_iter, norm(rhs-p.apply_A(R.dot(solcg))),
    print norm(exact-solcg)/norm(exact)
    
    return norm(exact-solcg)/norm(exact)


    

In [ ]:

    

def f(X):
    x = X[:,0]
    y = X[:,1]
    z = X[:,2]
    
    p = np.pi*2
    r = np.cos(p*x)*np.cos(p*y)*np.cos(p*z)
    #if not periodic:
    r += x
    return r

def f2(X):
    x = X[:,0]
    y = X[:,1]
    z = X[:,2]
    
    p = np.pi*2
    return np.cos(p*x)*np.cos(p*y)*np.cos(p*z)*3*(p)**2


    

In [ ]:

    

Na = [1,2,3]
Ea = [4,8,16]

erra = []
for N in Na:
    err = []
    for E in Ea:
        err += [solve_problem(f, f2, (E,E,E), N, False)]
    erra += [err]


    

In [ ]: