Hi Everyone,
I am trying to fetch one value (K1_ij) from an assembled stiffness matrix and assign the value to the other assembled stiffness matrix element (K0_ij). I am using functions GenericMatrix.get and GenericMatrix.set. The GenericMatrix.get works well and I can print out the K1_ij, which is extracted from the assembled stiffness K1 by K1.get. But after I assign K1_ij to the other stiffness matrix K0 with K0.set, something wrong happens to K0.
Attached below is a minimal working exmaple. It models a 2D in-plane elastic deformation. The top boundary of the model is fixed and the bottom boundary moves 0.1 to the right. The model size is 10 by 10 and it is discretized to 10 cells along each direction. K1 is twice of K0 and I am trying to replace the (5,5) element in K0 with the (5,5) element of K1. Before calling K0.set, I can solve the system with K0. But after the modification, I get an error as shown below when trying to solve the new system (other operations on K0 will also fail, such as applying bc to it or exporting K0 to a file). I am not sure what is going wrong. If there is other way to do this job, I would also like to try. The dolfin version is 2019.1.0 and I installed it via Anaconda.
Many thanks for your kind help!
Dunyu
RuntimeError Traceback (most recent call last)
in
64 K0.set(stiff_comp, i1, j1)
65
—> 66 dl.solve(K0, u_n.vector(), b ,‘lu’)
67 print(‘Solve after modification of K0’, u_n.vector())
~/software/anaconda3/envs/fenicsproject/lib/python3.9/site-packages/dolfin/fem/solving.py in solve(*args, **kwargs)
225 raise RuntimeError(“Not expecting keyword arguments when solving linear algebra problem.”)
226
→ 227 return dolfin.la.solver.solve(*args)
228
229
~/software/anaconda3/envs/fenicsproject/lib/python3.9/site-packages/dolfin/la/solver.py in solve(A, x, b, method, preconditioner)
70 “”"
71
—> 72 return cpp.la.solve(A, x, b, method, preconditioner)
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-support@googlegroups.com
*** Remember to include the error message listed below and, if possible,
*** include a minimal running example to reproduce the error.
*** -------------------------------------------------------------------------
*** Error: Unable to successfully call PETSc function ‘KSPSolve’.
*** Reason: PETSc error code is: 73 (Object is in wrong state).
*** Where: This error was encountered inside /home/conda/feedstock_root/build_artifacts/fenics-pkgs_1617882212586/work/dolfin/dolfin/la/PETScKrylovSolver.cpp.
*** Process: 0
*** DOLFIN version: 2019.1.0
*** Git changeset: 7f46aeb0b296da5bbb1fb0845822a72ab9b09c55
*** -------------------------------------------------------------------------
A MWE is attached below to show the problem.
#!/usr/bin/env python
import dolfin as dl
import numpy as np
mesh1 = dl.RectangleMesh(dl.Point(0.0, 0.0), dl.Point(10., 10.), 10, 10)
boundaries = dl.MeshFunction("size_t", mesh1, 1)
class Top_boundary(dl.SubDomain):
def inside(self, x, on_boundary):
return dl.near(x[1], 10., dl.DOLFIN_EPS) and on_boundary
class Bottom_boundary(dl.SubDomain):
def inside(self, x, on_boundary):
return dl.near(x[1], 0.0, dl.DOLFIN_EPS) and on_boundary
boundaries.set_all(0)
Top_boundary().mark(boundaries, 1)
Bottom_boundary().mark(boundaries, 2)
nu = 0.25
mu = 1
lam = mu*2.*nu/(1.-2.*nu)
# Define functions for displacement formulation
def epsilon(u):
eps = dl.sym(dl.nabla_grad(u))
return eps
def stress_E(u):
return 2.0*mu*epsilon(u) + lam*dl.tr(epsilon(u))*dl.Identity(2)
V = dl.VectorFunctionSpace(mesh1, 'CG', 1)
u_ = dl.TrialFunction(V)
w_ = dl.TestFunction(V)
u_n = dl.Function(V)
f = dl.Constant((0.0,0.0))
a0 = dl.inner(stress_E(u_), epsilon(w_))*dl.dx
K0= dl.assemble(a0)
K1= dl.assemble(2.*a0)
L = dl.inner(f, w_)*dl.dx
b = dl.assemble(L)
bc = [dl.DirichletBC(V.sub(0), dl.Constant(0.0), boundaries, 1),
dl.DirichletBC(V.sub(1), dl.Constant(0.0), boundaries, 1),
dl.DirichletBC(V.sub(0), dl.Constant(0.1), boundaries, 2)]
# Solve the system before the modification of K0
for ibc in bc:
ibc.apply(K0)
ibc.apply(b)
dl.solve(K0, u_n.vector(), b ,'lu')
print('Solve before modification of K0')
print('The solution is ', u_n.vector()[:])
i1 = [5]
j1 = [5]
stiff_comp = np.zeros((1,1))
K1.get(stiff_comp, i1, j1)
print('Get one component in K1', i1, j1, stiff_comp)
# Solve the system after the modification of K0
K0.set(stiff_comp, i1, j1)
dl.solve(K0, u_n.vector(), b ,'lu')
print('Solve after modification of K0', u_n.vector())