How to apply Dirichlet BCs to a matrix generated with petsc4py

Question

I am trying to understand why boundary conditions cannot be applied to a matrix originally built using petsc4py, even after wrapping it with dolfin's PETScMatrix wrapper. This code shows how to create the error.

import dolfin as dlf

from petsc4py import PETSc

mesh = dlf.UnitSquareMesh(2,2)
W = dlf.FunctionSpace(mesh, 'CG', 1)

def boundary(x, on_boundary):
    return on_boundary

bc = dlf.DirichletBC(W, dlf.Constant(0.0), boundary)

N = mesh.num_vertices()
v = PETSc.Vec()
v.create()
v.setSizes(N)
v.setType('standard')
v.setValues(range(N), [N]*N)

B_pet = PETSc.Mat()
B_pet.createAIJ([N,N], nnz=N)
B_pet.setDiagonal(v)
B_pet.assemblyBegin()
B_pet.assemblyEnd()

B = dlf.PETScMatrix(B_pet)
bc.apply(B) # error

The error given by this code is:

*** Error:   Unable to set given (local) rows to identity matrix.                                                                                             
*** Reason:  some diagonal elements not preallocated (try assembler option keep_diagonal).                                                                    
*** Where:   This error was encountered inside PETScMatrix.cpp.

Answer 1

Hi, consider

import dolfin as dlf

from petsc4py import PETSc

mesh = dlf.UnitSquareMesh(100, 100)
W = dlf.FunctionSpace(mesh, 'CG', 1)

def boundary(x, on_boundary):
    return on_boundary

bc = dlf.DirichletBC(W, dlf.Constant(0.0), boundary)

N = mesh.num_vertices()
v = PETSc.Vec()
v.create()
v.setSizes(N)
v.setType('standard')
v.setValues(range(N), [N]*N)

B_pet = PETSc.Mat()
B_pet.createAIJ([N,N], nnz=N)

lgmap = PETSc.LGMap().create(W.dofmap().dofs())
B_pet.setLGMap(lgmap, lgmap)

B_pet.setDiagonal(v)
B_pet.assemblyBegin()
B_pet.assemblyEnd()

B = dlf.PETScMatrix(B_pet)
bc.apply(B) # error 

Comments

That works, thank you! I do have one more question though, how would I get this to work in parallel? I'm guessing the dofmap will somehow have to be determined for the local process.. but I'm not sure how.

---

Hi, could you ask this as a separate question? Something like: Generating DOLFIN matrices with PETSC. I have the answer, but I would not like for it to get lost in this thread, which on the outside appears too specifc to boundary conditions.

---

Thank you! I have submitted the additional question here.