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
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm

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 [ ]:
from gmsh import MeshGmsh
from basis import LagrangeBasisHex
mnorm = lambda x: np.max(np.abs(x))

Setup mesh


In [ ]:
N = 4

periodic = False

In [ ]:
mesh = MeshGmsh()
#mesh.build("../dev/gmsh/cylinder.msh")
#mesh.build("../dev/gmsh/cube-cylinder-hole.msh")
mesh.build("../dev/gmsh/sphere.msh")


basis = LagrangeBasisHex(N)
mesh.add_basis(basis)

In [ ]:
# def f(X):
#     x = X[:,0]
#     y = X[:,1]
#     z = X[:,2]
    
#     return np.sin(np.pi*x)*np.sin(np.pi*y)*np.sin(np.pi*z)

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

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 [ ]:
lmap = LinearIsopMap()

etn = mesh.elem_to_node
Q, etd = mesh.Q, mesh.elem_to_dof
R = mesh.R
boundary_dofs = mesh.boundary_dofs

Poisson


In [ ]:
# Compatibility hacks
mesh.n_elem = mesh.n_elems

# Build the poisson problem
poisson = PoissonProblem(mesh, lmap)
poisson.build(mesh.nodes)

In [ ]:
p = poisson

nn = mesh.R.shape[0]
    
linOp = LinearOperator((nn, nn), matvec=p.apply_A)

Solve System


In [ ]:
p = poisson
dof_phys = p.dof_phys
fh  = f2(dof_phys)
fl = fh
rhs = p.B.dot(fl)
radj = np.zeros(p.n_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)

Solve with CG


In [ ]:
class CB(object):
    def __init__(self):
        self.n_iter = 0
    def __call__(self, x):
        self.n_iter += 1
        
cb = CB()
solcg, errc = sps.linalg.cg(linOp, rhs, tol=1e-14, 
                            maxiter=2000, callback=cb)

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
print norm(exact-solcg)/norm(exact)

In [ ]:
dp = p.dof_phys
err = (exact-solcg)[np.abs(dp[:,0])<1e-6]
dp = dp[np.abs(dp[:,0])<1e-6]

plt.scatter(dp[:,1],dp[:,2],
            c=err*100)


plt.show()
mnorm(err)

In [ ]:
dp = p.dof_phys
err = (exact-solcg)[np.abs(dp[:,2]-.5)<1e-6]
dp = dp[np.abs(dp[:,2]-.5)<1e-6]

plt.scatter(dp[:,0],dp[:,1],
            c=err*100)


plt.show()
mnorm(err)

In [ ]:
# fig = plt.figure()
# ax = fig.add_subplot(111, projection='3d')

# dp = p.dof_phys

# ax.scatter(dp[:,0], dp[:,1], dp[:,2])

# plt.show()

In [ ]:
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

dp = p.dof_phys[mesh.boundary_dofs]

ax.scatter(dp[:,0], dp[:,1], dp[:,2])

plt.show()

In [ ]: