set rows of rectangular petsc matrix to zero

Question

Hello Fenics users
I will get straight to the point:
Let's assume I have a rectangular petsc Matrix B constructed as follows

V : Finite Element space of velocity (P2)
Q : Finite Element space of pressure (P1)
v = TestFunction(V)
p = TestFunction(Q)
B = assemble(nabla_grad(p)*v*dx, tensor = PETScMatrix())

and

bc = DirichletBC(V, 0, boundary)  

The question is :
is there a way to set rows of the Matrix B corresponding to Dirichlet bcs to zero ?

thanx in advance

Answer 1

To get straight to the answer, not asking why you actually would like to set those matrix rows to zero:
you can manually access and change values in GenericMatrix() or GenericVector() instances, see this post. So yes, you can set a row in the GenericMatrix() object to zero as follows:

from dolfin import * 
import numpy as np
mesh = UnitSquareMesh(3,3)
V = FunctionSpace(mesh, 'CG', 1)
p,q = TrialFunction(V), TestFunction(V)

B = assemble(dot(grad(p),grad(q))*dx, tensor = PETScMatrix())

# Set the first row (idx_row) to zero
idx_row = np.array([0],dtype=np.intc)
B.zero(idx_row)
B.apply('insert')
print B.array()

Comments

Actually I want a link between bc defined as

bc = DirichletBC(V, 0, boundary)  

and the functions you provided. In such a way that I can us the bc object to vanish to Dirichlet rows of B.
because I don't know boundary nodes indices (the mesh I use is imported from an xml file).

---

So if I understand you correctly, you want to know the dof indices at the DirichletBC? Then have a look here, where the SubsetIterator is suggested to extract dof indices in a SubDomain.

---

Thank you this helped me a little bit, and I think I'm close to the answer
I followed the link you posted an tried to get boundary dofs
so I wrote something like that on my program

# Unique vertices of the line (by indices)
vs = []
for bd_id in Dirichlet_BCs_physical_ids():
    vs += list(set(sum((facets.entities(0).tolist() for facets in SubsetIterator(boundaries, bd_id)), [])))

# Get degrees of freedom associated to vertices (by indices)
v2d = V.dofmap().dofs()
d = v2d[vs]

print d

unfortunatly, the output does not contain all the dofs !!! knowing that I m using a Finite Element space of order 2, I m getting only dofs of 1st order !
any idea ?

---

My first thought was to use tabulate_facet_dofs, which should look something like:

fctlist = list(f.index() for f in SubsetIterator(facet_f, 1))
# Then print dofs associated to facets
for fct in fctlist: 
    print V.dofmap().tabulate_facet_dofs(fct)

This returns however not the correct result for me.

Therefore, I propose the following (and working!) workaround in order to get the dof indices on the boundary:

from dolfin import *
import numpy as np

# Define boundary function
def domain_boundary(x, on_boundary):
    return on_boundary

mesh = UnitSquareMesh(2, 2)
Q    = FunctionSpace(mesh, 'CG', 2)

# Define Dirichlet BC
bc = DirichletBC(Q, 1, domain_boundary)

# Make dummy function from Q
q = Function(Q)
bc.apply(q.vector())

# Printing the dof indices on boundary
print np.where(q.vector() == 1)

Hopefully this solves your problem?

---

Indeed yes, This with some modifications just solved the problem.

Thank you for the help.

Bests

Answer 2

As simple as

bc.zero(B)