[petsc-users] Including constrained dofs in a matrix/vector created with DMCreateXXXX()

[Jordan Wagner via petsc-users]

Hi,

Over the past few months, I have implemented dmplex/section in my 
preexisting finite element code. We already had our own implementations 
for things like the finite element space, boundary conditions, 
projections, etc, so I am mostly using dmplex/section at its lowest 
level; no DT/DS stuff for the moment. Everything is working fine at this 
level, and I am overall very happy with how it ties into the code that 
we already have.

I am currently doing my own preallocation for the system Jacobian but am 
wanting to use DMCreateMatrix() instead, as it seems pointless to keep 
using my own crappy function when this one is just sitting there 
waiting. However, I dislike how I am currently implementing this and was 
hoping to get some suggestions or some advice. My problem is handling 
nicely the Dirichlet dofs. My code takes the approach of first 
assembling the entire Jacobian matrix and load vector as if no 
constraints are imposed. We then have a function that applies the 
essential conditions after assembly and extracts a submatrix and 
subvector similar to the way ex49 
<https://www.mcs.anl.gov/petsc/petsc-current/src/ksp/ksp/examples/tutorials/ex49.c.html> 
is doing.

Since DMCreateMatrix() automatically partitions out the constraints 
specified in the section, I find myself having to create two nearly 
identical sections: one that has constraints and one that does not, and 
then I clone the DM and set the default section of each respectively. I 
can then use DMCreateMatrix() on the one dm/section with no constraints 
to preallocate the larger matrix that I am assembling. Then I use the 
dm/section with constraints to do all of my boundary condition and 
projection operations on the previously created matrix and vector. 
Essentially, I am using an unconstrained section only for allocation and 
a constrained section for most everything else.

This process actually works fine for me, but its pretty ugly, and I am 
sure there is a better way. I am wondering if I am missing something 
that could keep me from having to go through this process of creating 
two sections that differ only in the constraints. It seems to me if 
there were an option for DMCreateMatrix()/DMCreateVector() to either 
include or not include the dofs, that would completely solve my problem. 
That way, I can use the same section for both creation and allocation.

So, the main question I have is just whether there is a function or 
option that I am missing like DMCreateXXXX() but with the ability to 
retain the constrained rows and columns rather than compressing them out.

Thanks a lot, as always.