Efficient creation of an extractor

dolfinx
1

I have MixedElement space V=U\times P and I would like to build a PETSc.Mat object that could be multiplied by a vector associated to a Function from the V space and return a PETSc.Vec associated to a Function living on the U space.

This is what I call an extractor, it’s essentially a potentially very large matrix filled with 0 & 1s.

My first naive attempt was to create a projection bilinear form, something like :

v=dfx.Function(V)
u=dfx.Function(V.sub(0).collapse(collapsed_dofs=True)[0])
form=ufl.inner(u,v)*ufl.dx
B=dfx.fem.assemble_matrix(form)
B.assemble()

This throws an error though : RuntimeError: Cannot create sparsity pattern. Form is not a bilinear form though if the matrix is not square, it does look bilinear to me.

Right now I have a much less elegant iterative solution that involves creating a CSR matrix :

U=V.sub(0)
U_collapsed=U.collapse()
# Compute DoFs
dofs=dfx.fem.locate_dofs_geometrical((U, U_collapsed), lambda x: np.ones(x.shape[1]))[0]
n=dofs.size
m=V.dofmap.index_map.size_global * V.dofmap.index_map_bs
rows=np.zeros(m+1,dtype='int32')
for i in dofs: rows[i+1:]+=1 # Probably inefficient
cols=np.arange(n,dtype='int32')
vals=np.ones(n,dtype='int32')
B=pet.Mat().createAIJ([m, n],csr=[rows,cols,vals]) # AIJ represents sparse matrix
B.setUp()
B.assemble()

No surprise there, the line with a for loop is punitively slow. Anyone has an idea on doing this in a smarter manner ?

2

Not exactly sure what you’re trying to achieve here, but maybe consider this example from tIGAr.

3

I’m afraid I’m running out of words to clarify this… I’ll try a more formal approach.

I have a couple system of equations (resolvant Navier Stokes) :

\begin{cases} u\nabla u=-\nabla p+\nu\Delta u+f\\ \nabla\cdot u=0 \end{cases}\Leftrightarrow A\left(\begin{array}{c}u\\p\end{array}\right)=\left(\begin{array}{c}I\\0\end{array}\right)f

Now writing the left part of that equation in dolfinx is easy, and I can plug that into a petsc4py solver to get the eigenvalues I want. I’m having difficulty with the extractor on the right ; basically I want to construct that matrix efficiently and if possible in parallel. I guess you could think of it as the matrix equivalent of the reverse of a ufl.split.

Thanks for the link, it looks nice and parallel, but I see you guys also went with a for loop on the dofs. I guess I could just take your approach and put 1 instead of nodesAndEvals.

4

Just a quick follow-up ; why does x_nodes = self.V.tabulate_dof_coordinates().reshape((-1,self.mesh.geometry().dim())) is included in the label loop ? That looks like a bug to me.