On Tue, 28 Jul 2015 07:33:18 +0200 Julian Andrej <j...@tf.uni-kiel.de> wrote:
> Hello, > > in a minimal code example i get the error [0] Matrix is missing > diagonal entry in row 25 if i use the PETSc method zeroRows() to apply > my boundary conditions. I use this formulation because i split my > formulation for the operators of the Stokes operators on a mixed > function in order to decouple the pressure. > The documentation states, that keep_diagonal preserves all diagonal > entries, which is maybe not the case here? It is the case. The matrix has shape (162, 25) and all diagonal entries are indeed there. Just try M.mat().view() > > from dolfin import * > from petsc4py import PETSc > > mesh = UnitSquareMesh(4, 4) > V = VectorFunctionSpace(mesh, "CG", 2) > Q = FunctionSpace(mesh, "CG", 1) > > u = TrialFunction(V) > v = TestFunction(V) > p = TrialFunction(Q) > q = TestFunction(Q) > > a = div(v) * p * dx > > # Define boundary condition > u0 = Constant((0.0, 0.0)) > bc = DirichletBC(V, u0, 'on_boundary') > > M = PETScMatrix() > assemble(a, tensor=M, keep_diagonal=True) > > # This operation succeeds > bc.zero(M) > > # Extract the boundary dofs > bcdict = bc.get_boundary_values() > bcinds = bcdict.keys() > > # This one fails > M.mat().zeroRows(bcinds) The actual problem is that PETSc refuses to change allocated zeros on diagonal into nonzeros because of MAT_KEEP_NONZERO_PATTERN flag we set in DOLFIN. You can do M.mat().setOption(PETSc.Mat.Option.KEEP_NONZERO_PATTERN, False) to workaround this. But are you sure that this is what you want to do? Shouldn't it be rather M.mat().zeroRows([i for i in bcinds if i < M.size(1)]) M.mat().zeroRows([i for i in bcinds if i >= M.size(1)], diag=0.0) Jan > _______________________________________________ > fenics mailing list > fenics@fenicsproject.org > http://fenicsproject.org/mailman/listinfo/fenics _______________________________________________ fenics mailing list fenics@fenicsproject.org http://fenicsproject.org/mailman/listinfo/fenics
