I have defined a function Continuum_relaxation (code attached) which solved a linear elasticity problem with some periodic and Dirichlet BCs. However the solver machinery is more general Non-linear one. It has Periodic BCs on the x and y directions. One of the z-boundary is fixed to zero displacement. Whereas on the other z-boundary every node has a Dirichlet BC.
When I run this code, I get the following error:
File "AC_coupling.py", line 456, in Continuum_relaxation
solver.solve()
RuntimeError:
*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
*** fenics@fenicsproject.org
*** Remember to include the error message listed below and, if possible,
*** include a minimal running example to reproduce the error.
*** -------------------------------------------------------------------------
*** Error: Unable to access vector of degrees of freedom.
*** Reason: Cannot access a non-const vector from a subfunction.
*** Where: This error was encountered inside Function.cpp.
*** Process: 0
*** DOLFIN version: 1.6.0
*** Git changeset: unknown
I am learning fenics and have no idea how to solve this issue. I will very grateful if anyone can help me out.
Thank you very much.
Hope to hear from some of you soon.
Sincerely,
Shankha
def Continuum_relaxation(interface_atom_disp, n_pad_atoms, nodes, orient, xlo, xhi, ylo, yhi, z_end_continuum, mu, lmbda, alpha):
""" Create mesh and define function space """
# Importing the mesh
mesh = Mesh("mesh.xml")
mesh.init(2) # Compute faces
mesh.init(0, 0) # Compute vertex neighbors for each vertex
mesh.init(1, 1) # Compute edge neighbors for each edge
""" Defining Sub Domains for Boundary Conditions"""
"""Sub domain for Periodic boundary condition"""
class PeriodicDomain(SubDomain):
#Defining the "target domain" G
def inside(self, x, on_boundary):
return bool((near(x[0], xhi) or near(x[1], yhi)) and (not ((near(x[0], xlo) and near(x[1], yhi)) or (near(x[0], xhi) and near(x[1], ylo)))) and on_boundary)
#Mapping domain H to "target domain" G
def map(self, x, y):
if near(x[0], xlo) and near(x[1], ylo):
y[0] = x[0] + (xhi - xlo)
y[1] = x[1] + (yhi - ylo)
y[2] = x[2]
elif near(x[0], xlo):
y[0] = x[0] + (xhi - xlo)
y[1] = x[1]
y[2] = x[2]
elif near(x[1], ylo):
y[0] = x[0]
y[1] = x[1] + (yhi - ylo)
y[2] = x[2]
else:
y[0] = -1000
y[1] = -1000
y[2] = -1000
"""Sub domains for every interface node"""
interface = [CompiledSubDomain("x[0] == xi && x[1] == yi && x[2] == zi", xi = pos[0], yi = pos[1], zi = pos[2]) for pos in nodes[np.array([list(nodes[:,0]).index(i) for i in interface_atom_disp[:,0]]), 1:]]
"""Sub domain for free surface boundary opposite to the interface"""
def free_surface(x, on_boundary):
return bool(near(x[2], z_end_continuum) and (not ((near(x[2], z_end_continuum) and near(x[0], xhi)) or (near(x[2], z_end_continuum) and near(x[1], yhi)))) and on_boundary)
""" Function space with linear functions(deg=1) """
V = VectorFunctionSpace(mesh, "Lagrange", 1, constrained_domain=PeriodicDomain())
""" Assigning Dirichlet Boundary Conditions """
u0_interface = [Constant(disp) for disp in interface_atom_disp[:,1:]]
u0_free_surface = Constant((0.0, 0.0, 0.0))
bc_interface = [DirichletBC(V, u0_interface[count] , interface[count]) for count in range(0, np.shape(interface_atom_disp)[0])]
bc_free_surface = [DirichletBC(V, u0_free_surface, free_surface)]
bcs = bc_interface + bc_free_surface
""" Define variational problem """
du = TrialFunction(V) # Incremental displacement
v = TestFunction(V) # Test function
u = Function(V) # Displacement from previous iteration
# Define material model
def material_model(u):
# Kinematics
epsilon = sym(grad(u))
Id = Identity(u.geometric_dimension())
nu1 = np.array([np.dot(orient['x'], np.array([1, 0, 0])) / np.linalg.norm(orient['x']), np.dot(orient['y'], np.array([1, 0, 0])) / np.linalg.norm(orient['y']), np.dot(orient['z'], np.array([1, 0, 0])) / np.linalg.norm(orient['z'])])
nu2 = np.array([np.dot(orient['x'], np.array([0, 1, 0])) / np.linalg.norm(orient['x']), np.dot(orient['y'], np.array([0, 1, 0])) / np.linalg.norm(orient['y']), np.dot(orient['z'], np.array([0, 1, 0])) / np.linalg.norm(orient['z'])])
nu3 = np.array([np.dot(orient['x'], np.array([0, 0, 1])) / np.linalg.norm(orient['x']), np.dot(orient['y'], np.array([0, 0, 1])) / np.linalg.norm(orient['y']), np.dot(orient['z'], np.array([0, 0, 1])) / np.linalg.norm(orient['z'])])
# Material law
i,j = indices(2)
O1 = as_tensor(Constant(nu1)[i] * Constant(nu1)[j], (i,j))
O2 = as_tensor(Constant(nu2)[i] * Constant(nu2)[j], (i,j))
O3 = as_tensor(Constant(nu3)[i] * Constant(nu3)[j], (i,j))
sigma = (2.0*mu*epsilon) + (lmbda * tr(epsilon) * Id) + (alpha*((inner(epsilon,O1)*O1) + (inner(epsilon,O2)*O2) + (inner(epsilon,O3)*O3)))
jacobian = 0
return (sigma,jacobian)
# Bilinear form definitions
Sigma , J = material_model(u)
L = inner(Sigma, nabla_grad(v))*dx
traction = Constant((0.0, 0.0, 0.0)) # Traction force on the boundary
F = dot(traction, v)*ds
B = L - F
J = derivative(B, u, du)
""" Compute solution """
# Set up the problem and solver
problem = NonlinearVariationalProblem(B, u, bcs, J)
solver = NonlinearVariationalSolver(problem)
# Set up solver parameters
prm = solver.parameters
prm['newton_solver']['absolute_tolerance'] = 1E-8
prm['newton_solver']['relative_tolerance'] = 1E-7
prm['newton_solver']['maximum_iterations'] = 10
prm['newton_solver']['relaxation_parameter'] = 1.0
# Solve the problem
solver.solve()
""" Output """
disp = np.vstack(np.split(u.vector().array()[vertex_to_dof_map()], mesh.num_vertices()))
pad_atom_disp = np.column_stack((np.array(range(0, n_pad_atoms)), disp[0:n_pad_atoms,:]))
id_coor_disp = np.column_stack((np.array(range(0, mesh.num_vertices())), mesh.coordinates(), disp))
return pad_atom_disp, id_coor_disp'
