I have a free energy density equation given by the formula:
f = \frac{K}{2}[(\frac{\partial{\phi}}{\partial{x}})^2 + (\frac{\partial{\phi}}{\partial{z}} - q_0)^2]
where K and q0 are constants and \phi(x,z) is a scalar field which outputs an angle defined on each point in the domain of x and z.
The total energy given by
F =\int_{0}^{L}\int_{0}^{d} f \,dx\,dz
L and d are the dimensions of the rectangular domain of where the function f is defined.
Now I want to find \phi(x,z) that satisfy the solution to the equation \delta F = 0 .
There are two Dirichlet boundary conditions given below:
\phi(x,0) =\frac{\pi x}{\Lambda} lambda is constant
and
\phi(x,L) =(2u+1) \frac{\pi}{2}
where u is an integer.
My first attempt was the following
from fenics import *
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
pi =3.141
# Step 1: Define the Problem (Geometry and Material Parameters)
L = 14e-6 # Length of the domain in the y direction
d = 14e-6 # Length of the domain in the x direction
lam = 60e-6
K = Constant(7.49e-12) # Splay elastic constant
K2 = Constant(3.81e-12) # Twist elastic constant
q0 = Constant(0)
# Step 2: Define Function Space
mesh = RectangleMesh(Point(0, 0), Point(d, L), 15, 15) # Define the rectangular mesh
V = FunctionSpace(mesh, "CG", 1) # Vector function space for the director field
# Step 3: Define Trial and Test Functions
n = Function(V) # The director field
#v = TestFunction(V)
# Step 4: Define the Free Energy Functional
F = 0.5 * K * inner(grad(n)[0], grad(n)[0]) * dx +0.5*K2*inner(grad(n)[1]-q0, grad(n)[1]-q0) * dx
# Step 5: Compute the Derivative of the Free Energy Functional (Variational Problem)
dFdn = derivative(F, n)
# Step 6: Apply Boundary Conditions
# Boundary condition at y=0 (linear with respect to x)
bx0 = Expression("pi*x[0]/lam", degree=1,lam=lam,pi=pi)
bc_x0 = DirichletBC(V, bx0, "near(x[1], 0)")
#Step 7 : Apply Boundary Conditions
# Boundary condition at y=L (takes only odd multiples of pi/2)
u = Constant(0)
theta_L = Expression("pi/2*(2*u+1)", pi=pi,u=u, degree=1)
bc_L = DirichletBC(V, theta_L, "near(x[1], 14e-6)")
# Combine the boundary conditions
bc = [bc_x0, bc_L]
# Step 7: Compute the Equilibrium Director Field
#n_eq = TestFunction(V) # The equilibrium director field
#var_form = inner(F,n_eq) * dx
solve(dFdn == 0, n, bc)
# Step 8: Post-process and visualize the results (e.g., plot the director field)
# Plot solution and mesh
plot(n)
#plot(mesh)
# Save solution to file in VTK format
vtkfile = File('poisson/solution.pvd')
vtkfile << n
# Convert the FEniCS Function to a numpy array
n_array = n.vector().get_local()
n_array= n_array.reshape((16,16))
# Create the 3D plot
fig = plt.figure()
ax = plt.axes(projection='3d')
ax.view_init(60,30)
# Data for a three-dimensional line
Z = n_array
x_line = np.linspace(0,d,16)
y_line = np.linspace(0,L,16)
X,Y = np.meshgrid(x_line,y_line)
ax.plot_surface(X, Y, 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.plot(Z[:,0])
#2d Image render
plt.figure()
plt.quiver(X,Y,np.cos(Z),0)
I am pretty sure teh solutions are wrong. Any help in solving this problem would be appreciated