<div dir="ltr"><div dir="ltr">On Tue, Feb 19, 2019 at 1:13 AM Jordan Wagner via petsc-users <<a href="mailto:petsc-users@mcs.anl.gov">petsc-users@mcs.anl.gov</a>> wrote:<br></div><div class="gmail_quote"><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">
  
    
  
  <div bgcolor="#FFFFFF">
    <p>Hi, <br>
    </p>
    <p>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. <br>
    </p>
    <p>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 <a href="https://www.mcs.anl.gov/petsc/petsc-current/src/ksp/ksp/examples/tutorials/ex49.c.html" target="_blank">ex49</a>
      is doing. <br>
    </p>
    <p>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.<br>
    </p>
    <p>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.</p>
    <p>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.</p>
    <p></p></div></blockquote><div>I do not completely understand the above, but here is how Section works. A _local_ Section contains all dofs and is intended for assembly using local Vecs. To communicate with the solver, you make a _global_ Section using PetscSectionCreateGlobalSection(). This has an option to include or exclude constraints. For example, I make a global Section with constraints when I output for visualization. Maybe you can use that to do what you want.</div><div><br></div><div>What I do not understand is why you do not just use a global Section without constraints to assemble your matrix, rather than eliminating the rows/columns later. The global Section gives back negative indices for unknown that are constrained or not owned. These are ignored by MatSetValues() and VecSetValues(). This is how the default Plex assembly makes operators without constrained variables.</div><div><br></div><div>  Thanks,</div><div><br></div><div>     Matt</div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div bgcolor="#FFFFFF"><p>Thanks a lot, as always.<br>
    </p>
  </div>

</blockquote></div><br clear="all"><div><br></div>-- <br><div dir="ltr" class="gmail_signature"><div dir="ltr"><div><div dir="ltr"><div><div dir="ltr"><div>What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br>-- Norbert Wiener</div><div><br></div><div><a href="http://www.cse.buffalo.edu/~knepley/" target="_blank">https://www.cse.buffalo.edu/~knepley/</a><br></div></div></div></div></div></div></div></div>