How to print the tangent stiffness matrix for a given problem?
jiachenguoNU opened this issue · 0 comments
Hi:
I'm new to fenics and I want to check the tangent stiffness matrix for a given problem. Is there anyway I can obtain this? I'm running the following simple linear elasticity code.
from fenics import *
from ufl import nabla_div
import numpy as np
Scaled variables
L = 1; W = 1
E = 100.
nu = 0.3
mu = E/(2.(1. + nu))
lambda_ = Enu/((1+nu)(1-2nu))
rho = 1
delta = W/L
gamma = 0.4*delta**2
beta = 1.25
g = 10.
Create mesh and define function space
mesh = BoxMesh(Point(0, 0, 0), Point(L, W, W), 1, 1, 1)
V = VectorFunctionSpace(mesh, 'P', 1)
Define boundary condition
tol = 1E-14
def clamped_boundary(x, on_boundary):
return on_boundary and x[0] < tol
def right_boundary(x, on_boundary):
return on_boundary and near(x[0], 1., tol)
bc = DirichletBC(V, Constant((0, 0, 0)), clamped_boundary)
bc2 = DirichletBC(V.sub(1), Constant((0.01)), right_boundary)
Define strain and stress
def epsilon(u):
return 0.5*(nabla_grad(u) + nabla_grad(u).T)
#return sym(nabla_grad(u))
def sigma(u):
return lambda_nabla_div(u)Identity(d) + 2muepsilon(u)
Define variational problem
u = TrialFunction(V)
d = u.geometric_dimension() # space dimension
v = TestFunction(V)
f = Constant((0, 0, 0))
T = Constant((0, 0, 0))
a = inner(sigma(u), epsilon(v))*dx
L = dot(f, v)*dx + dot(T, v)*ds
Compute solution
u = Function(V)
A,b = assemble_system(a,L,[bc, bc2])
solve(a == L, u, [bc, bc2])
Plot solution
plot(u, title='Displacement', mode='displacement')
Plot stress
s = sigma(u) - (1./3)*tr(sigma(u))Identity(d) # deviatoric stress
von_Mises = sqrt(3./2inner(s, s))
V = FunctionSpace(mesh, 'P', 1)
von_Mises = project(von_Mises, V)
plot(von_Mises, title='Stress intensity')
Compute magnitude of displacement
u_magnitude = sqrt(dot(u, u))
u_magnitude = project(u_magnitude, V)
plot(u_magnitude, 'Displacement magnitude')
print('min/max u:',
u_magnitude.vector().get_local().min(),
u_magnitude.vector().get_local().max())