Hey, i am making small steps in the mixed element approach, however i can provide the iterative version, it works in terms of not resulting in erroers, however the results are not stable. It is not possible to reproduce and the chances are high, that the array, which contains the values is filled with nan entries.
The complete Code:
### General Declarations ###
from __future__ import print_function
from fenics import *
import matplotlib.pyplot as plt
import numpy as np
parameters['allow_extrapolation'] = True
### declare mesh and submesh ###
height = 1.0
width = 1.0
length = 1.0
mesh = BoxMesh(Point(0.0, 0.0, 0.0), Point(length,width,height), 9, 9, 9)
## Create submesh ##
marker = MeshFunction("size_t", mesh, mesh.topology().dim() - 1, 0)
for f in facets(mesh):
marker[f] = height - DOLFIN_EPS < f.midpoint().z() < height + DOLFIN_EPS
submesh = MeshView.create(marker, 1)
### Computation of Eigenvalues ###
X = FiniteElement('CG', submesh.ufl_cell(), 2)
Y = FiniteElement('CG', submesh.ufl_cell(), 2)
Z = X * Y
S = FunctionSpace(submesh, Z)
#Computation Eigenvalues
u = TrialFunctions(S)
phi = TestFunctions(S)
a = inner(grad(u[1]), grad(phi[0]))*dx + inner(grad(u[0]), grad(phi[1]))*dx
def boundary(x, on_boundary):
return on_boundary
bc = DirichletBC(S, Constant((0,0)), boundary)
L_dummy = Constant(0)*phi[0]*dx
m = u[0]*phi[0]*dx + u[1]*phi[1]*dx
A, _ = assemble_system(a,L_dummy,bc)
B = assemble(m)
eigensolver = SLEPcEigenSolver(as_backend_type(A), as_backend_type(B))
eigensolver.parameters['spectrum'] = 'target real'
eigensolver.parameters['tolerance'] = 1e-6
eigensolver.solve()
r, c, rx, cx = eigensolver.get_eigenpair(0)
print("Eigenvalue: ", r)
ef = Function(S)
ef.vector()[:] = rx
#%% ### The Coupled Problem ###
num_steps = 20 # number of time steps
dt = 0.005 # time step size
T = dt * num_steps # final time
t = 0
################ Functions/TrialFunctions/TestFunctions of the Plate ################
u_m_derivative_Plate = Function(S)
u_n_derivative_Plate = Function(S)
u_m_Plate = Function(S)
u_n_Plate = Function(S)
u_Plate = TrialFunctions(S)
phi_Plate = TestFunctions(S)
normalDerivative_Damping = Function(S)
################ FunctionSpace Function/TrialFunctions/TestFunctions damping ################
V_Damping = FunctionSpace(mesh, 'CG', 2)
v_Damping = TrialFunction(V_Damping)
psi_Damping = TestFunction(V_Damping)
boundaryFunction_Damping = Function(V_Damping)
Solution_Damping = Function(V_Damping)
#for plotting
valueArray_damping = np.array([])
t_plot = np.linspace(0, T, num = num_steps)
############## Boundaries Damping ################
def Top(x, on_boundary): return x[2] > height - DOLFIN_EPS
def Container(x, on_boundary):
return x[0] > length - DOLFIN_EPS or x[0] < DOLFIN_EPS or x[1] > width - DOLFIN_EPS or x[1] < DOLFIN_EPS
boundaryValue_Damping = Constant(0.0)
bcTop = DirichletBC(V_Damping, boundaryValue_Damping, Top)
bcContainer = DirichletBC(V_Damping, boundaryValue_Damping, Container)
bcs_Damping = [bcTop, bcContainer]
v_initial_Damping = interpolate(boundaryValue_Damping, V_Damping)
############## Boundary Plate ################
def boundary(x, on_boundary):
return on_boundary
bc_Plate = DirichletBC(S, Constant((0,0)), boundary)
############## bilinear/linear Form Plate and load ################
a_Plate = -dt*dt*(inner(grad(u_Plate[0]),grad(phi_Plate[0])) + inner(grad(u_Plate[1]),grad(phi_Plate[1])) + u_Plate[1]*phi_Plate[0])*dx + (u_Plate[0]*phi_Plate[1])*dx
f_Plate = Constant(0.0)
L_Plate = (2*u_n_Plate[0]*phi_Plate[1] - u_m_Plate[0]*phi_Plate[1])*dx + dt*dt*(inner(f_Plate*ef[0],phi_Plate[1]))*dx
solution_Plate=Function(S)
l_derivative_Plate = Function(S)
############## bilinear/linear Form Damping and load ################
f_Damping = Constant(0.0)
a_Damping = (inner(v_Damping,psi_Damping) + 2*dt * inner(grad(v_Damping), grad(psi_Damping)))*dx
L_Damping = (inner(v_initial_Damping,psi_Damping) + dt*inner(f_Damping,psi_Damping))*dx
############# Start Calculations ################
for n in range(num_steps):
#compute solution Damping
solve(a_Damping == L_Damping, Solution_Damping, bcs_Damping, solver_parameters={"linear_solver":"direct"})
v_initial_Damping.assign(Solution_Damping)
Solution_DampingZDerivative = Solution_Damping.dx(2)
plt.figure()
plot(Solution_DampingZDerivative)
normalDerivative_Damping = project(Solution_DampingZDerivative, S.sub(0).collapse(), solver_type="mumps")
# Compute solution Plate
solve(a_Plate == L_Plate, solution_Plate, bc_Plate, solver_parameters={"linear_solver":"direct"})
#save to plot
pointValue = solution_Plate(0.5,0.5,height)
valueArray_damping = np.append(valueArray_damping, pointValue[0])
u_m_Plate.assign(u_n_Plate)
u_m_derivative_Plate.assign(u_m_Plate)
#Update previous solution
u_n_Plate.assign(solution_Plate)
u_n_derivative_Plate.assign(solution_Plate)
l_derivative_Plate.assign((u_n_derivative_Plate - u_m_derivative_Plate)/2*dt)
speed, speedAppendix = l_derivative_Plate.split()
boundaryFunction_Damping = interpolate(speed, V_Damping)
bcTop = DirichletBC(V_Damping, boundaryFunction_Damping, Top)
bcContainer = DirichletBC(V_Damping, boundaryValue_Damping, Container)
bcs_Damping = [bcTop, bcContainer]
#Update current time
t += dt
f_Plate = Constant(2*r*sin(r*t))
L_Plate = (2*u_n_Plate[0]*phi_Plate[1] - u_m_Plate[0]*phi_Plate[1])*dx + dt*dt*(inner(f_Plate*ef[0],phi_Plate[1]) - inner(normalDerivative_Damping, phi_Plate[0]))*dx
#Update time of the expression
f_Plate.t = t
############# Plotting ################
plt.figure()
plt.plot(t_plot, valueArray_damping, 'b-')
plt.grid(axis='y')
plt.xlabel('t')
plt.ylabel('u',rotation=0)
plt.legend(['Entwicklung Auslenkung'], loc='upper left')
plt.savefig('GekoppeltesSystem.pdf')
