The answer is not as it should be. Please help!!
from dolfin import *
import fenics as fe
from fenics import *
import matplotlib.pyplot as plt
from ufl import nabla_grad
from ufl import nabla_div
import numpy as np
import pygmsh
import gmsh
import meshio
import pandas as pd
import os
import sympy as sp
import math
# Form compiler options
parameters["form_compiler"]["cpp_optimize"] = True
parameters["form_compiler"]["optimize"] = True
## Defining Boundaries
def bottom(x, on_boundary):
return (on_boundary and fe.near(x[1], 0.0, 1.0e-5))
def top(x, on_boundary):
return (on_boundary and fe.near(x[1], 1.0, 1.0e-5))
def left(x, on_boundary):
return (on_boundary and fe.near(x[0], 0.0, 1.0e-5))
def right(x, on_boundary):
return (on_boundary and fe.near(x[0], 1.0, 1.0e-5))
# Elastic parameters
E = 210.0
nu = 0.3
mu = Constant(E / (2.0*(1.0 + nu)))
lmbda = Constant(E*nu / ((1.0 + nu)*(1.0 - 2.0*nu)))
# Mass density
rho = Constant(1.0)
# Rayleigh damping coefficients
eta_m = Constant(0.)
eta_k = Constant(0.)
# Generalized-alpha method parameters
alpha_m = Constant(0.2)
alpha_f = Constant(0.4)
gamma = Constant(0.5+alpha_f-alpha_m)
beta = Constant((gamma+0.5)**2/4.)
# Time-stepping parameters
T = 1.0
Nsteps = 10
Rxlb = np.zeros(Nsteps+1)
Rxrb = np.zeros(Nsteps+1)
Rybb = np.zeros(Nsteps+1)
Rytb = np.zeros(Nsteps+1)
dt = Constant(T/Nsteps)
## Dimensions
L_x = 1.0
L_y = 1.0
## Meshing
gmsh.initialize()
with pygmsh.occ.Geometry() as geom:
geom.characteristic_length_min = 0.0
geom.characteristic_length_max = 0.03
rectangle = geom.add_rectangle([0, 0, 0], L_x, L_y)
mesh = geom.generate_mesh()
triangle_cells = None
for cell in mesh.cells:
if cell.type == "triangle":
# Collect the individual meshes
if triangle_cells is None:
triangle_cells = cell.data
else:
triangle_cells = np.vstack((triangle_cells, cell.data))
triangle_data = None
triangle_mesh = meshio.Mesh(points=mesh.points, cells={"triangle": triangle_cells})
nodes = mesh.points[:,0:2]
cells = triangle_cells
mesh = fe.Mesh()
editor = fe.MeshEditor()
editor.open(mesh,"triangle", 2, 2)
editor.init_vertices(len(nodes))
editor.init_cells(len(cells))
[editor.add_vertex(i,n) for i,n in enumerate(nodes)]
[editor.add_cell(i,n) for i,n in enumerate(cells)]
editor.close()
# Define function space for displacement, velocity and acceleration
V = VectorFunctionSpace(mesh, "Lagrange", 1)
# Current (unknown) displacement
u = Function(V, name="displacement")
# Fields from previous time step (displacement, velocity, acceleration)
u_old = Function(V)
v_old = Function(V)
a_old = Function(V)
# Test and trial functions
du = TrialFunction(V)
u_ = TestFunction(V)
I = Identity(len(u))
# Create mesh function over the cell facets
boundary_subdomains = fe.MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
boundary_subdomains.set_all(0)
AutoSubDomain(left).mark(boundary_subdomains, 1)
AutoSubDomain(right).mark(boundary_subdomains, 2)
AutoSubDomain(bottom).mark(boundary_subdomains, 3)
AutoSubDomain(top).mark(boundary_subdomains, 4)
# Define measure for boundary condition integral
dss = ds(subdomain_data=boundary_subdomains)
bc_l = DirichletBC(V.sub(0), 0.0, left)
left_bdry_dofs = np.fromiter(bc_l.get_boundary_values().keys(), dtype = int)
bc_r = DirichletBC(V.sub(0), 0.05, right)
right_bdry_dofs = np.fromiter(bc_r.get_boundary_values().keys(), dtype = int)
bc_b = DirichletBC(V.sub(1), 0.0, bottom)
bottom_bdry_dofs = np.fromiter(bc_b.get_boundary_values().keys(), dtype = int)
bc_t = DirichletBC(V.sub(1), 0.1, top)
top_bdry_dofs = np.fromiter(bc_t.get_boundary_values().keys(), dtype = int)
bc = [bc_l,bc_r,bc_b,bc_t]
def sigma(r):
return 2.0*mu*sym(grad(r)) + lmbda*tr(sym(grad(r)))*Identity(len(r))
# Mass form
def m(u, u_):
return rho*inner(u, u_)*dx
# Elastic stiffness form
def k(u, u_):
return inner(sigma(u), sym(grad(u_)))*dx
# Rayleigh damping form
def c(u, u_):
return eta_m*m(u, u_) + eta_k*k(u, u_)
# Work of external forces
def Wext(u_):
return dot(u_, p)*dss(3)
def P(dWdF,u_):
return inner(dWdF,grad(u_))*dx
def update_a(u, u_old, v_old, a_old, ufl=True):
if ufl:
dt_ = dt
beta_ = beta
else:
dt_ = float(dt)
beta_ = float(beta)
return (u-u_old-dt_*v_old)/beta_/dt_**2 - (1-2*beta_)/2/beta_*a_old
# Update formula for velocity
# v = dt * ((1-gamma)*a0 + gamma*a) + v0
def update_v(a, u_old, v_old, a_old, ufl=True):
if ufl:
dt_ = dt
gamma_ = gamma
else:
dt_ = float(dt)
gamma_ = float(gamma)
return v_old + dt_*((1-gamma_)*a_old + gamma_*a)
def update_fields(u, u_old, v_old, a_old):
"""Update fields at the end of each time step."""
# Get vectors (references)
u_vec, u0_vec = u.vector(), u_old.vector()
v0_vec, a0_vec = v_old.vector(), a_old.vector()
# use update functions using vector arguments
a_vec = update_a(u_vec, u0_vec, v0_vec, a0_vec, ufl=False)
v_vec = update_v(a_vec, u0_vec, v0_vec, a0_vec, ufl=False)
# Update (u_old <- u)
v_old.vector()[:], a_old.vector()[:] = v_vec, a_vec
u_old.vector()[:] = u.vector()
def avg(x_old, x_new, alpha):
return alpha*x_old + (1-alpha)*x_new
F = I + grad(u) # Deformation gradient
J = det(F)
C = F.T*F # Right Cauchy-Green tensor
B = F*F.T
# Invariants of deformation tensors
I1 = tr(C) + 1
I2 = 0.5*((tr(C)+1)**2 - (tr(C*C)+1))
I3 = J**2
I1_bar = J**(-2/3)*I1
I2_bar = J**(-4/3)*I2
F_inv = inv(I) - (1/(1+tr(grad(u)*inv(I))))*(inv(I)*grad(u)*inv(I))
# # Stored strain energy density
# W = 0.5*(I1_bar - 3) + 1.5*((J-1)**2)
# Calculating stresses
dWdF = (J*F_inv)*((((2*(0.5))/J**(5/3))*B+((2*0.0)/J**(5/3))*I1_bar*(B))-(2/(3*J))*(I1_bar*0.5+2*I2_bar*0.0)*I-((2/J**(7/3))*0.0*(B*B.T))+(3*J-3.0)*I)
# Residual
a_new = update_a(du, u_old, v_old, a_old, ufl=True)
v_new = update_v(a_new, u_old, v_old, a_old, ufl=True)
# res = m(avg(a_old, a_new, alpha_m), u_) + P(avg(u_old, du, alpha_f),I,u_)
res = m(avg(a_old, a_new, alpha_m), u_) + P(dWdF,u_)
# res = m(avg(a_old, a_new, alpha_m), u_) + P(dWdF,u_) + c(avg(v_old, v_new, alpha_f), u_) + k(avg(u_old, du, alpha_f), u_)
a_form = lhs(res)
L_form = rhs(res)
# Define solver for reusing factorization
K, res = assemble_system(a_form, L_form, bc)
solver = LUSolver(K, "mumps")
solver.parameters["symmetric"] = True
# Time-stepping
time = np.linspace(0, T, Nsteps+1)
for (i, dt) in enumerate(np.diff(time)):
t = time[i+1]
print("Time: ", t)
# Solve for new displacement
K = assemble(a_form)
res = assemble(L_form)
bc_l.apply(K,res)
bc_r.apply(K,res)
bc_t.apply(K,res)
bc_b.apply(K,res)
fe.solve(K, u.vector(), res)
# Update old fields with new quantities
update_fields(u, u_old, v_old, a_old)
f_ext_unknown = assemble(L_form)
## Sorting x and y degrees of freedom
x_dofs = V.sub(0).dofmap().dofs()
y_dofs = V.sub(1).dofmap().dofs()
## x component of residuals of nodes present in at left boundary
Rxlb_step = [0 for element in range(len(left_bdry_dofs))]
k = 0
for i in left_bdry_dofs:
Rxlb_step[k] = f_ext_unknown[i]
k +=1
## x component of residuals of nodes present in right boundary
Rxrb_step = [0 for element in range(len(right_bdry_dofs))]
k = 0
for i in right_bdry_dofs:
Rxrb_step[k] = f_ext_unknown[i]
k +=1
## y component of residuals of nodes present in at bottom boundary
Rybb_step = [0 for element in range(len(bottom_bdry_dofs))]
k = 0
for i in bottom_bdry_dofs:
Rybb_step[k] = f_ext_unknown[i]
k +=1
## y component of residuals of nodes present in top boundary
Rytb_step = [0 for element in range(len(top_bdry_dofs))]
k = 0
for i in top_bdry_dofs:
Rytb_step[k] = f_ext_unknown[i]
k +=1
Rxlb[i] = sum(Rxlb_step)
Rxrb[i] = sum(Rxrb_step)
Rybb[i] = sum(Rybb_step)
Rytb[i] = sum(Rytb_step)
reactions = [sum(Rxlb),sum(Rxrb),sum(Rybb),sum(Rytb)]
fe.plot(u, mode="displacement")
# fe.plot(a_old)
a_x_t = np.zeros(len(nodes))
a_y_t = np.zeros(len(nodes))
for n in range(0,len(nodes)):
a_x_t[n] = a_old(nodes[n,0],nodes[n,1])[0]
a_y_t[n] = a_old(nodes[n,0],nodes[n,1])[1]
print(reactions)