Here is complete working example. Thanks to Miro who provided the core code.
from mshr import *
from dolfin import *
# Test for PETSc
if not has_linear_algebra_backend("PETSc"):
print("DOLFIN has not been configured with PETSc. Exiting.")
exit()
##############################
### Geometry, Domain, Mesh ###
##############################
class FreeBoundary(SubDomain):
def inside(self, x, on_boundary):
return on_boundary
class LeftBoundary(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and abs(x[0] + 1) < DOLFIN_EPS
class RightBoundary(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and abs(x[0] - 1) < DOLFIN_EPS
# Create geo and mesh
res = 100
a = 1 # oter square
b = 0.5 # inner sqare - hole a plate
r1 = Rectangle(Point(-a, -a), Point(a, a))
r2 = Rectangle(Point(-b, -b), Point(b, b))
geo = r1 - r2 # Square plate with square hole
mesh = generate_mesh(geo,res)
###########################################
### 2D Linear Elasticity Material Model ###
###########################################
# Elasticity parameters
E, nu = 10.0, 0.3 # Young's module and Poisson ratio
mu_, lambda_ = Constant(E/(2*(1 + nu))), Constant(E*nu/((1 + nu)*(1 - 2*nu)))# Lame coefficients
# Strain
epsilon = lambda u: sym(grad(u))
# Stress
sigma = lambda u: 2*mu_*epsilon(u) + lambda_*tr(epsilon(u))*Identity(2)
# Nullspace of 2D rigid motions
# two translations + one rotation
nullspace = [Constant((1, 0)), Constant((0, 1)), Expression(('-x[1]', 'x[0]'), degree=1)]
#
########################################################################
####### Preliminary check for zero eigenvalues ###########
########################################################################
#
# Let's first show that the energy form without constraints is not pos.def
V = VectorFunctionSpace(mesh, 'CG', 1) # Space for displacement
u, v = TrialFunction(V), TestFunction(V)
# Energy
a = inner(sigma(u), epsilon(v))*dx
A = PETScMatrix()
assemble(a, A)
# Every z is eigenv with eigenvalue 0
print '1. Eigenvalues:'
for z in nullspace:
v = interpolate(z, V).vector()
null = v.copy()
print null.norm('l2'),
A.mult(v, null)
print null.norm('l2')
print
########################################################################
# Confirm that A has three zero eigenvalue with slepc
eigensolver = SLEPcEigenSolver(A)
eigensolver.parameters['problem_type'] = 'hermitian'
eigensolver.parameters['spectrum'] = 'smallest magnitude'
eigensolver.solve(5)
assert eigensolver.get_number_converged() > 3
print '2. Eigenvalues:'
for i in range(4):
w = eigensolver.get_eigenvalue(i)
print w
print
########################################################################
# Get rid of singularity
M = VectorFunctionSpace(mesh, 'R', 0, 3) # Space for all multipliers
W = MixedFunctionSpace([V, M])
u, mus = TrialFunctions(W)
v, nus = TestFunctions(W)
# Energy
a = inner(sigma(u), epsilon(v))*dx
# Lagrange multipliers contrib to a
for i, e in enumerate(nullspace):
mu = mus[i]
nu = nus[i]
a += mu*inner(v, e)*dx + nu*inner(u, e)*dx
# See that there are now not zero eigenvalues
A = PETScMatrix()
assemble(a, A)
eigensolver = SLEPcEigenSolver(A)
eigensolver.parameters['problem_type'] = 'hermitian'
eigensolver.parameters['spectrum'] = 'smallest magnitude'
eigensolver.solve(5)
assert eigensolver.get_number_converged() > 1
print '3. Eigenvalues:'
for i in range(1):
w = eigensolver.get_eigenvalue(i)
print w
########################################################################
#
########################################################################
##### Solve the problem #####
########################################################################
#
# Energy
a = inner(sigma(u), epsilon(v))*dx
# Lagrange multipliers contrib to a
for i, e in enumerate(nullspace):
mu = mus[i]
nu = nus[i]
a += mu*inner(v, e)*dx + nu*inner(u, e)*dx
#Initialize subdomain instances
free_b = FreeBoundary()
left_b = LeftBoundary()
right_b = RightBoundary()
# Initialize mesh function for boundary domains
boundaries = FacetFunction("size_t", mesh)
boundaries.set_all(0)
free_b.mark(boundaries, 1)
left_b.mark(boundaries, 2)
right_b.mark(boundaries, 3)
# Boundary conditions - Neumann (stress on L&R side) - Traction
gL = Constant((-1.0, 0.0))
gR = Constant((1.0, 0.0))
# Define new measures associated with the exterior boundaries
ds = Measure('ds', domain=mesh, subdomain_data=boundaries)
# Natural BC
L = dot(gL, v)*ds(2) + dot(gR, v)*ds(3)
A, b = assemble_system(a, L)
print "Norm(b)=",b.norm("l2") # Check that b is not zero
print
print
# Solve
w_h = Function(W)
# Solve with direct solver
problem = LinearVariationalProblem(a, L, w_h)
solver = LinearVariationalSolver(problem)
solver.parameters["linear_solver"] = "direct"
solver.parameters["symmetric"] = True
solver.solve()
# Solution split
u_h, mus_h = w_h.split(deepcopy=True)
u_x, u_y = u_h.split()
# Post processing displacement
plot(u_h, mode = "displacement", interactive=True, scalarbar=True, title="u")
File("u_x.pvd", "compressed") << u_x
File("u_y.pvd", "compressed") << u_y
# Post processing stress
Fs = FunctionSpace(mesh, 'CG', 1)
Ts = TensorFunctionSpace(mesh, 'DG', 0, symmetry=True)
ss = sigma(u_h)
ssT = project(ss, Ts)
sxx = project(ssT[0,0], Fs)
syy = project(ssT[1,1], Fs)
sxy = project(ssT[0,1], Fs)
tau_max = project(sqrt( ( (sxx - syy)/2 )**2 + sxy**2 ), Fs) #Tau_max
# Plot similar to photoelastic experimetnt
plot(tau_max, interactive=True, rescale=False, scalarbar=True, title="Tau_max")
File("sxx.pvd", "compressed") << sxx
File("syy.pvd", "compressed") << syy
File("sxy.pvd", "compressed") << sxy
File("tau_max.pvd", "compressed") << tau_max
