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))
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
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)
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)
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 [ ]: