[petsc-users] "Exact-ish" preallocation for (mostly) unstructured mesh with periodic boundary conditions?

[Ramsey, James J CIV (US)]

I'm trying to get a better preallocation scheme. I noticed an earlier e-mail <https://lists.mcs.anl.gov/mailman/htdig/petsc-users/2010-September/006886.html> mentioning that "for finite elements with a normal cell division or finite differences with a normal vertex division, you need only count the nonzeros for rows you own since the partitions overlap on the boundary." Unfortunately, this doesn't quite work for me. Here's the latest attempt at a routine, with some details skipped:

std::vector<PetscInt> d_nnz(nLocalRows), o_nnz(nLocalRows);
std::vector<std::set<PetscInt> > d_nzInds(nLocalRows), o_nzInds(nLocalRows);

for ( ... loop over elements ...) {
    
   wrapNodeIds(*mesh, nIdVec);  // Accounts for periodic boundary conditions
   mesh->nodes2globalDOF(nIdVec, globalDOFsPerElem);

   std::vector<PetscInt>::iterator gdBegin = globalDOFsPerElem.begin();
   std::vector<PetscInt>::iterator gdEnd = globalDOFsPerElem.end();

   for (std::vector<PetscInt>::iterator itrI = gdBegin; itrI != gdEnd; ++itrI) {

       PetscInt I = *itrI;

       if ((I >= minDOFId) && (I <= maxDOFId)) {
          PetscInt Ioffset = I - minDOFId;
        
          for (std::vector<PetscInt>::iterator itrJ = gdBegin; itrJ != gdEnd; ++itrJ) {

            PetscInt J = *itrJ;
          
            if ((J >= minDOFId) && (J <= maxDOFId)) {
              // J will not be inserted if it is already in d_nzInds[Ioffset]. This avoids duplicates.
              d_nzInds[Ioffset].insert(J);
            }
            else {
              o_nzInds[Ioffset].insert(J);
            }

          }

     }
}

for (PetscInt i = 0; i < nLocalRows; ++i) {
      d_nnz[i] = d_nzInds[i].size();
      o_nnz[i] = o_nzInds[i].size();
}

For nodes in the interior, this seems to work well, but at the boundaries of the domain, the number of non-zeros is underestimated. It seems like parallel communication is unavoidable here, but I'm not sure how to do it, and I'm unsure if I really want to mess with relatively low-level MPI routines to get what I want.

Any ideas?