I am trying to plot the solution to my PDE problem and when I turn it into a numpy array using solution_array.vector().get_local() and plot it again I see a different solution being plotted.
Can anyone point out the fault ?
from fenics import *
import numpy as np
from dolfin import *
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
# Step 1: Define the Problem (Geometry and Material Parameters)
L_k = 1e-6 # Length of the domain in the y direction
lam = 1e-6
K1 = Constant(3.81e-12) # Splay elastic constant
K2 = Constant(3.81e-10) # Twist elastic constant
q0 = Constant(0)
# Create mesh and define function space
dim_x = 16
dim_y = 16
# Sub domain for Periodic boundary condition
class PeriodicBoundary(SubDomain):
# Left boundary is "target domain" G
def inside(self, x, on_boundary):
return bool(x[0] < DOLFIN_EPS and on_boundary)
# Map right boundary (H) to left boundary (G)
def map(self, x, y):
y[0] = x[0]-lam
y[1] = x[1]
# Create mesh and finite element
mesh = RectangleMesh(Point(0,0), Point(lam, L_k),dim_x-1,dim_y-1) # Define the rectangular mesh
#1d Experiment
#mesh_1D = IntervalMesh(dim_x,0,L_k)
V = FunctionSpace(mesh, "CG", 1)#constrained_domain=PeriodicBoundary())
# Create Dirichlet boundary condition at phi(x,0)
theta0 = Constant(0) #Expression("DOLFIN_PI", degree=1,lam=lam)
def bottom_boundary(x, on_boundary):
return near(x[0], 0.0) and on_boundary
bc_bottom = DirichletBC(V,theta0,bottom_boundary)
# Boundary condition at y=L (takes only odd multiples of pi/2)
u_k = Constant(1)
def top_boundary(x, on_boundary):
return near(x[0],L_k) and on_boundary
theta_L = Constant(DOLFIN_PI/2)#Expression("(2*u_k+1)*DOLFIN_PI/2",degree=1,u_k = u_k)
bc_top = DirichletBC(V, theta_L,top_boundary)
# Collect boundary conditions
bcs = [bc_bottom,bc_top]
# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
f = Constant(0)
a = K1 * dot(grad(u)[0], grad(v)[0]) * dx + K2 * dot(grad(u)[1], grad(v)[1]) * dx
#a = dot(grad(u), grad(v)) * dx
L = f*v*dx
# Compute solution
phi = Function(V)
solve(a == L, phi, bcs)
# Plot solutionT
#plot(mesh,scalarbar = False,title="Finite Element Mesh")
plot(phi,scalarbar = False,interactive=True)
# Convert the FEniCS Function to a numpy array
xyz = V.tabulate_dof_coordinates()
#print(xyz[:,0])
(xyz)
n_array = phi.vector().get_local()
n_array= n_array.reshape((dim_x,dim_y))
print("Here is n_array")
# Create the 3D plot
fig = plt.figure()
ax = plt.axes(projection='3d')
ax.view_init(30,60)
# Data for a three-dimensional line
Z = np.fliplr(n_array.T)
x_line = np.linspace(0,lam,dim_x)
y_line = np.linspace(0,L_k,dim_y)
X,Y= np.meshgrid(y_line,x_line)
ax.plot_surface(X, Y, np.cos(Z), rstride=1, cstride=1,
cmap='viridis', edgecolor='none')
ax.contour3D(X,Y,Z, 50, cmap='binary')
ax.set_title('surface');
#Line-plot
plt.figure()
#plt.xlabel("x")
#plt.ylabel("Director $\phi(x,z)$ ")
plt.plot(x_line,Z[:,dim_x//2],color = 'r') # Bottom phi = pi*x/LAMBDA
plt.plot(x_line,Z[:,dim_x-1],color = 'g') # top phi = (2*u+1)/2*pi
plt.figure()
plt.plot(y_line,Z[dim_x//2,:],color = 'r') # Bottom phi = pi*x/LAMBDA
#plt.plot(y_line,Z[dim_x-1,:],color = 'g') # top phi = (2*u+1)/2*pi
#plt.imshow(np.cos(Z))
#2d Image render
plt.figure()
#quiveropts = dict(pivot='middle')
plt.imshow(Z)
plt.figure()
plt.quiver(Y,X,np.cos(Z),np.zeros_like(Z),pivot='mid')