# PETsc error code is : 63

I am trying to solve three pde’s with mixedElement. I suspect i did somethin wrong in applying boundary condition but i am new to Fenics so would really appreciate any kind of help.

from dolfin import *
from mshr import *
import numpy as np
import matplotlib.pyplot as plt
#%matplotlib notebook

L = 1.
R = 0.1
N = 50  # mesh density

domain = Rectangle(Point(0., 0.), Point(L, L)) - Circle(Point(0., 0.), R, 100)
mesh = generate_mesh(domain, N)

T0 = Constant(293.)
DThole = Constant(10.)
E = 70e3
nu = 0.3
lmbda = Constant(E*nu/((1+nu)*(1-2*nu)))
mu = Constant(E/2/(1+nu))
rho = Constant(2700.)     # density
alpha = 2.31e-5  # thermal expansion coefficient
kappa = Constant(alpha*(2*mu + 3*lmbda))
cV = Constant(910e-6)*rho  # specific heat per unit volume at constant strain
k = Constant(237e-6)  # thermal conductivity
T = 2.0            # final time
num_steps = 50     # number of time steps
dt = T / num_steps # time step size
A1 = 40
A2 = 0.5
E1 = 30000
E2 = 60000
m = 0.7
n = 1.2
Rg = 8.31
T0 = Constant(283.)

Vue = VectorElement('CG', mesh.ufl_cell(), 2) # displacement finite element
Vte = FiniteElement('CG', mesh.ufl_cell(), 1) # temperature finite element
Vce = FiniteElement('CG', mesh.ufl_cell(), 1) # Curing finite element
VT = FunctionSpace(mesh, MixedElement([Vue, Vte, Vce]))

def inner_boundary(x, on_boundary):
return near(x[0]**2+x[1]**2, R**2, 1e-3) and on_boundary
def bottom(x, on_boundary):
return near(x[1], 0) and on_boundary
def left(x, on_boundary):
return near(x[0], 0) and on_boundary

bc1 = DirichletBC(VT.sub(0).sub(1), Constant(0.), inner_boundary)
bc2 = DirichletBC(VT.sub(0).sub(1), Constant(0.), bottom)
bc3 = DirichletBC(VT.sub(0).sub(1), Constant(0.), left)

bc4 = DirichletBC(VT.sub(1), Constant(1000.), inner_boundary)
bc5 = DirichletBC(VT.sub(1), Constant(400.), bottom)
bc6 = DirichletBC(VT.sub(1), Constant(400.), left)
bc7 = DirichletBC(VT.sub(2), Constant(1.), inner_boundary)
bc8 = DirichletBC(VT.sub(2), Constant(0.1), bottom)
bc9 = DirichletBC(VT.sub(2), Constant(0.1), left)
#bc4 = DirichletBC(VT.sub(1), Constant(100.), boundary)
bcT = [bc1, bc2, bc3, bc4, bc5, bc6, bc7, bc8, bc9]

u = Function(VT)
du, dTheta, dalpha = split(u)
v = TestFunction(VT)
u_, Theta_, Alpha_ = split(v)
Uold = Function(VT)
uold, Theta_old, Alpha_old = split(Uold)
f = Constant(5)
Theta = dTheta - T0
def eps(v):

def sigma(v, Theta):
return (lmbda*tr(eps(v)) - kappa*Theta)*Identity(2) + 2*mu*eps(v)

def K_1(dTheta):
return A1*exp(-E1/(Rg*dTheta))
def K_2(dTheta):
return A2*exp(-E2/(Rg*dTheta))
def S_c(dalpha, dTheta):
return (K_1(dTheta) + (K_2(dTheta)*pow(dalpha, m)))*pow((1-dalpha),n)

dt = Constant(0.)

mech_form = inner(sigma(du, dTheta), eps(u_))*dx
#cure_form = ((dalpha - Alpha_old)/dt)*Alpha_*dx - S_c(Alpha_old)*Alpha_*dx
cure_form = dalpha*Alpha_*dx - Alpha_old*Alpha_*dx - dt*S_c(Alpha_old, Theta_old)*Alpha_*dx
form = mech_form + therm_form + cure_form

#Pushp, Pushp_to_VT

Theta_o = Expression('293.', degree=1)
Alpha_o = Expression('0.0293', degree=2)
#uold, Theta_old, Alpha_old = Uold.split()

Theta_old_tmp = VT.sub(1).collapse()
Alpha_old_tmp = VT.sub(2).collapse()

Theta_old = interpolate(Theta_o, Theta_old_tmp)
Alpha_old = interpolate(Alpha_o, Alpha_old_tmp)
parameters['allow_extrapolation'] = True # TODO
Nincr = 100
t = np.logspace(1, 4, Nincr+1)
#Nx = 100
#x = np.linspace(-5, 10, Nx)
#T_res = np.zeros((Nx, Nincr+1))
C = Function(VT)
#u = Function(VT)
for (i, dti) in enumerate(np.diff(t)):
print("Increment " + str(i+1))
dt.assign(dti)
solve(form == 0, C, bcT)
print ("vdfv")

Uold.assign(u)
#T_res[:, i+1] = [C(xi, 0.)[1] for xi in x]
#vtkfile << (T_res, t)


I am getting the following error:

Error:   Unable to successfully call PETSc function 'MatSetValuesLocal'.
*** Reason:  PETSc error code is: 63 (Argument out of range)"


Please edit the post so that the code is formatted properly, using “” instead of the blockquote. Furthermore, the code must be complete, i.e. include all required imports.

Remove the definition of C here. The input to

should be the function you used to define your residual, i.e.

It worked. Thanks. but i was trying to follow the logic and procedure given in the following link which solves coupled thermo-mechanical eqns.

U = Function(V)`