text:
   # assume communication to 6 neighbors by default
using MPI
function  get_neighbors(domain::Domain)
   rowMin = !(domain.m_rowLoc == 0)
   rowMax = !(domain.m_rowLoc == domain.m_tp - 1)
   colMin = !(domain.m_colLoc == 0)
   colMax = !(domain.m_colLoc == domain.m_tp - 1)
   planeMin = !(domain.m_planeLoc == 0)
   planeMax = !(domain.m_planeLoc == domain.m_tp - 1)

   return rowMin, rowMax, colMin, colMax, planeMin, planeMax
end

copyto_zero!(dest, doffs, src, soffs, nelems) = copyto!(dest, doffs+1, src, soffs+1, nelems)

# Assume 128 byte coherence
# Assume Float64 is an "integral power of 2" bytes wide
const CACHE_COHERENCE_PAD_REAL = div(128, sizeof(Float64))

function commRecv(domain::Domain, msgType, xferFields, dx, dy, dz, doRecv, planeOnly)
   comm = domain.comm
   if comm === nothing
       return
   end
   comm = comm::MPI.Comm

   # post recieve buffers for all incoming messages
   maxPlaneComm = xferFields * domain.maxPlaneSize
   maxEdgeComm  = xferFields * domain.maxEdgeSize
   pmsg = 0 # plane comm msg
   emsg = 0 # edge comm msg
   cmsg = 0 # corner comm msg

   baseType = MPI.Datatype(Float64) # TODO support Float32

   # assume communication to 6 neighbors by default
   rowMin, rowMax, colMin, colMax, planeMin, planeMax = get_neighbors(domain)

   fill!(domain.recvRequest, MPI.Request())

   myRank = MPI.Comm_rank(comm)

   # post receives
   function irecv!(fromProc, offset, recvCount)
      idx = offset + 1
      data = view(domain.commDataRecv, idx:(idx+recvCount))
      return MPI.Irecv!(MPI.Buffer(data), fromProc, msgType, domain.comm::MPI.Comm)
   end

   # receive data from neighboring domain faces
   if planeMin && doRecv
      # contiguous memory
      fromProc = myRank - domain.m_tp^2
      recvCount = dx * dy * xferFields
      req = irecv!(fromProc, pmsg*maxPlaneComm, recvCount)
      domain.recvRequest[pmsg+1] = req
      pmsg += 1
   end

   if planeMax
      # contiguous memory
      fromProc = myRank + domain.m_tp^2
      recvCount = dx * dy * xferFields
      req = irecv!(fromProc, pmsg*maxPlaneComm, recvCount)
      domain.recvRequest[pmsg+1] = req
      pmsg += 1
   end

   if rowMin && doRecv
      # semi-contiguous memory
      fromProc = myRank - domain.m_tp
      recvCount = dx * dz * xferFields
      req = irecv!(fromProc, pmsg*maxPlaneComm, recvCount)
      domain.recvRequest[pmsg+1] = req
      pmsg += 1
   end

   if rowMax
      # semi-contiguous memory
      fromProc = myRank + domain.m_tp
      recvCount = dx * dz * xferFields
      req = irecv!(fromProc, pmsg*maxPlaneComm, recvCount)
      domain.recvRequest[pmsg+1] = req
      pmsg += 1
   end

   if colMin && doRecv
      # scattered memory
      fromProc = myRank - 1
      recvCount = dy * dz * xferFields
      req = irecv!(fromProc, pmsg*maxPlaneComm, recvCount)
      domain.recvRequest[pmsg+1] = req
      pmsg += 1
   end

   if colMax
      # scattered memory
      fromProc = myRank + 1
      recvCount = dy * dz * xferFields
      req = irecv!(fromProc, pmsg*maxPlaneComm, recvCount)
      domain.recvRequest[pmsg+1] = req
      pmsg += 1
   end

   if !planeOnly
      # receive data from domains connected only by an edge
      if rowMin && colMin && doRecv
         fromProc = myRank - domain.m_tp - 1
         recvCount = dz * xferFields
         offset = pmsg * maxPlaneComm + emsg * maxEdgeComm
         req = irecv!(fromProc, offset, recvCount)
         domain.recvRequest[pmsg+emsg+1] = req
         emsg += 1
      end

      if rowMin && planeMin && doRecv
         fromProc = myRank - domain.m_tp^2 - domain.m_tp
         recvCount = dx * xferFields
         offset = pmsg * maxPlaneComm + emsg * maxEdgeComm
         req = irecv!(fromProc, offset, recvCount)
         domain.recvRequest[pmsg+emsg+1] = req
         emsg += 1
      end

      if colMin && planeMin && doRecv
         fromProc = myRank - domain.m_tp^2 - 1
         recvCount = dy * xferFields
         offset = pmsg * maxPlaneComm + emsg * maxEdgeComm
         req = irecv!(fromProc, offset, recvCount)
         domain.recvRequest[pmsg+emsg+1] = req
         emsg += 1
      end

      if rowMax && colMax
         fromProc = myRank + domain.m_tp + 1
         recvCount = dz * xferFields
         offset = pmsg * maxPlaneComm + emsg * maxEdgeComm
         req = irecv!(fromProc, offset, recvCount)
         domain.recvRequest[pmsg+emsg+1] = req
         emsg += 1
      end

      if rowMax && planeMax
         fromProc = myRank + domain.m_tp^2 + domain.m_tp
         recvCount = dx * xferFields
         offset = pmsg * maxPlaneComm + emsg * maxEdgeComm
         req = irecv!(fromProc, offset, recvCount)
         domain.recvRequest[pmsg+emsg+1] = req
         emsg += 1
      end

      if colMax && planeMax
         fromProc = myRank + domain.m_tp^2 + 1
         recvCount = dy * xferFields
         offset = pmsg * maxPlaneComm + emsg * maxEdgeComm
         req = irecv!(fromProc, offset, recvCount)
         domain.recvRequest[pmsg+emsg+1] = req
         emsg += 1
      end

      if rowMax && colMin
         fromProc = myRank + domain.m_tp - 1
         recvCount = dz * xferFields
         offset = pmsg * maxPlaneComm + emsg * maxEdgeComm
         req = irecv!(fromProc, offset, recvCount)
         domain.recvRequest[pmsg+emsg+1] = req
         emsg += 1
      end

      if rowMin && planeMax
         fromProc = myRank + domain.m_tp^2 - domain.m_tp
         recvCount = dx * xferFields
         offset = pmsg * maxPlaneComm + emsg * maxEdgeComm
         req = irecv!(fromProc, offset, recvCount)
         domain.recvRequest[pmsg+emsg+1] = req
         emsg += 1
      end

      if colMin && planeMax
         fromProc = myRank + domain.m_tp^2 - 1
         recvCount = dy * xferFields
         offset = pmsg * maxPlaneComm + emsg * maxEdgeComm
         req = irecv!(fromProc, offset, recvCount)
         domain.recvRequest[pmsg+emsg+1] = req
         emsg += 1
      end

      if rowMin && colMax && doRecv
         fromProc = myRank - domain.m_tp + 1
         recvCount = dz * xferFields
         offset = pmsg * maxPlaneComm + emsg * maxEdgeComm
         req = irecv!(fromProc, offset, recvCount)
         domain.recvRequest[pmsg+emsg+1] = req
         emsg += 1
      end

      if rowMax && planeMin && doRecv
         fromProc = myRank - domain.m_tp^2 + domain.m_tp
         recvCount = dx * xferFields
         offset = pmsg * maxPlaneComm + emsg * maxEdgeComm
         req = irecv!(fromProc, offset, recvCount)
         domain.recvRequest[pmsg+emsg+1] = req
         emsg += 1
      end

      if colMax && planeMin && doRecv
         fromProc = myRank - domain.m_tp^2 + 1
         recvCount = dy * xferFields
         offset = pmsg * maxPlaneComm + emsg * maxEdgeComm
         req = irecv!(fromProc, offset, recvCount)
         domain.recvRequest[pmsg+emsg+1] = req
         emsg += 1
      end

      # receive data from domains connected only by a corner
      if rowMin && colMin && planeMin && doRecv
         # corner at domain logical coord (0, 0, 0)
         fromProc = myRank - domain.m_tp^2 - domain.m_tp - 1
         recvCount = xferFields
         offset = pmsg * maxPlaneComm + emsg * maxEdgeComm + cmsg * CACHE_COHERENCE_PAD_REAL
         req = irecv!(fromProc, offset, recvCount)
         domain.recvRequest[pmsg+emsg+cmsg+1] = req
         cmsg += 1
      end

      if rowMin && colMin && planeMax
         # corner at domain logical coord (0, 0, 1)
         fromProc = myRank + domain.m_tp^2 - domain.m_tp - 1
         recvCount = xferFields
         offset = pmsg * maxPlaneComm + emsg * maxEdgeComm + cmsg * CACHE_COHERENCE_PAD_REAL
         req = irecv!(fromProc, offset, recvCount)
         domain.recvRequest[pmsg+emsg+cmsg+1] = req
         cmsg += 1
      end

      if rowMin && colMax && planeMin && doRecv
         # corner at domain logical coord (1, 0, 0)
         fromProc = myRank - domain.m_tp^2 - domain.m_tp + 1
         recvCount = xferFields
         offset = pmsg * maxPlaneComm + emsg * maxEdgeComm + cmsg * CACHE_COHERENCE_PAD_REAL
         req = irecv!(fromProc, offset, recvCount)
         domain.recvRequest[pmsg+emsg+cmsg+1] = req
         cmsg += 1
      end

      if rowMin && colMax && planeMax
         # corner at domain logical coord (1, 0, 1)
         fromProc = myRank + domain.m_tp^2 - domain.m_tp + 1
         recvCount = xferFields
         offset = pmsg * maxPlaneComm + emsg * maxEdgeComm + cmsg * CACHE_COHERENCE_PAD_REAL
         req = irecv!(fromProc, offset, recvCount)
         domain.recvRequest[pmsg+emsg+cmsg+1] = req
         cmsg += 1
      end

      if rowMax && colMin && planeMin && doRecv
         # corner at domain logical coord (0, 1, 0)
         fromProc = myRank - domain.m_tp^2 + domain.m_tp - 1
         recvCount = xferFields
         offset = pmsg * maxPlaneComm + emsg * maxEdgeComm + cmsg * CACHE_COHERENCE_PAD_REAL
         req = irecv!(fromProc, offset, recvCount)
         domain.recvRequest[pmsg+emsg+cmsg+1] = req
         cmsg += 1
      end

      if rowMax && colMin && planeMax
         # corner at domain logical coord (0, 1, 1)
         fromProc = myRank + domain.m_tp^2 + domain.m_tp - 1
         recvCount = xferFields
         offset = pmsg * maxPlaneComm + emsg * maxEdgeComm + cmsg * CACHE_COHERENCE_PAD_REAL
         req = irecv!(fromProc, offset, recvCount)
         domain.recvRequest[pmsg+emsg+cmsg+1] = req
         cmsg += 1
      end

      if rowMax && colMax && planeMin && doRecv
         # corner at domain logical coord (1, 1, 0)
         fromProc = myRank - domain.m_tp^2 + domain.m_tp + 1
         recvCount = xferFields
         offset = pmsg * maxPlaneComm + emsg * maxEdgeComm + cmsg * CACHE_COHERENCE_PAD_REAL
         req = irecv!(fromProc, offset, recvCount)
         domain.recvRequest[pmsg+emsg+cmsg+1] = req
         cmsg += 1
      end

      if rowMax && colMax && planeMax
         # corner at domain logical coord (1, 1, 1)
         fromProc = myRank + domain.m_tp^2 + domain.m_tp + 1
         recvCount = xferFields
         offset = pmsg * maxPlaneComm + emsg * maxEdgeComm + cmsg * CACHE_COHERENCE_PAD_REAL
         req = irecv!(fromProc, offset, recvCount)
         domain.recvRequest[pmsg+emsg+cmsg+1] = req
         cmsg += 1
      end
   end
end

function commSend(domain::Domain, msgType, fields,
                  dx, dy, dz, doSend, planeOnly)

   comm = domain.comm
   if comm === nothing
      return
   end

   xferFields = length(fields)

   # post recieve buffers for all incoming messages
   maxPlaneComm = xferFields * domain.maxPlaneSize
   maxEdgeComm  = xferFields * domain.maxEdgeSize
   pmsg = 0 # plane comm msg
   emsg = 0 # edge comm msg
   cmsg = 0 # corner comm msg

   # MPI_Status status[26] ;

   # assume communication to 6 neighbors by default
   rowMin,rowMax, colMin, colMax, planeMin, planeMax = get_neighbors(domain)

   fill!(domain.sendRequest, MPI.Request())
   myRank = MPI.Comm_rank(comm)

   # post sends
   if planeMin || planeMax
      # ASSUMING ONE DOMAIN PER RANK, CONSTANT BLOCK SIZE HERE
      sendCount = dx * dy

      if planeMin
         # contiguous memory
         srcOffset = 0
         offset = pmsg * maxPlaneComm
         for field in fields
            copyto_zero!(domain.commDataSend, offset, field, srcOffset, sendCount)
            offset += sendCount
         end
         idx = pmsg * maxPlaneComm + 1
         src = view(domain.commDataSend, idx:(idx+(xferFields * sendCount)))

         otherRank = myRank - domain.m_tp^2
         req = MPI.Isend(MPI.Buffer(src), otherRank, msgType, comm)
         domain.sendRequest[pmsg+1] = req
         pmsg += 1
      end

      if planeMax && doSend
         # contiguous memory
         srcOffset = dx*dy*(dz - 1)
         offset = pmsg * maxPlaneComm
         for field in fields
            copyto_zero!(domain.commDataSend, offset, field, srcOffset, sendCount)
            offset += sendCount
         end
         idx = pmsg * maxPlaneComm + 1
         src = view(domain.commDataSend, idx:(idx+(xferFields * sendCount)))

         otherRank = myRank + domain.m_tp^2
         req = MPI.Isend(MPI.Buffer(src), otherRank, msgType, comm)
         domain.sendRequest[pmsg+1] = req
         pmsg += 1
      end
   end

   if rowMin | rowMax
      # ASSUMING ONE DOMAIN PER RANK, CONSTANT BLOCK SIZE HERE
      sendCount = dx * dz

      if rowMin
         # contiguous memory
         srcOffset = 0
         offset = pmsg * maxPlaneComm
         for field in fields
            for i in 0:(dz-1)
               copyto_zero!(domain.commDataSend, offset+i*dx, field, srcOffset+i*dx*dy, dx)
            end
            offset += sendCount
         end
         idx = pmsg * maxPlaneComm + 1
         src = view(domain.commDataSend, idx:(idx+(xferFields * sendCount)))

         otherRank = myRank - domain.m_tp
         req = MPI.Isend(MPI.Buffer(src), otherRank, msgType, comm)
         domain.sendRequest[pmsg+1] = req
         pmsg += 1
      end
      if rowMax && doSend
         # contiguous memory
         srcOffset = dx*(dy - 1)
         offset = pmsg * maxPlaneComm
         for field in fields
            for i in 0:(dz-1)
               copyto_zero!(domain.commDataSend, offset+i*dx , field, srcOffset+i*dx*dy, dx)
            end
            offset += sendCount
         end
         idx = pmsg * maxPlaneComm + 1
         src = view(domain.commDataSend, idx:(idx+(xferFields * sendCount)))

         otherRank = myRank + domain.m_tp
         req = MPI.Isend(MPI.Buffer(src), otherRank, msgType, comm)
         domain.sendRequest[pmsg+1] = req
         pmsg += 1
      end
   end

   if colMin | colMax
      # ASSUMING ONE DOMAIN PER RANK, CONSTANT BLOCK SIZE HERE
      sendCount = dy * dz

      if colMin
         # contiguous memory
         srcOffset = 0
         offset = pmsg * maxPlaneComm
         for field in fields
            for i in 0:(dz-1)
               for j in 0:(dy-1)
                  domain.commDataSend[offset+i*dy+j + 1] = field[srcOffset+i*dx*dy+j*dx + 1]
               end
            end
            offset += sendCount
         end
         idx = pmsg * maxPlaneComm + 1
         src = view(domain.commDataSend, idx:(idx+(xferFields * sendCount)))

         otherRank = myRank - 1
         req = MPI.Isend(MPI.Buffer(src), otherRank, msgType, comm)
         domain.sendRequest[pmsg+1] = req
         pmsg += 1
      end

      if colMax && doSend
         # contiguous memory
         srcOffset = dx - 1
         offset = pmsg * maxPlaneComm
         for field in fields
            for i in 0:(dz-1)
               for j in 0:(dy-1)
                  domain.commDataSend[offset+i*dy+j + 1] = field[srcOffset+j*dx + 1]
               end
            end
            offset += sendCount
         end
         idx = pmsg * maxPlaneComm + 1
         src = view(domain.commDataSend, idx:(idx+(xferFields * sendCount)))
         otherRank = myRank + 1
         req = MPI.Isend(MPI.Buffer(src), otherRank, msgType, comm)
         domain.sendRequest[pmsg+1] = req
         pmsg += 1
      end
   end
   if !planeOnly
      if rowMin && colMin
         offset = pmsg * maxPlaneComm + emsg * maxEdgeComm
         for field in fields
            for i in 0:(dz-1)
               domain.commDataSend[offset+i + 1] = field[i*dx*dy + 1]
            end
            offset += dz
         end
         idx = pmsg * maxPlaneComm + emsg * maxEdgeComm + 1
         src = view(domain.commDataSend, idx:(idx+(xferFields * dz)))
         otherRank = myRank - domain.m_tp - 1

         req = MPI.Isend(MPI.Buffer(src), otherRank, msgType, comm)
         domain.sendRequest[pmsg+emsg+1] = req
         emsg += 1
      end

      if rowMin && planeMin
         offset = pmsg * maxPlaneComm + emsg * maxEdgeComm
         for field in fields
            for i in 0:(dx-1)
               domain.commDataSend[offset+i + 1] = field[i + 1]
            end
            offset += dx
         end
         idx = pmsg * maxPlaneComm + emsg * maxEdgeComm + 1
         src = view(domain.commDataSend, idx:(idx+(xferFields * dx)))
         otherRank = myRank - domain.m_tp^2 - domain.m_tp
         req = MPI.Isend(MPI.Buffer(src), otherRank, msgType, comm)
         domain.sendRequest[pmsg+emsg+1] = req
         emsg += 1
      end

      if colMin && planeMin
         offset = pmsg * maxPlaneComm + emsg * maxEdgeComm
         for field in fields
            for i in 0:(dy-1)
               domain.commDataSend[offset+i + 1] = field[i*dx + 1]
            end
            offset += dy
         end
         idx = pmsg * maxPlaneComm + emsg * maxEdgeComm + 1
         src = view(domain.commDataSend, idx:(idx+(xferFields * dy)))
         otherRank = myRank - domain.m_tp^2 - 1
         req = MPI.Isend(MPI.Buffer(src), otherRank, msgType, comm)
         domain.sendRequest[pmsg+emsg+1] = req
         emsg += 1
      end

      if rowMax && colMax && doSend
         offset = pmsg * maxPlaneComm + emsg * maxEdgeComm
         srcOffset = dx*dy - 1
         for field in fields
            for i in 0:(dz-1)
               domain.commDataSend[offset+i + 1] = field[srcOffset+i*dx*dy + 1]
            end
            offset += dz
         end
         idx = pmsg * maxPlaneComm + emsg * maxEdgeComm + 1
         src = view(domain.commDataSend, idx:(idx+(xferFields * dz)))
         otherRank = myRank + domain.m_tp + 1
         req = MPI.Isend(MPI.Buffer(src), otherRank, msgType, comm)
         domain.sendRequest[pmsg+emsg+1] = req
         emsg += 1
      end

      if rowMax && planeMax && doSend
         offset = pmsg * maxPlaneComm + emsg * maxEdgeComm
         srcOffset = dx*(dy-1) + dx*dy*(dz-1)
         for field in fields
            for i in 0:(dx-1)
               domain.commDataSend[offset+i + 1] = field[srcOffset+i + 1]
            end
            offset += dx
         end
         idx = pmsg * maxPlaneComm + emsg * maxEdgeComm + 1
         src = view(domain.commDataSend, idx:(idx+(xferFields * dx)))
         otherRank = myRank + domain.m_tp^2 + domain.m_tp
         req = MPI.Isend(MPI.Buffer(src), otherRank, msgType, comm)
         domain.sendRequest[pmsg+emsg+1] = req
         emsg += 1
      end

      if colMax && planeMax && doSend
         offset = pmsg * maxPlaneComm + emsg * maxEdgeComm
         srcOffset = dx*dy*(dz-1) + dx - 1
         for field in fields
            for i in 0:(dy-1)
               domain.commDataSend[offset+i + 1] = field[srcOffset+i*dx + 1]
            end
            offset += dy
         end
         idx = pmsg * maxPlaneComm + emsg * maxEdgeComm + 1
         src = view(domain.commDataSend, idx:(idx+(xferFields * dy)))
         otherRank = myRank + domain.m_tp^2 + 1
         req = MPI.Isend(MPI.Buffer(src), otherRank, msgType, comm)
         domain.sendRequest[pmsg+emsg+1] = req
         emsg += 1
      end

      if rowMax && colMin && doSend
         offset = pmsg * maxPlaneComm + emsg * maxEdgeComm
         srcOffset = dx*(dy-1)
         for field in fields
            for i in 0:(dz-1)
               domain.commDataSend[offset+i + 1] = field[srcOffset+i*dx*dy + 1]
            end
            offset += dz
         end
         idx = pmsg * maxPlaneComm + emsg * maxEdgeComm + 1
         src = view(domain.commDataSend, idx:(idx+(xferFields * dz)))
         otherRank = myRank + domain.m_tp - 1
         req = MPI.Isend(MPI.Buffer(src), otherRank, msgType, comm)
         domain.sendRequest[pmsg+emsg+1] = req
         emsg += 1
      end

      if rowMin && planeMax && doSend
         offset = pmsg * maxPlaneComm + emsg * maxEdgeComm
         srcOffset = dx*dy*(dz-1)
         for field in fields
            for i in 0:(dx-1)
               domain.commDataSend[offset+i + 1] = field[srcOffset+i + 1]
            end
            offset += dx
         end
         idx = pmsg * maxPlaneComm + emsg * maxEdgeComm + 1
         src = view(domain.commDataSend, idx:(idx+(xferFields * dx)))
         otherRank = myRank + domain.m_tp^2 - domain.m_tp
         req = MPI.Isend(MPI.Buffer(src), otherRank, msgType, comm)
         domain.sendRequest[pmsg+emsg+1] = req
         emsg += 1
      end

      if colMin && planeMax && doSend
         offset = pmsg * maxPlaneComm + emsg * maxEdgeComm
         srcOffset = dx*dy*(dz-1)
         for field in fields
            for i in 0:(dy-1)
               domain.commDataSend[offset+i + 1] = field[srcOffset+i*dx + 1]
            end
            offset += dy
         end
         idx = pmsg * maxPlaneComm + emsg * maxEdgeComm + 1
         src = view(domain.commDataSend, idx:(idx+(xferFields * dy)))
         otherRank = myRank + domain.m_tp^2 - 1
         req = MPI.Isend(MPI.Buffer(src), otherRank, msgType, comm)
         domain.sendRequest[pmsg+emsg+1] = req
         emsg += 1
      end

      if rowMin && colMax
         offset = pmsg * maxPlaneComm + emsg * maxEdgeComm
         srcOffset = dx - 1
         for field in fields
            for i in 0:(dz-1)
               domain.commDataSend[offset+i + 1] = field[srcOffset+i*dx*dy + 1]
            end
            offset += dz
         end
         idx = pmsg * maxPlaneComm + emsg * maxEdgeComm + 1
         src = view(domain.commDataSend, idx:(idx+(xferFields * dz)))
         otherRank = myRank - domain.m_tp + 1
         req = MPI.Isend(MPI.Buffer(src), otherRank, msgType, comm)
         domain.sendRequest[pmsg+emsg+1] = req
         emsg += 1
      end

      if rowMax && planeMin
         offset = pmsg * maxPlaneComm + emsg * maxEdgeComm
         srcOffset = dx*(dy - 1)
         for field in fields
            for i in 0:(dx-1)
               domain.commDataSend[offset+i + 1] = field[srcOffset+i + 1]
            end
            offset += dx
         end
         idx = pmsg * maxPlaneComm + emsg * maxEdgeComm + 1
         src = view(domain.commDataSend, idx:(idx+(xferFields * dx)))
         otherRank = myRank - domain.m_tp^2 + domain.m_tp
         req = MPI.Isend(MPI.Buffer(src), otherRank, msgType, comm)
         domain.sendRequest[pmsg+emsg+1] = req
         emsg += 1
      end

      if colMax && planeMin
         offset = pmsg * maxPlaneComm + emsg * maxEdgeComm
         srcOffset = dx - 1
         for field in fields
            for i in 0:(dy-1)
               domain.commDataSend[offset+i + 1] = field[srcOffset+i*dx + 1]
            end
            offset += dy
         end
         idx = pmsg * maxPlaneComm + emsg * maxEdgeComm + 1
         src = view(domain.commDataSend, idx:(idx+(xferFields * dy)))
         otherRank = myRank - domain.m_tp^2 + 1
         req = MPI.Isend(MPI.Buffer(src), otherRank, msgType, comm)
         domain.sendRequest[pmsg+emsg+1] = req
         emsg += 1
      end

      if rowMin && colMin && planeMin
         # corner at domain logical coord (0, 0, 0) */
         offset = pmsg * maxPlaneComm + emsg * maxEdgeComm +  cmsg * CACHE_COHERENCE_PAD_REAL
         for (fi, field) in enumerate(fields)
            # fi is one based already
            domain.commDataSend[offset+fi] = field[1]
         end
         idx = pmsg * maxPlaneComm + emsg * maxEdgeComm + cmsg * CACHE_COHERENCE_PAD_REAL + 1
         src = view(domain.commDataSend, idx:(idx+xferFields))
         otherRank = myRank - domain.m_tp^2 - domain.m_tp - 1
         req = MPI.Isend(MPI.Buffer(src), otherRank, msgType, comm)
         domain.sendRequest[pmsg+emsg+cmsg+1] = req
         cmsg += 1
      end

      if rowMin && colMin && planeMax && doSend
         # corner at domain logical coord (0, 0, 1)
         offset = pmsg * maxPlaneComm + emsg * maxEdgeComm + cmsg * CACHE_COHERENCE_PAD_REAL
         for (fi, field) in enumerate(fields)
            # fi is one based already
            domain.commDataSend[offset+fi] = field[dx*dy*(dz-1) + 1]
         end
         idx = pmsg * maxPlaneComm + emsg * maxEdgeComm + cmsg * CACHE_COHERENCE_PAD_REAL + 1
         src = view(domain.commDataSend, idx:(idx+xferFields))
         otherRank = myRank + domain.m_tp^2 - domain.m_tp - 1
         req = MPI.Isend(MPI.Buffer(src), otherRank, msgType, comm)
         domain.sendRequest[pmsg+emsg+cmsg+1] = req
         cmsg += 1
      end

      if rowMin && colMax && planeMin
         # corner at domain logical coord (1, 0, 0)
         offset = pmsg * maxPlaneComm + emsg * maxEdgeComm + cmsg * CACHE_COHERENCE_PAD_REAL
         for (fi, field) in enumerate(fields)
            # fi is one based already
            domain.commDataSend[offset+fi] = field[(dx - 1) + 1]
         end
         idx = pmsg * maxPlaneComm + emsg * maxEdgeComm + cmsg * CACHE_COHERENCE_PAD_REAL + 1
         src = view(domain.commDataSend, idx:(idx+xferFields))
         otherRank = myRank - domain.m_tp^2 - domain.m_tp + 1
         req = MPI.Isend(MPI.Buffer(src), otherRank, msgType, comm)
         domain.sendRequest[pmsg+emsg+cmsg+1] = req
         cmsg += 1
      end

      if rowMin && colMax && planeMax && doSend
         # corner at domain logical coord (1, 0, 1)
         offset = pmsg * maxPlaneComm + emsg * maxEdgeComm + cmsg * CACHE_COHERENCE_PAD_REAL
         for (fi, field) in enumerate(fields)
            # fi is one based already
            domain.commDataSend[offset+fi] = field[dx*dy*(dz - 1) + (dx - 1) + 1]
         end
         idx = pmsg * maxPlaneComm + emsg * maxEdgeComm + cmsg * CACHE_COHERENCE_PAD_REAL + 1
         src = view(domain.commDataSend, idx:(idx+xferFields))
         otherRank = myRank + domain.m_tp^2 - domain.m_tp + 1
         req = MPI.Isend(MPI.Buffer(src), otherRank, msgType, comm)
         domain.sendRequest[pmsg+emsg+cmsg+1] = req
         cmsg += 1
      end

      if rowMax && colMin && planeMin
         # corner at domain logical coord (0, 1, 0)
         offset = pmsg * maxPlaneComm + emsg * maxEdgeComm + cmsg * CACHE_COHERENCE_PAD_REAL
         for (fi, field) in enumerate(fields)
            # fi is one based already
            domain.commDataSend[offset+fi] = field[dx*(dy - 1) + 1]
         end
         idx = pmsg * maxPlaneComm + emsg * maxEdgeComm + cmsg * CACHE_COHERENCE_PAD_REAL + 1
         src = view(domain.commDataSend, idx:(idx+xferFields))
         otherRank = myRank - domain.m_tp^2 + domain.m_tp - 1
         req = MPI.Isend(MPI.Buffer(src), otherRank, msgType, comm)
         domain.sendRequest[pmsg+emsg+cmsg+1] = req
         cmsg += 1
      end

      if rowMax && colMin && planeMax && doSend
         # corner at domain logical coord (0, 1, 1)
         offset = pmsg * maxPlaneComm + emsg * maxEdgeComm + cmsg * CACHE_COHERENCE_PAD_REAL
         for (fi, field) in enumerate(fields)
            # fi is one based already
            domain.commDataSend[offset+fi] = field[dx*dy*(dz - 1) + dx*(dy - 1)+ 1]
         end
         idx = pmsg * maxPlaneComm + emsg * maxEdgeComm + cmsg * CACHE_COHERENCE_PAD_REAL + 1
         src = view(domain.commDataSend, idx:(idx+xferFields))
         otherRank = myRank + domain.m_tp^2 + domain.m_tp - 1
         req = MPI.Isend(MPI.Buffer(src), otherRank, msgType, comm)
         domain.sendRequest[pmsg+emsg+cmsg+1] = req
         cmsg += 1
      end

      if rowMax && colMax && planeMin
         # corner at domain logical coord (1, 1, 0)
         offset = pmsg * maxPlaneComm + emsg * maxEdgeComm + cmsg * CACHE_COHERENCE_PAD_REAL
         for (fi, field) in enumerate(fields)
            # fi is one based already
            domain.commDataSend[offset+fi] = field[dx*dy - 1 + 1]
         end
         idx = pmsg * maxPlaneComm + emsg * maxEdgeComm + cmsg * CACHE_COHERENCE_PAD_REAL + 1
         src = view(domain.commDataSend, idx:(idx+xferFields))
         otherRank = myRank - domain.m_tp^2 + domain.m_tp + 1
         req = MPI.Isend(MPI.Buffer(src), otherRank, msgType, comm)
         domain.sendRequest[pmsg+emsg+cmsg+1] = req
         cmsg += 1
      end

      if rowMax && colMax && planeMax && doSend
         # corner at domain logical coord (1, 1, 1)
         offset = pmsg * maxPlaneComm + emsg * maxEdgeComm + cmsg * CACHE_COHERENCE_PAD_REAL
         for (fi, field) in enumerate(fields)
            # fi is one based already
            domain.commDataSend[offset+fi] = field[dx*dy*dz - 1 + 1]
         end
         idx = pmsg * maxPlaneComm + emsg * maxEdgeComm + cmsg * CACHE_COHERENCE_PAD_REAL + 1
         src = view(domain.commDataSend, idx:(idx+xferFields))
         otherRank = myRank + domain.m_tp^2 + domain.m_tp + 1
         req = MPI.Isend(MPI.Buffer(src), otherRank, msgType, comm)
         domain.sendRequest[pmsg+emsg+cmsg+1] = req
         cmsg += 1
      end
   end

   for i in 1:(pmsg+emsg+cmsg)
      MPI.Wait!(domain.sendRequest[i])
   end
   # MPI.Waitall!(domain.sendRequest)
end

function commSBN(domain::Domain, fields)
   comm = domain.comm
   if comm === nothing
      return
   end

   xferFields = length(fields)

   maxPlaneComm = xferFields * domain.maxPlaneSize
   maxEdgeComm  = xferFields * domain.maxEdgeSize
   pmsg = 0 # plane comm msg
   emsg = 0 # edge comm msg
   cmsg = 0 # corner comm msg

   dx = domain.sizeX + 1
   dy = domain.sizeY + 1
   dz = domain.sizeZ + 1

   rowMin, rowMax, colMin, colMax, planeMin, planeMax = get_neighbors(domain)

   myRank = MPI.Comm_rank(comm)

   if planeMin || planeMax
      opCount = dx * dy ;

      if planeMin
         srcaddr = pmsg * maxPlaneComm;
         status = MPI.Wait!(domain.recvRequest[pmsg + 1]);
         for field in fields
            for i in 1:opCount
               field[i] += domain.commDataRecv[srcaddr+i]
            end
            srcaddr += opCount
         end
         pmsg += 1
      end
      if planeMax
         srcaddr = pmsg * maxPlaneComm;
         status = MPI.Wait!(domain.recvRequest[pmsg + 1]);
         for field in fields
            for i in 1:opCount
               field[dx*dy*(dz-1) + i] += domain.commDataRecv[srcaddr+i]
            end
            srcaddr += opCount
         end
         pmsg += 1
      end
   end

   if rowMin || rowMax
      opCount = dx * dz

      if rowMin
         srcaddr = pmsg * maxPlaneComm;
         status = MPI.Wait!(domain.recvRequest[pmsg + 1]);
         for field in fields
            for i in 0:(dz-1)
               for j in 0:(dx-1)
                  field[i*dx*dy + j + 1] += domain.commDataRecv[srcaddr + i*dx + j + 1]
               end
            end
            srcaddr += opCount
         end
         pmsg += 1
      end

      if rowMax
         srcaddr = pmsg * maxPlaneComm;
         status = MPI.Wait!(domain.recvRequest[pmsg + 1]);
         for field in fields
            for i in 0:(dz-1)
               for j in 0:(dx-1)
                  field[dx*(dy-1) + i*dx*dy + j + 1] += domain.commDataRecv[srcaddr + i*dx + j + 1]
               end
            end
            srcaddr += opCount
         end
         pmsg += 1
      end
   end

   if colMin || colMax
      opCount = dy * dz ;

      if colMin
         srcaddr = pmsg * maxPlaneComm;
         status = MPI.Wait!(domain.recvRequest[pmsg + 1]);
         for field in fields
            for i in 0:(dz-1)
               for j in 0:(dy-1)
                  field[i*dx*dy + j*dx + 1] += domain.commDataRecv[srcaddr + i*dy + j + 1]
               end
            end
            srcaddr += opCount
         end
         pmsg += 1
      end

      if colMax
         srcaddr = pmsg * maxPlaneComm;
         status = MPI.Wait!(domain.recvRequest[pmsg + 1]);
         for field in fields
            for i in 0:(dz-1)
               for j in 0:(dy-1)
                  field[dx - 1 + i*dx*dy + j*dx + 1] += domain.commDataRecv[srcaddr + i*dy + j + 1]
               end
            end
            srcaddr += opCount
         end
         pmsg += 1
      end
   end

   if rowMin && colMin
      srcaddr = pmsg * maxPlaneComm + emsg * maxEdgeComm
      status = MPI.Wait!(domain.recvRequest[pmsg+emsg + 1])
      for field in fields
         for i in 0:(dz-1)
            # (domain.*dest)(i*dx*dy) += srcaddr[i] ;
            field[i*dx*dy + 1] += domain.commDataRecv[srcaddr + i + 1]
         end
         srcaddr += dz
      end
      emsg += 1
   end

   if rowMin && planeMin
      srcaddr = pmsg * maxPlaneComm + emsg * maxEdgeComm
      status = MPI.Wait!(domain.recvRequest[pmsg+emsg + 1])
      for field in fields
         for i in 0:(dx-1)
            # (domain.*dest)(i) += srcaddr[i] ;
            field[i + 1] += domain.commDataRecv[srcaddr + i + 1]
         end
         srcaddr += dx
      end
      emsg += 1
   end

   if colMin && planeMin
      srcaddr = pmsg * maxPlaneComm + emsg * maxEdgeComm
      status = MPI.Wait!(domain.recvRequest[pmsg+emsg + 1])
      for field in fields
         for i in 0:(dy-1)
            # (domain.*dest)(i*dx) += srcaddr[i] ;
            field[i*dx + 1] += domain.commDataRecv[srcaddr + i + 1]
         end
         srcaddr += dy
      end
      emsg += 1
   end

   if rowMax && colMax
      srcaddr = pmsg * maxPlaneComm + emsg * maxEdgeComm
      status = MPI.Wait!(domain.recvRequest[pmsg+emsg + 1])
      for field in fields
         for i in 0:(dz-1)
            # (domain.*dest)(dx*dy - 1 + i*dx*dy) += srcaddr[i] ;
            field[dx*dy - 1 + i*dx*dy + 1] += domain.commDataRecv[srcaddr + i + 1]
         end
         srcaddr += dz
      end
      emsg += 1
   end

   if rowMax && planeMax
      srcaddr = pmsg * maxPlaneComm + emsg * maxEdgeComm
      status = MPI.Wait!(domain.recvRequest[pmsg+emsg + 1])
      for field in fields
         for i in 0:(dx-1)
            # (domain.*dest)(dx*(dy-1) + dx*dy*(dz-1) + i) += srcaddr[i] ;
            field[dx*(dy-1) + dx*dy*(dz-1) + i + 1] += domain.commDataRecv[srcaddr + i + 1]
         end
         srcaddr += dx
      end
      emsg += 1
   end

   if colMax && planeMax
      srcaddr = pmsg * maxPlaneComm + emsg * maxEdgeComm
      status = MPI.Wait!(domain.recvRequest[pmsg+emsg + 1])
      for field in fields
         for i in 0:(dy-1)
            # (domain.*dest)(dx*dy*(dz-1) + dx - 1 + i*dx) += srcaddr[i] ;
            field[dx*dy*(dz-1) + dx - 1 + i*dx + 1] += domain.commDataRecv[srcaddr + i + 1]
         end
         srcaddr += dy
      end
      emsg += 1
   end

   if rowMax && colMin
      srcaddr = pmsg * maxPlaneComm + emsg * maxEdgeComm
      status = MPI.Wait!(domain.recvRequest[pmsg+emsg + 1])
      for field in fields
         for i in 0:(dz-1)
            # (domain.*dest)(dx*(dy-1) + i*dx*dy) += srcaddr[i] ;
            field[dx*(dy-1) + i*dx*dy + 1] += domain.commDataRecv[srcaddr + i + 1]
         end
         srcaddr += dz
      end
      emsg += 1
   end

   if rowMin && planeMax
      srcaddr = pmsg * maxPlaneComm + emsg * maxEdgeComm
      status = MPI.Wait!(domain.recvRequest[pmsg+emsg + 1])
      for field in fields
         for i in 0:(dx-1)
            # (domain.*dest)(dx*dy*(dz-1) + i) += srcaddr[i] ;
            field[dx*dy*(dz-1) + i + 1] += domain.commDataRecv[srcaddr + i + 1]
         end
         srcaddr += dx
      end
      emsg += 1
   end

   if colMin && planeMax
      srcaddr = pmsg * maxPlaneComm + emsg * maxEdgeComm
      status = MPI.Wait!(domain.recvRequest[pmsg+emsg + 1])
      for field in fields
         for i in 0:(dy-1)
            # (domain.*dest)(dx*dy*(dz-1) + i*dx) += srcaddr[i] ;
            field[dx*dy*(dz-1) + i*dx + 1] += domain.commDataRecv[srcaddr + i + 1]
         end
         srcaddr += dy
      end
      emsg += 1
   end

   if rowMin && colMax
      srcaddr = pmsg * maxPlaneComm + emsg * maxEdgeComm
      status = MPI.Wait!(domain.recvRequest[pmsg+emsg + 1])
      for field in fields
         for i in 0:(dz-1)
            # (domain.*dest)(dx - 1 + i*dx*dy) += srcaddr[i] ;
            field[dx - 1 + i*dx*dy + 1] += domain.commDataRecv[srcaddr + i + 1]
         end
         srcaddr += dz
      end
      emsg += 1
   end

   if rowMax && planeMin
      srcaddr = pmsg * maxPlaneComm + emsg * maxEdgeComm
      status = MPI.Wait!(domain.recvRequest[pmsg+emsg + 1])
      for field in fields
         for i in 0:(dx-1)
            # (domain.*dest)(dx*(dy - 1) + i) += srcaddr[i] ;
            field[dx*(dy-1) + i + 1] += domain.commDataRecv[srcaddr + i + 1]
         end
         srcaddr += dx
      end
      emsg += 1
   end

   if colMax && planeMin
      srcaddr = pmsg * maxPlaneComm + emsg * maxEdgeComm
      status = MPI.Wait!(domain.recvRequest[pmsg+emsg + 1])
      for field in fields
         for i in 0:(dy-1)
            # (domain.*dest)(dx - 1 + i*dx) += srcaddr[i] ;
            field[dx - 1 + i*dx + 1] += domain.commDataRecv[srcaddr + i + 1]
         end
         srcaddr += dy
      end
      emsg += 1
   end

   if rowMin && colMin && planeMin
      idx = 0
      comBuf = pmsg * maxPlaneComm + emsg * maxEdgeComm + cmsg * CACHE_COHERENCE_PAD_REAL
      status = MPI.Wait!(domain.recvRequest[pmsg+emsg+cmsg + 1])
      for (fi, field) in enumerate(fields)
         field[idx+1] += domain.commDataRecv[comBuf + fi] # No 1 here as already taken care of by fi
      end
      cmsg += 1
   end

   if rowMin && colMin && planeMax
      idx = dx*dy*(dz - 1)
      comBuf = pmsg * maxPlaneComm + emsg * maxEdgeComm + cmsg * CACHE_COHERENCE_PAD_REAL
      status = MPI.Wait!(domain.recvRequest[pmsg+emsg+cmsg + 1])
      for (fi, field) in enumerate(fields)
         field[idx+1] += domain.commDataRecv[comBuf + fi] # No 1 here as already taken care of by fi
      end
      cmsg += 1
   end

   if rowMin && colMax && planeMin
      idx = dx - 1
      comBuf = pmsg * maxPlaneComm + emsg * maxEdgeComm + cmsg * CACHE_COHERENCE_PAD_REAL
      status = MPI.Wait!(domain.recvRequest[pmsg+emsg+cmsg + 1])
      for (fi, field) in enumerate(fields)
         field[idx+1] += domain.commDataRecv[comBuf + fi] # No 1 here as already taken care of by fi
      end
      cmsg += 1
   end

   if rowMin && colMax && planeMax
      idx = dx*dy*(dz - 1) + (dx - 1)
      comBuf = pmsg * maxPlaneComm + emsg * maxEdgeComm + cmsg * CACHE_COHERENCE_PAD_REAL
      status = MPI.Wait!(domain.recvRequest[pmsg+emsg+cmsg + 1])
      for (fi, field) in enumerate(fields)
         field[idx+1] += domain.commDataRecv[comBuf + fi] # No 1 here as already taken care of by fi
      end
      cmsg += 1
   end

   if rowMax && colMin && planeMin
      idx = dx*(dy - 1)
      comBuf = pmsg * maxPlaneComm + emsg * maxEdgeComm + cmsg * CACHE_COHERENCE_PAD_REAL
      status = MPI.Wait!(domain.recvRequest[pmsg+emsg+cmsg + 1])
      for (fi, field) in enumerate(fields)
         field[idx+1] += domain.commDataRecv[comBuf + fi] # No 1 here as already taken care of by fi
      end
      cmsg += 1
   end

   if rowMax && colMin && planeMax
      idx = dx*dy*(dz - 1) + dx*(dy - 1)
      comBuf = pmsg * maxPlaneComm + emsg * maxEdgeComm + cmsg * CACHE_COHERENCE_PAD_REAL
      status = MPI.Wait!(domain.recvRequest[pmsg+emsg+cmsg + 1])
      for (fi, field) in enumerate(fields)
         field[idx+1] += domain.commDataRecv[comBuf + fi] # No 1 here as already taken care of by fi
      end
      cmsg += 1
   end

   if rowMax && colMax && planeMin
      idx = dx*dy - 1
      comBuf = pmsg * maxPlaneComm + emsg * maxEdgeComm + cmsg * CACHE_COHERENCE_PAD_REAL
      status = MPI.Wait!(domain.recvRequest[pmsg+emsg+cmsg + 1])
      for (fi, field) in enumerate(fields)
         field[idx+1] += domain.commDataRecv[comBuf + fi] # No 1 here as already taken care of by fi
      end
      cmsg += 1
   end

   if rowMax && colMax && planeMax
      idx = dx*dy*dz - 1
      comBuf = pmsg * maxPlaneComm + emsg * maxEdgeComm + cmsg * CACHE_COHERENCE_PAD_REAL
      status = MPI.Wait!(domain.recvRequest[pmsg+emsg+cmsg + 1])
      for (fi, field) in enumerate(fields)
         field[idx+1] += domain.commDataRecv[comBuf + fi] # No 1 here as already taken care of by fi
      end
      cmsg += 1
   end
end

function commMonoQ(domain::Domain)
   comm = domain.comm
   if comm === nothing
       return
   end

   xferFields = 3 # delv_xi, delv_eta, delv_zeta
   fields = (domain.delv_xi, domain.delv_eta, domain.delv_zeta)
   maxPlaneComm = xferFields * domain.maxPlaneSize

   pmsg = 0 # plane comm msg
   dx = domain.sizeX
   dy = domain.sizeY
   dz = domain.sizeZ

   # assume communication to 6 neighbors by default
   rowMin, rowMax, colMin, colMax, planeMin, planeMax = get_neighbors(domain)

   numElem = domain.numElem

   # point into ghost data area
   fieldOffsets = [numElem, numElem, numElem]

   myRank = MPI.Comm_rank(comm)

   if planeMin || planeMax
      # ASSUMING ONE DOMAIN PER RANK, CONSTANT BLOCK SIZE HERE */
      opCount = dx * dy

      if planeMin
         # contiguous memory
         MPI.Wait!(domain.recvRequest[pmsg+1])

         offset = pmsg * maxPlaneComm
         for (fi, field) in enumerate(fields)
            fOffset = fieldOffsets[fi]
            copyto_zero!(field, fOffset, domain.commDataRecv, offset, opCount)

            offset += opCount
            fieldOffsets[fi] += opCount
         end
         pmsg += 1
      end
      if planeMax
         # contiguous memory
         MPI.Wait!(domain.recvRequest[pmsg+1])
         offset = pmsg * maxPlaneComm
         for (fi, field) in enumerate(fields)
            fOffset = fieldOffsets[fi]
            copyto_zero!(field, fOffset, domain.commDataRecv, offset, opCount)

            offset += opCount
            fieldOffsets[fi] += opCount
         end
         pmsg += 1
      end
   end

   if rowMin || rowMax
      # ASSUMING ONE DOMAIN PER RANK, CONSTANT BLOCK SIZE HERE
      opCount = dx * dz

      if rowMin
         # contiguous memory
         MPI.Wait!(domain.recvRequest[pmsg+1])
         offset = pmsg * maxPlaneComm
         for (fi, field) in enumerate(fields)
            fOffset = fieldOffsets[fi]
            copyto_zero!(field, fOffset, domain.commDataRecv, offset, opCount)

            offset += opCount
            fieldOffsets[fi] += opCount
         end
         pmsg += 1
      end

      if rowMax
         # contiguous memory
         MPI.Wait!(domain.recvRequest[pmsg+1])
         offset = pmsg * maxPlaneComm
         for (fi, field) in enumerate(fields)
            fOffset = fieldOffsets[fi]
            copyto_zero!(field, fOffset, domain.commDataRecv, offset, opCount)
            offset += opCount
            fieldOffsets[fi] += opCount
         end
         pmsg += 1
      end
   end
   if colMin || colMax
      # ASSUMING ONE DOMAIN PER RANK, CONSTANT BLOCK SIZE HERE
      opCount = dy * dz

      if colMin
         # contiguous memory
         MPI.Wait!(domain.recvRequest[pmsg+1])
         offset = pmsg * maxPlaneComm
         for (fi, field) in enumerate(fields)
            fOffset = fieldOffsets[fi]
            copyto_zero!(field, fOffset, domain.commDataRecv, offset, opCount)

            offset += opCount
            fieldOffsets[fi] += opCount
         end
         pmsg += 1
      end
      if colMax
         # contiguous memory
         MPI.Wait!(domain.recvRequest[pmsg+1])
         offset = pmsg * maxPlaneComm
         for (fi, field) in enumerate(fields)
            fOffset = fieldOffsets[fi]
            copyto_zero!(field, fOffset, domain.commDataRecv, offset, opCount)

            offset += opCount
            fieldOffsets[fi] += opCount
         end
         pmsg += 1
      end
   end
   return nothing
end

function commSyncPosVel(domain::Domain)
   comm = domain.comm
   if comm === nothing
      return
   end
   doRecv = false
   xferFields = 6 ; # x, y, z, xd, yd, zd
   fields = (domain.x, domain.y, domain.z, domain.xd, domain.yd, domain.zd)
   maxPlaneComm = xferFields * domain.maxPlaneSize
   maxEdgeComm  = xferFields * domain.maxEdgeSize
   pmsg = 0 # plane comm msg
   emsg = 0 # edge comm msg
   cmsg = 0 # corner comm msg
   dx = domain.sizeX + 1
   dy = domain.sizeY + 1
   dz = domain.sizeZ + 1

   # assume communication to 6 neighbors by default
   rowMin, rowMax, colMin, colMax, planeMin, planeMax = get_neighbors(domain)

   myRank = MPI.Comm_rank(comm)

   if planeMin | planeMax
      # ASSUMING ONE DOMAIN PER RANK, CONSTANT BLOCK SIZE HERE
      opCount = dx * dy

      if planeMin && doRecv
         # contiguous memory
         MPI.Wait!(domain.recvRequest[pmsg+1])
         offset = pmsg * maxPlaneComm
         for field in fields
            copyto_zero!(field, 0, domain.commDataRecv, offset, opCount)
            offset += opCount
         end
         pmsg += 1
      end

      if planeMax
         # contiguous memory
         MPI.Wait!(domain.recvRequest[pmsg+1])

         destOffset = dx*dy*(dz - 1)
         offset = pmsg * maxPlaneComm
         for field in fields
            copyto_zero!(field, destOffset, domain.commDataRecv, offset, opCount)
            offset += opCount
         end
         pmsg += 1
      end

   end
   if rowMin | rowMax
      # ASSUMING ONE DOMAIN PER RANK, CONSTANT BLOCK SIZE HERE
      opCount = dx * dz

      if rowMin && doRecv
         # contiguous memory
         MPI.Wait!(domain.recvRequest[pmsg+1])
         offset = pmsg * maxPlaneComm
         for field in fields
            for i in 0:(dz-1)
               destOffset = i*dx*dy
               copyto_zero!(field, destOffset, domain.commDataRecv, offset + i*dx, dx)
            end
            offset += opCount
         end
         pmsg += 1
      end
      if rowMax
         # contiguous memory
         MPI.Wait!(domain.recvRequest[pmsg+1])
         offset = pmsg * maxPlaneComm
         for field in fields
            for i in 0:(dz-1)
               destOffset = dx*(dy - 1) + i*dx*dy
               for j in 0:(dx-1)
                  field[destOffset + j + 1] = domain.commDataRecv[offset + i*dx + j + 1]
               end
            end
            offset += opCount
         end
         pmsg += 1
      end
   end
   if colMin | colMax
      # ASSUMING ONE DOMAIN PER RANK, CONSTANT BLOCK SIZE HERE
      opCount = dy * dz

      if colMin && doRecv
         # contiguous memory
         MPI.Wait!(domain.recvRequest[pmsg+1])
         offset = pmsg * maxPlaneComm
         for field in fields
            for i in 0:(dz-1)
               for j in 0:(dy-1)
                  field[i*dx*dy+j*dx + 1] = domain.commDataRecv[offset + i*dy + j + 1]
               end
            end
            offset += opCount
         end
         pmsg += 1
      end
      if colMax
         # contiguous memory
         MPI.Wait!(domain.recvRequest[pmsg+1])
         offset = pmsg * maxPlaneComm
         destOffset = dx-1
         for field in fields
            for i in 0:(dz-1)
               for j in 0:(dy-1)
                  field[destOffset + i*dx*dy + j*dx + 1] = domain.commDataRecv[offset + i*dy + j + 1]
               end
            end
            offset += opCount
         end
         pmsg += 1
      end
   end

   if rowMin && colMin && doRecv
      MPI.Wait!(domain.recvRequest[pmsg+emsg+1])
      offset = pmsg * maxPlaneComm + emsg * maxEdgeComm
      for field in fields
         for i in 0:(dz-1)
            field[i*dx*dy + 1] = domain.commDataRecv[offset + i + 1]
         end
         offset += dz
      end
      emsg += 1
   end

   if rowMin && planeMin && doRecv
      MPI.Wait!(domain.recvRequest[pmsg+emsg+1])
      offset = pmsg * maxPlaneComm + emsg * maxEdgeComm
      for field in fields
         copyto_zero!(field, 0, domain.commDataRecv, offset, dx)
         offset += dx
      end
      emsg += 1
   end

   if colMin && planeMin && doRecv
      MPI.Wait!(domain.recvRequest[pmsg+emsg+1])
      offset = pmsg * maxPlaneComm + emsg * maxEdgeComm
      for field in fields
         for i in 0:(dy-1)
            field[i*dx + 1] = domain.commDataRecv[offset + i + 1]
         end
         offset += dy
      end
      emsg += 1
   end

   if rowMax && colMax
      MPI.Wait!(domain.recvRequest[pmsg+emsg+1])
      offset = pmsg * maxPlaneComm + emsg * maxEdgeComm
      destOffset = dx*dy - 1
      for field in fields
         for i in 0:(dz-1)
            field[destOffset + i*dx*dy + 1] = domain.commDataRecv[offset + i + 1]
         end
         offset += dz
      end
      emsg += 1
   end

   if rowMax && planeMax
      MPI.Wait!(domain.recvRequest[pmsg+emsg+1])
      offset = pmsg * maxPlaneComm + emsg * maxEdgeComm
      destOffset = dx*(dy-1) + dx*dy*(dz-1)
      for field in fields
         copyto_zero!(field, destOffset, domain.commDataRecv, offset, dx)
         offset += dx
      end
      emsg += 1
   end

   if colMax && planeMax
      MPI.Wait!(domain.recvRequest[pmsg+emsg+1])
      offset = pmsg * maxPlaneComm + emsg * maxEdgeComm
      destOffset = dx*dy*(dz-1) + dx - 1
      for field in fields
         for i in 0:(dy-1)
            field[destOffset + i*dx + 1] = domain.commDataRecv[offset + i + 1]
         end
         offset += dy
      end
      emsg += 1
   end

   if rowMax && colMin
      MPI.Wait!(domain.recvRequest[pmsg+emsg+1])
      offset = pmsg * maxPlaneComm + emsg * maxEdgeComm
      destOffset = dx*(dy-1)
      for field in fields
         for i in 0:(dz-1)
            field[destOffset + i*dx*dy + 1] = domain.commDataRecv[offset + i + 1]
         end
         offset += dz
      end
      emsg += 1
   end

   if rowMin && planeMax
      MPI.Wait!(domain.recvRequest[pmsg+emsg+1])
      offset = pmsg * maxPlaneComm + emsg * maxEdgeComm
      destOffset = dx*dy*(dz-1)
      for field in fields
         copyto_zero!(field, destOffset, domain.commDataRecv, offset, dx)
         offset += dx
      end
      emsg += 1
   end

   if colMin && planeMax
      MPI.Wait!(domain.recvRequest[pmsg+emsg+1])
      offset = pmsg * maxPlaneComm + emsg * maxEdgeComm
      destOffset = dx*dy*(dz-1)
      for field in fields
         for i in 0:(dy-1)
            field[destOffset + i*dx + 1] = domain.commDataRecv[offset + i + 1]
         end
         offset += dy
      end
      emsg += 1
   end

   if rowMin && colMax && doRecv
      MPI.Wait!(domain.recvRequest[pmsg+emsg+1])
      offset = pmsg * maxPlaneComm + emsg * maxEdgeComm
      destOffset = dx-1
      for field in fields
         for i in 0:(dz-1)
            field[destOffset + i*dx*dy + 1] = domain.commDataRecv[offset + i + 1]
         end
         offset += dz
      end
      emsg += 1
   end

   if rowMax && planeMin && doRecv
      MPI.Wait!(domain.recvRequest[pmsg+emsg+1])
      offset = pmsg * maxPlaneComm + emsg * maxEdgeComm
      destOffset = dx*(dy - 1)
      for field in fields
         copyto_zero!(field, destOffset, domain.commDataRecv, offset, dx)
         offset += dx
      end
      emsg += 1
   end

   if colMax && planeMin && doRecv
      MPI.Wait!(domain.recvRequest[pmsg+emsg+1])
      offset = pmsg * maxPlaneComm + emsg * maxEdgeComm
      destOffset = dx - 1
      for field in fields
         for i in 0:(dy-1)
            field[destOffset + i*dx + 1] = domain.commDataRecv[offset + i + 1]
         end
         offset += dy
      end
      emsg += 1
   end

   if rowMin && colMin && planeMin && doRecv
      # corner at domain logical coord (0, 0, 0)
      MPI.Wait!(domain.recvRequest[pmsg+emsg+cmsg+1])
      offset = pmsg * maxPlaneComm + emsg * maxEdgeComm + cmsg * CACHE_COHERENCE_PAD_REAL
      for (fi, field) in enumerate(fields)
         # fi is one based
         field[1] = domain.commDataRecv[offset + fi]
      end
      cmsg += 1
   end

   if rowMin && colMin && planeMax
      # corner at domain logical coord (0, 0, 1)
      MPI.Wait!(domain.recvRequest[pmsg+emsg+cmsg+1])
      offset = pmsg * maxPlaneComm + emsg * maxEdgeComm + cmsg * CACHE_COHERENCE_PAD_REAL
      idx = dx*dy*(dz - 1)
      for (fi, field) in enumerate(fields)
         # fi is one based
         field[idx + 1] = domain.commDataRecv[offset + fi]
      end
      cmsg += 1
   end

   if rowMin && colMax && planeMin && doRecv
      # corner at domain logical coord (1, 0, 0)
      MPI.Wait!(domain.recvRequest[pmsg+emsg+cmsg+1])
      offset = pmsg * maxPlaneComm + emsg * maxEdgeComm + cmsg * CACHE_COHERENCE_PAD_REAL
      idx = dx - 1
      for (fi, field) in enumerate(fields)
         # fi is one based
         field[idx + 1] = domain.commDataRecv[offset + fi]
      end
      cmsg += 1
   end

   if rowMin && colMax && planeMax
      # corner at domain logical coord (1, 0, 1)
      MPI.Wait!(domain.recvRequest[pmsg+emsg+cmsg+1])
      offset = pmsg * maxPlaneComm + emsg * maxEdgeComm + cmsg * CACHE_COHERENCE_PAD_REAL
      idx = dx*dy*(dz - 1) + (dx - 1)
      for (fi, field) in enumerate(fields)
         # fi is one based
         field[idx + 1] = domain.commDataRecv[offset + fi]
      end
      cmsg += 1
   end

   if rowMax && colMin && planeMin && doRecv
      # corner at domain logical coord (0, 1, 0)
      MPI.Wait!(domain.recvRequest[pmsg+emsg+cmsg+1])
      offset = pmsg * maxPlaneComm + emsg * maxEdgeComm + cmsg * CACHE_COHERENCE_PAD_REAL
      idx = dx*(dy - 1)
      for (fi, field) in enumerate(fields)
         # fi is one based
         field[idx + 1] = domain.commDataRecv[offset + fi]
      end
      cmsg += 1
   end

   if rowMax && colMin && planeMax
      # corner at domain logical coord (0, 1, 1)
      MPI.Wait!(domain.recvRequest[pmsg+emsg+cmsg+1])
      offset = pmsg * maxPlaneComm + emsg * maxEdgeComm + cmsg * CACHE_COHERENCE_PAD_REAL
      idx = dx*dy*(dz - 1) + dx*(dy - 1)
      for (fi, field) in enumerate(fields)
         # fi is one based
         field[idx + 1] = domain.commDataRecv[offset + fi]
      end
      cmsg += 1
   end

   if rowMax && colMax && planeMin && doRecv
      # corner at domain logical coord (1, 1, 0)
      MPI.Wait!(domain.recvRequest[pmsg+emsg+cmsg+1])
      offset = pmsg * maxPlaneComm + emsg * maxEdgeComm + cmsg * CACHE_COHERENCE_PAD_REAL
      idx = dx*dy - 1
      for (fi, field) in enumerate(fields)
         # fi is one based
         field[idx + 1] = domain.commDataRecv[offset + fi]
      end
      cmsg += 1
   end

   if rowMax && colMax && planeMax
      # corner at domain logical coord (1, 1, 1)
      MPI.Wait!(domain.recvRequest[pmsg+emsg+cmsg+1])
      offset = pmsg * maxPlaneComm + emsg * maxEdgeComm + cmsg * CACHE_COHERENCE_PAD_REAL
      idx = dx*dy*dz - 1
      for (fi, field) in enumerate(fields)
         # fi is one based
         field[idx + 1] = domain.commDataRecv[offset + fi]
      end
      cmsg += 1
   end

   return nothing
end


text:
using Pkg

Pkg.add("Glob")
Pkg.add("HTTP")
Pkg.add("Gumbo")
Pkg.add("Cascadia")

Pkg.add("SparseArrays")
Pkg.add("LinearAlgebra")
Pkg.add("SharedArrays")
Pkg.add("StaticArrays")
Pkg.add("Rotations")
Pkg.add("IterativeSolvers")
Pkg.add("Arpack")
Pkg.add("LinearMaps")

Pkg.add("PyPlot")
Pkg.add("Printf")
Pkg.add("LaTeXStrings")

Pkg.add("HDF5")
Pkg.add("JLD")
Pkg.add("JSON")

Pkg.add("BenchmarkTools")

ext: jl

max_stars_repo_name: wilsonkl/sfmflex-release

max_issues_repo_name: wilsonkl/sfmflex-release

max_forks_repo_name: wilsonkl/sfmflex-release


text:
""" Adds one layer of depth to MIOTree."""
function add_depth(mt::MIOTree)
    all_nodes = allnodes(mt)
    all_leaves = [nd for nd in all_nodes if is_leaf(nd)]
    max_idx = maximum([nd.idx for nd in all_nodes])
    all_idxs = max_idx .+ collect(1:2*length(all_leaves))
    new_leaves = BinaryNode[]
    for leaf in all_leaves
        leftchild(leaf, BinaryNode(popfirst!(all_idxs)))
        push!(new_leaves, leaf.left)
        rightchild(leaf, BinaryNode(popfirst!(all_idxs)))
        push!(new_leaves, leaf.right)
    end
    return
end


text:
@testset "NALU" begin
    nalu = NALU(2, 3)
    x = rand(2, 10) .+ 1
    z = nalu(x)

    @test size(z) == (3,10)

    ps = params(nalu)
    @test length(ps) == 4

    loss() = Flux.mse(nalu(x), reshape(repeat(x[1,:] .* x[2,:], inner=(3,1)),3,10))
    opt = ADAM(1.)
    train_data = Iterators.repeated((), 2000)
    Flux.train!(loss, ps, train_data, opt)

    @test all(isapprox.(weights(nalu.nac), [1.0 1.0; 1.0 1.0; 1.0 1.0], atol=1e-2))
    @test all(isapprox.(gate(nalu, x[:,1]), [0.0, 0.0, 0.0], atol=1e-3))
end


text:
import Base.==

#Base.iszero(t::Term) = iszero(t.α)
Base.iszero(p::Polynomial) = isempty(p)

# TODO This should be in Base with T instead of PolyVar{C}.
# See https://github.com/blegat/MultivariatePolynomials.jl/issues/3
function (==)(x::Vector{PolyVar{C}}, y::Vector{PolyVar{C}}) where C
    if length(x) != length(y)
        false
    else
        #for (xi, yi) in zip(x, y)
        for i in 1:length(x)
            if x[i] != y[i]
                return false
            end
        end
        true
    end
end

# Comparison of PolyVar

function (==)(x::PolyVar{C}, y::PolyVar{C}) where C
    x.id == y.id
end

Base.isless(x::PolyVar{C}, y::PolyVar{C}) where C = isless(y.id, x.id)

# Comparison of Monomial

# graded lex ordering
function mycomp(x::Monomial{C}, y::Monomial{C}) where C
    degx = degree(x)
    degy = degree(y)
    if degx != degy
        degx - degy
    else
        i = j = 1
        # since they have the same degree,
        # if we get j > nvariables(y), the rest in x.z should be zeros
        while i <= nvariables(x) && j <= nvariables(y)
            if x.vars[i] > y.vars[j]
                if x.z[i] == 0
                    i += 1
                else
                    return 1
                end
            elseif x.vars[i] < y.vars[j]
                if y.z[j] == 0
                    j += 1
                else
                    return -1
                end
            elseif x.z[i] != y.z[j]
                return x.z[i] - y.z[j]
            else
                i += 1
                j += 1
            end
        end
        0
    end
end

function (==)(x::Monomial{C}, y::Monomial{C}) where C
    mycomp(x, y) == 0
end
(==)(x::PolyVar{C}, y::Monomial{C}) where C = convert(Monomial{C}, x) == y

# graded lex ordering
function Base.isless(x::Monomial{C}, y::Monomial{C}) where C
    mycomp(x, y) < 0
end
Base.isless(x::Monomial{C}, y::PolyVar{C}) where C = isless(x, convert(Monomial{C}, y))
Base.isless(x::PolyVar{C}, y::Monomial{C}) where C = isless(convert(Monomial{C}, x), y)

# Comparison of MonomialVector
function (==)(x::MonomialVector{C}, y::MonomialVector{C}) where C
    if length(x.Z) != length(y.Z)
        return false
    end
    allvars, maps = mergevars([_vars(x), _vars(y)])
    # Should be sorted in the same order since the non-common
    # polyvar should have exponent 0
    for (a, b) in zip(x.Z, y.Z)
        A = zeros(length(allvars))
        B = zeros(length(allvars))
        A[maps[1]] = a
        B[maps[2]] = b
        if A != B
            return false
        end
    end
    return true
end
(==)(mv::AbstractVector, x::MonomialVector) = monovec(mv) == x
(==)(x::MonomialVector, mv::AbstractVector) = x == monovec(mv)

# Comparison of Term
function (==)(p::Polynomial{C}, q::Polynomial{C}) where {C}
    # terms should be sorted and without zeros
    if length(p) != length(q)
        return false
    end
    for i in 1:length(p)
        if p.x[i] != q.x[i]
            # There should not be zero terms
            @assert p.a[i] != 0
            @assert q.a[i] != 0
            return false
        end
        if p.a[i] != q.a[i]
            return false
        end
    end
    return true
end

function grlex(x::Vector{Int}, y::Vector{Int})
    @assert length(x) == length(y)
    degx = sum(x)
    degy = sum(y)
    if degx != degy
        degx < degy
    else
        for (a, b) in zip(x, y)
            if a < b
                return true
            elseif a > b
                return false
            end
        end
        false
    end
end

function Base.isapprox(p::Polynomial{C, S}, q::Polynomial{C, T};
                       rtol::Real=Base.rtoldefault(S, T, 0), atol::Real=0,
                       ztol::Real=iszero(atol) ? Base.rtoldefault(S, T, 0) : atol) where {C, S, T}
    i = j = 1
    while i <= length(p.x) || j <= length(q.x)
        if i > length(p.x) || (j <= length(q.x) && q.x[j] > p.x[i])
            if !isapproxzero(q.a[j], ztol=ztol)
                return false
            end
            j += 1
        elseif j > length(q.x) || p.x[i] > q.x[j]
            if !isapproxzero(p.a[i], ztol=ztol)
                return false
            end
            i += 1
        else
            if !isapprox(p.a[i], q.a[j], rtol=rtol, atol=atol)
                return false
            end
            i += 1
            j += 1
        end
    end
    true
end


text:
#!/usr/bin/env julia6
using ReadGlobal

function makep!(d11::T,d22::T,d33::T,d12::T,d13::T,d23::T,
    w12::T,w13::T,w23::T,
    p11::T,p22::T,p33::T,p12::T,p13::T,p23::T) where {T<:Vector{<:Number}}

    Threads.@threads for i in 1:length(d11)
        @inbounds p11[i] = -2*(d12[i]*w12[i] + d13[i]*w13[i])
        @inbounds p22[i] = 2*(d12[i]*w12[i] - d23[i]*w23[i])
        @inbounds p33[i] = 2*(d13[i]*w13[i] + d23[i]*w23[i])
        @inbounds p12[i] = w12[i]*(d11[i]-d22[i]) - d13[i]*w23[i] - d23[i]*w13[i]
        @inbounds p13[i] = w13[i]*(d11[i]-d33[i]) + d12[i]*w23[i] - d23[i]*w12[i]
        @inbounds p23[i] = w23[i]*(d22[i]-d33[i]) + d12[i]*w13[i] + d13[i]*w12[i]
    end
    return 0
end

function makeP()
    N, Fil = getnfilter()

    d11 = "$(Fil)D11_N$N"
    d22 = "$(Fil)D22_N$N"
    d33 = "$(Fil)D33_N$N"
    d12 = "$(Fil)D12_N$N"
    d13 = "$(Fil)D13_N$N"
    d23 = "$(Fil)D23_N$N"

    w12 = "$(Fil)W12_N$N"
    w13 = "$(Fil)W13_N$N"
    w23 = "$(Fil)W23_N$N"

    p11 = "$(Fil)P11_N$N"
    p22 = "$(Fil)P22_N$N"
    p33 = "$(Fil)P33_N$N"
    p12 = "$(Fil)P12_N$N"
    p13 = "$(Fil)P13_N$N"
    p23 = "$(Fil)P23_N$N"

    doinchunks(makep!,input=(d11,d22,d33,d12,d13,d23,w12,w13,w23),output=(p11,p22,p33,p12,p13,p23))
end

makeP()


text:
struct ImageWidget <: Widget
    center::Point{2,Float64}
    dims::Point{2,Float64}
    uv_scale::Point{2,Float64}
end

Base.show(io::IO, w::ImageWidget) = print(io, "Image(width=", w.dims[1], ", height=", w.dims[2], ")")

vertex_data_type(::Type{ImageWidget}) = PosUV{Point2f,Point2f}
mesh_encoding_type(T::Type{ImageWidget}) = MeshVertexEncoding{TriangleStrip, vertex_data_type(T)}

function GeometryExperiments.PointSet(img::ImageWidget)
    set = PointSet(HyperCube, Point2f)
    (Translation(img.center) ∘ Scaling(img.dims) ∘ Scaling(1 / 2, 1 / 2))(set)
end

function GeometryExperiments.MeshVertexEncoding(img::ImageWidget)
    pos = (Translation(-1.0, -1.0) ∘ inv(Scaling(1920 / 2, 1080 / 2)))(PointSet(img))
    uv = (Scaling(img.uv_scale) ∘ Scaling(0.5, 0.5) ∘ Translation(1.0, 1.0))(PointSet(HyperCube, Point2f))
    MeshVertexEncoding(collect(vertex_data_type(img).(pos.points, uv.points)), Triangle)
end

nvertices(::Type{ImageWidget}) = 4
resource_types(::Type{ImageWidget}) = (SampledImage,)


text:
# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.5.2
#   kernelspec:
#     display_name: Julia 1.5.0
#     language: julia
#     name: julia-1.5
# ---

# # Clang
#
# clang is an open-source compiler built on the LLVM framework and targeting C, C++, and Objective-C. LLVM is also the JIT backend for Julia. Due to a highly modular design, Clang has in recent years become the core of a growing number of projects utilizing pieces of the compiler, such as tools for source-to-souce translation, static analysis and security evaluation, and editor tools for code completion, formatting, etc.
#
# While LLVM and Clang are written in C++, the Clang project maintains a C-exported interface called "libclang" which provides access to the abstract syntax tree and type representations. Thanks to the ubiquity of support for C calling conventions, a number of languages have utilized libclang as a basis for tooling related ot C and C++.
#
# The Clang.jl Julia package wraps libclang, provides a small convenience API for Julia-style programming, and provides a C-to-Julia wrapper generator built on libclang functionality.

using Clang
using Clang_jll

# +
example_h = raw"""
struct ExStruct {
    int    kind;
    char*  name;
    double* data;
};

void setvalue(struct ExStruct* e, int k, char* c, double* data);
struct ExStruct create(int k, char* c, double* data);
"""

example_c = raw"""
#include "example.h"

void setvalue(struct ExStruct* e, int k, char* c, double* data){
	e->kind = k;
	e->name = c;
	e->data = data;
}

struct ExStruct create(int k, char* c, double* data){
	struct ExStruct e;
	e.kind = k;
	e.name = c;
	e.data = data;
	return e;
}
"""

main_c = raw"""
#include <stdio.h>
#include "example.h"

int main(){
	struct ExStruct e;
	double a[3];
	a[0]=-0.0;
	a[1]=-1.0;
	a[2]=-2.0;
	e = create(42, "Hello", a);
	printf("%d\n", e.kind);
	printf("%s\n", e.name);
	printf("%f\n", e.data[0]);
	printf("%f\n", e.data[1]);
	printf("%f\n", e.data[2]);

	setvalue(&e, 42, "Hello", a);
	printf("%d\n", e.kind);
	printf("%s\n", e.name);
	printf("%f\n", e.data[0]);
	printf("%f\n", e.data[1]);
	printf("%f\n", e.data[2]);
	return 0;
}
"""

headerfile = "example.h"
cfile = "example.c"
mainfile = "main.c"

open(headerfile, "w") do f
    print(f, example_h)
end

open(cfile, "w") do f
    print(f, example_c)
end

open(mainfile, "w") do f
    print(f, main_c)
end

trans_unit = parse_header(headerfile)

# +

const LIBCLANG_INCLUDE = joinpath(dirname(Clang_jll.libclang_path), "..", "include", "clang-c") |> normpath
const TARGET_HEADERS = ["example.h"]

wrapcontext = init(;
    headers = TARGET_HEADERS,
    output_file = joinpath(@__DIR__, "example_api.jl"),
    common_file = joinpath(@__DIR__, "example_common.jl"),
    clang_includes = vcat(LIBCLANG_INCLUDE, CLANG_INCLUDE),
    clang_args = ["-I", joinpath(LIBCLANG_INCLUDE, "..")],
    header_wrapped = (root, current)->root == current,
    header_library = x->"example",
    clang_diagnostics = true,
)

wrapcontext |> run
# -

run(`gcc -shared -fPIC -o example.so example.c`)
run(`gcc main.c example.c`)
run(`./a.out`)

# +
const example="./example.so" # should be set const
include("example_common.jl")
include("example_api.jl")

e = create(42,"Hello",[1.,2.,3.])
refe = Ref(e)

@show e.kind
@show e.name |> unsafe_string
@show unsafe_load.(Ref(e.data), [1,2,3])

setvalue(refe,1,"Hi",[-1.,-2.,-3.])

@show refe.x.kind
@show refe.x.name |> unsafe_string
@show unsafe_load.(Ref(refe.x.data), [1,2,3])
# -

# # Clean up

# +
rm("example.c")
rm("main.c")
rm("example.h")
rm("example.so")
rm("a.out")

rm("example_common.jl")
rm("example_api.jl")

rm("LibTemplate.jl")
rm("clang.jl")
rm("ctypes.jl")


text:
function initializeAgent()
    m_s_t_min = s_0 # Set initial brain state prior
    v_s_t_min = 0.01

    # Initialize messages and marginals dictionary
    messages = Vector{Message}(undef, 59)
    marginals = Dict{Symbol, ProbabilityDistribution}()
    function infer(u_t::Float64, o_t::Float64)
        data = Dict(:o_t       => o_t,
                    :u_t       => u_t,
                    :m_o       => zeros(T-1),
                    :v_o       => sigma*ones(T-1),
                    :m_s_t_min => m_s_t_min,
                    :v_s_t_min => v_s_t_min)

        step!(data, marginals, messages)
    end

    marginals[:u_2] = ProbabilityDistribution(Univariate, GaussianMeanVariance, m=0.0, v=tiny) # Register initial action
    act() = mode(marginals[:u_2]) # Choose the mode of the current control state as action

    function slide(slide_msg_idx=7)
        (m_s_t_min, V_s_t_min) = ForneyLab.unsafeMeanCov(messages[slide_msg_idx].dist) # Reset prior state statistics
    end

    return (infer, act, slide)
end

max_stars_repo_name: Yikeru/BMLIP

max_issues_repo_name: Yikeru/BMLIP

max_forks_repo_name: Yikeru/BMLIP


text: Dict{String, Any}("f" => Dict{String, Any}("value" => "false", "type" => "bool"), "t" => Dict{String, Any}("value" => "true", "type" => "bool"))


text:
## (c) 2016 David A. van Leeuwen

## Unit tests for ../src/constructors.jl

## This code is licensed under the MIT license
## See the file LICENSE.md in this distribution

using DataStructures
using Test

print("construction, ")
## constructors
a = rand(2,3)
n1 = @inferred NamedArray(a, (OrderedDict(:a => 1, :b=>2), OrderedDict(:c=>3, :d=>2, :e=>1)), (:dim1, :dim2))
println(n1, " ", n1[1,1], " ", n1[:a, :e])
n2 = @inferred NamedArray(a, (OrderedDict(:a => 1, :b=>2), OrderedDict(:c=>3, :d=>2, :e=>1)))
println(n2)
n3 = @inferred NamedArray(a, ([:a, :b], [:c, :d, :e]), (:dim1, :dim2))
println(n3)
n4 = @inferred NamedArray(a, ([:a, :b], [:c, :d, :e]))
println(n4)
n5 = @inferred NamedArray(a, [[Char(64+i) for i in 1:d] for d in 2:3], ["a", "b"])
n6 = @inferred NamedArray(a, [[Char(64+i) for i in 1:d] for d in 2:3])
n7 = @inferred NamedArray(a)

n1 = @inferred NamedArray(ComplexF32, 5, 8)
n2 = @inferred NamedArray(a, (["s", "t"],[:a, :b, :c]), ("string", :symbol))

n = @inferred NamedArray(rand(2,4))
@inferred setnames!(n, ["one", "two"], 1)
@inferred setnames!(n, ["a", "b", "c", "d"], 2)

a = [1 2 3; 4 5 6]
n3 = @inferred NamedArray(a, (["a","b"],["C","D","E"]))
n4 = @inferred NamedArray(a, (OrderedDict("a"=>1,"b"=>2), OrderedDict("C"=>1,"D"=>2,"E"=>3)))

@test n3.array == n4.array == a
@test dimnames(n3) == dimnames(n4) == Any[:A,:B]
@test names(n3, 1) == names(n4, 1) == ["a","b"]
@test names(n3, 2) == names(n4, 2) == ["C","D","E"]

## 0-dim case #21
n0 = @inferred NamedArray(Array{Int}(undef))
@test size(n0) == ()
@test n0[1] == n0.array[1]

## Calling constructors through convert, #38
@test convert(NamedArray, n.array) == NamedArray(n.array)
@test convert(NamedVector, [1, 2, 3]) == NamedArray([1, 2, 3])

## new uninitialized construction syntax, #46
for i in 1:5
	dims = 2*fill(1, i)
	global n1 = @inferred NamedArray(Int, dims...) # implicit assignment to global variable `n1`
	global n2 = @inferred NamedArray{Int}(dims...) # implicit assignment to global variable `n2`
end

## repeated indices #63
@test_throws ErrorException @inferred n[["one", "one"], :]


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "XML2Builder"
version = v"2.9.9"

# Collection of sources required to build XML2Builder
sources = [
    "https://github.com/GNOME/libxml2/archive/v$(version).tar.gz" =>
    "d673f0284cec867ee00872a8152e0c3c09852f17fd9aa93f07579a37534f0bfe",
]

# Bash recipe for building across all platforms
script = raw"""
cd ${WORKSPACE}/srcdir/libxml2-*
./autogen.sh
./configure --prefix=${prefix} --host=${target} --without-python --with-zlib=${prefix}
make -j${nproc} install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

# The products that we will ensure are always built
products(prefix) = [
    LibraryProduct(prefix, "libxml2", :libxml2)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    "https://github.com/bicycle1885/ZlibBuilder/releases/download/v1.0.4/build_Zlib.v1.2.11.jl"
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)


text:
using OptimJulia
using SymPy

x1,x2 = Sym("x1,x2")

vars = [x1,x2]
f = -log(1 - x1 - x2) - log(x1) - log(x2)
x0 = [0.85,0.05]

Newton(f, vars, x0)


text:
using DynamicPricingExamples
using Test


text:
# TODO write tests with dummy structs and descriptors
# The code is currently tested in a PowerSystems unit test.


text:
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.8.0-DEV
#     language: julia
#     name: julia-1.8
# ---

# %% [markdown]
# https://discourse.julialang.org/t/function-compiles-every-time/63273

# %%
using StaticArrays
using AstrodynamicalModels
using DifferentialEquations

# %%
using StaticArrays
using AstrodynamicalModels
using DifferentialEquations

defaults = (; reltol=1e-14, abstol=1e-14, save_everystep=false)
kwargs   = (;)
options  = merge(defaults, kwargs)

r = randn(3) * 1e4
v = randn(3) * 1e3
T = 10_000.0

p1 = ODEProblem(R2BP, MVector{6}(vcat(r,v)), (0.0, T), MVector{1}(398600.435))
T1 = typeof(p1)

p2 = ODEProblem(R2BP, MVector{6}(vcat(r,v)), (0.0, T), MVector{1}(398600.435))
T2 = typeof(p2)

T1 === T2 #> false

# %%
using StaticArrays
using AstrodynamicalModels
using DifferentialEquations

function propagate(r, v, T; kwargs...)

  defaults = (; reltol=1e-14, abstol=1e-14, save_everystep=false)
  options  = merge(defaults, kwargs)

  problem   = ODEProblem(R2BP, MVector{6}(vcat(r,v)), (0.0, T), MVector{1}(398600.435))
  solutions = solve(problem; options...)

end

r = randn(3) * 1e4
v = randn(3) * 1e3
T = 10_000.0

@time propagate(r,v,T) #  28.067676 seconds (93.98 M allocations: 5.404 GiB, 3.61% gc time)
@time propagate(r,v,T) #   7.187065 seconds (13.05 M allocations: 704.718 MiB, 1.34% gc time, 99.70% compilation time)
@time propagate(r,v,T) #   7.154281 seconds (13.06 M allocations: 705.808 MiB, 1.79% gc time, 99.43% compilation time)

# %%
using StaticArrays
using AstrodynamicalModels
using DifferentialEquations

defaults = (; reltol=1e-14, abstol=1e-14, save_everystep=false)
kwargs   = (;)
options  = merge(defaults, kwargs)

r = randn(3) * 1e4
v = randn(3) * 1e3
T = 10_000.0

problem   = ODEProblem(R2BP, MVector{6}(vcat(r,v)), (0.0, T), MVector{1}(398600.435))

@time solutions = solve(problem; options...) #  15.728376 seconds (50.96 M allocations: 2.763 GiB, 4.17% gc time, 2.34% compilation time)
@time solutions = solve(problem; options...) #   0.012221 seconds (540.66 k allocations: 9.171 MiB)
@time solutions = solve(problem; options...) #   0.012453 seconds (540.66 k allocations: 9.171 MiB)

# %%
using StaticArrays
using AstrodynamicalModels
using DifferentialEquations

function propagate_remake(problem_default, r, v, T; kwargs...)

  defaults = (; reltol=1e-14, abstol=1e-14, save_everystep=false)
  options  = merge(defaults, kwargs)

  problem_remaked = remake(problem;
        u0 = MVector{6}(vcat(r,v)), tspan = (0.0, T), p = MVector{1}(398600.435))
  solutions = solve(problem_remaked; options...)

end

r = randn(3) * 1e4
v = randn(3) * 1e3
T = 10_000.0
problem = ODEProblem(R2BP, MVector{6}(vcat(r,v)), (0.0, T), MVector{1}(398600.435))

r = randn(3) * 1e4
v = randn(3) * 1e3
T = 10_000.0
@time sol1 = propagate_remake(problem, r, v, T)

r = randn(3) * 1e4
v = randn(3) * 1e3
T = 10_000.0
@time sol2 = propagate_remake(problem, r, v, T)

r = randn(3) * 1e4
v = randn(3) * 1e3
T = 10_000.0
@time sol3 = propagate_remake(problem, r, v, T);

# 9.130384 seconds (13.13 M allocations: 715.117 MiB, 1.53% gc time, 0.16% compilation time)
# 0.014155 seconds (452.63 k allocations: 7.765 MiB)
# 0.013558 seconds (448.20 k allocations: 7.689 MiB)

# %%
using Plots
PP = []
for sol in (sol1, sol2, sol3)
    t = range(sol.prob.tspan...; length=1000)
    y = hcat((sol.(t; idxs=i) for i in 1:3)...)
    P = plot(t, y; label="")
    push!(PP, P)
end
plot(PP...; size=(720, 150), layout=(1, 3))

# %%

max_stars_repo_path: 0006/Each time you create problem, the type of problem changes.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_name: genkuroki/public

max_forks_repo_name: genkuroki/public


text:
## Banded Cholesky decomposition

struct BandedCholesky{T} <: Factorization{T}
    data::BandedMatrix{T}  # symmetric banded matrix
end

# conversion
convert(::Type{BandedCholesky{T}}, B::BandedCholesky{S}) where {T<:Number, S<:Number} =
    BandedCholesky{T}(convert(BandedMatrix{T}, B.data))

# size of the parent array
@inline size(A::BandedCholesky) = size(A.data)
@inline size(A::BandedCholesky, k::Integer) = size(A.data, k)
@inline bandwidth(A::BandedCholesky, k) = bandwidth(A.data, k)

# Cholesky factorisation.
function cholesky!(A::Symmetric{T,<:BandedMatrix}) where {T<:Number}
    P = parent(A)
    pbtrf!('U', size(A, 1), bandwidth(A), bandeddata(parent(P)))
    BandedCholesky{T}(P)
end
cholesky(A::Symmetric{<:Any, <:BandedMatrix}) = cholesky!(copy(A))
cholesky(F::BandedCholesky) = F # no op

function getindex(F::BandedCholesky, d::Symbol)
    d == :U && return F.data # UpperTriangular(F.data)
    d == :L && return F.data' # LowerTriangular(F.data)
    throw(KeyError(d))
end

## Utilities

# check if matrix is square before solution of linear system or before converting
checksquare(A::BandedCholesky) = checksquare(A.data)

## Conversion/Promotion

# Converts A and b to the narrowest blas type
for typ in [BandedCholesky]
    @eval function _convert_to_blas_type(A::$typ, B::AbstractVecOrMat{T}) where {T<:Number}
        TS = _promote_to_blas_type(eltype(A), eltype(B))
        AA = convert($typ{TS}, A)
        BB = convert(Array{TS, ndims(B)}, B)
        AA, BB # one of these two might make a copy
    end
end


text:
using Docile, Docile.Interface, Lexicon, DataFrames

myfilter(x::Module; files = [""], categories = [:comment, :module, :function, :method, :type, :typealias, :macro, :global]) =
    filter(metadata(x), files = files, categories = categories)

myfilter(x::Metadata; files = [""], categories = [:comment, :module, :function, :method, :type, :typealias, :macro, :global]) =
    filter(x, files = files, categories = categories)

# Stuff from Lexicon.jl:
writeobj(any) = string(any)
writeobj(m::Method) = first(split(string(m), "("))
# from base/methodshow.jl
function url(m)
    line, file = m
    try
        d = dirname(file)
        u = Pkg.Git.readchomp(`config remote.origin.url`, dir=d)
        u = match(Pkg.Git.GITHUB_REGEX,u).captures[1]
        root = cd(d) do # dir=d confuses --show-toplevel, apparently
            Pkg.Git.readchomp(`rev-parse --show-toplevel`)
        end
        if beginswith(file, root)
            commit = Pkg.Git.readchomp(`rev-parse HEAD`, dir=d)
            return "https://github.com/$u/tree/$commit/"*file[length(root)+2:end]*"#L$line"
        else
            return Base.fileurl(file)
        end
    catch
        return Base.fileurl(file)
    end
end

function mysave(file::String, m::Module, order = [:source])
    mysave(file, documentation(m), order)
end
function mysave(file::String, docs::Metadata, order = [:source])
    isfile(file) || mkpath(dirname(file))
    open(file, "w") do io
        info("writing documentation to $(file)")
        println(io)
        for (k,v) in EachEntry(docs, order = order)
            name = writeobj(k)
            source = v.data[:source]
            catgory = category(v)
            comment = catgory == :comment
            println(io)
            println(io)
            !comment && println(io, "## $name")
            println(io)
            println(io, v.docs.data)
            path = last(split(source[2], r"v[\d\.]+(/|\\)"))
            !comment && println(io, "[$(path):$(source[1])]($(url(source)))")
            println(io)
        end
    end
end

mysave("api/maintypes.md", myfilter(DataFrames, files = ["abstractdataframe.jl", "dataframe.jl", "subdataframe.jl"],
                                    categories = [:comment, :type, :typealias]))
mysave("api/utilities.md", myfilter(DataFrames, files = ["abstractdataframe.jl", "dataframe.jl", "subdataframe.jl"],
                                    categories = [:comment, :function, :method, :macro, :global]))
mysave("api/manipulation.md", myfilter(DataFrames, files = ["reshape.jl", "join.jl"]))


text:
using Documenter, StokesParticles

push!(LOAD_PATH,"../src/")

makedocs(;
    modules=[StokesParticles],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/pseastham/StokesParticles.jl/blob/{commit}{path}#L{line}",
    sitename="StokesParticles.jl",
    authors="Patrick Eastham",
    assets=String[],
)

deploydocs(;
    repo="github.com/pseastham/StokesParticles.jl",
)


text:
module SlurmClusterManager

export SlurmManager, launch, manage

using Distributed
import Distributed: launch, manage

include("slurmmanager.jl")

end # module


text:
using StaticArrays
using ClimateMachine.VariableTemplates
using ClimateMachine.BalanceLaws:
    BalanceLaw,
    Prognostic,
    Auxiliary,
    Gradient,
    GradientFlux,
    GradientLaplacian,
    Hyperdiffusive

import ClimateMachine.BalanceLaws:
    vars_state,
    number_states,
    flux_first_order!,
    flux_second_order!,
    source!,
    compute_gradient_argument!,
    compute_gradient_flux!,
    nodal_init_state_auxiliary!,
    update_auxiliary_state!,
    init_state_prognostic!,
    boundary_conditions,
    boundary_state!,
    wavespeed,
    transform_post_gradient_laplacian!

using ClimateMachine.Mesh.Geometry: LocalGeometry
using ClimateMachine.DGMethods.NumericalFluxes:
    NumericalFluxFirstOrder,
    NumericalFluxSecondOrder,
    NumericalFluxGradient,
    GradNumericalFlux,
    DivNumericalPenalty

import ClimateMachine.DGMethods.NumericalFluxes:
    numerical_flux_first_order!, boundary_flux_second_order!

using CLIMAParameters
struct EarthParameterSet <: AbstractEarthParameterSet end
const param_set = EarthParameterSet()

struct Advection{N} <: BalanceLaw end
struct NoAdvection <: BalanceLaw end

struct Diffusion{N} <: BalanceLaw end
struct NoDiffusion <: BalanceLaw end

struct HyperDiffusion{N} <: BalanceLaw end
struct NoHyperDiffusion <: BalanceLaw end

abstract type AdvectionDiffusionProblem end

"""
    AdvectionDiffusion{N} <: BalanceLaw

A balance law describing a system of `N` advection-diffusion-hyperdiffusion
equations:

```
∂ρ
-- = - ∇ • (u ρ - σ + η)
∂t

σ = D ∇ ρ
η = H ∇ Δρ
```
Where

 - `ρ` is the solution vector
 - `u` is the advection velocity
 - `σ` is the DG diffusion auxiliary variable
 - `D` is the diffusion tensor
 - `η` is the DG hyperdiffusion auxiliary variable
 - `H` is the hyperdiffusion tensor
"""
struct AdvectionDiffusion{N, dim, P, fluxBC, A, D, HD} <: BalanceLaw
    problem::P
    advection::A
    diffusion::D
    hyperdiffusion::HD

    function AdvectionDiffusion{dim}(
        problem::P;
        num_equations = 1,
        flux_bc = false,
        advection::Bool = true,
        diffusion::Bool = true,
        hyperdiffusion::Bool = false,
    ) where {dim, P <: AdvectionDiffusionProblem}
        N = num_equations
        adv = advection ? Advection{N}() : NoAdvection()
        A = typeof(adv)
        diff = diffusion ? Diffusion{N}() : NoDiffusion()
        D = typeof(diff)
        hyperdiff = hyperdiffusion ? HyperDiffusion{N}() : NoHyperDiffusion()
        HD = typeof(hyperdiff)
        new{N, dim, P, flux_bc, A, D, HD}(problem, adv, diff, hyperdiff)
    end
end

# Auxiliary variables, always store
# `coord` coordinate points (needed for BCs)
function vars_state(m::AdvectionDiffusion, st::Auxiliary, FT)
    @vars begin
        coord::SVector{3, FT}
        advection::vars_state(m.advection, st, FT)
        diffusion::vars_state(m.diffusion, st, FT)
        hyperdiffusion::vars_state(m.hyperdiffusion, st, FT)
    end
end

#   `u` advection velocity
vars_state(::Advection{1}, ::Auxiliary, FT) = @vars(u::SVector{3, FT})
vars_state(::Advection{N}, ::Auxiliary, FT) where {N} =
    @vars(u::SMatrix{3, N, FT, 3N})
#   `D` diffusion tensor
vars_state(::Diffusion{1}, ::Auxiliary, FT) = @vars(D::SMatrix{3, 3, FT, 9})
vars_state(::Diffusion{N}, ::Auxiliary, FT) where {N} =
    @vars(D::SArray{Tuple{3, 3, N}, FT, 3, 9N})
#   `H` hyperdiffusion tensor
vars_state(::HyperDiffusion{1}, ::Auxiliary, FT) =
    @vars(H::SMatrix{3, 3, FT, 9})
vars_state(::HyperDiffusion{N}, ::Auxiliary, FT) where {N} =
    @vars(H::SArray{Tuple{3, 3, N}, FT, 3, 9N})

# Density `ρ` is the only state
vars_state(::AdvectionDiffusion{1}, ::Prognostic, FT) = @vars(ρ::FT)
vars_state(::AdvectionDiffusion{N}, ::Prognostic, FT) where {N} =
    @vars(ρ::SVector{N, FT})

function vars_state(m::AdvectionDiffusion{N}, ::Gradient, FT) where {N}
    # For pure advection we don't need gradients
    if m.diffusion isa NoDiffusion && m.hyperdiffusion isa NoHyperDiffusion
        return @vars()
    else  # Take the gradient of density
        return N == 1 ? @vars(ρ::FT) : @vars(ρ::SVector{N, FT})
    end
end

# Take the gradient of laplacian of density ρ
vars_state(::HyperDiffusion{1}, ::GradientLaplacian, FT) = @vars(ρ::FT)
vars_state(::HyperDiffusion{N}, ::GradientLaplacian, FT) where {N} =
    @vars(ρ::SVector{N, FT})
vars_state(m::AdvectionDiffusion, st::GradientLaplacian, FT) =
    vars_state(m.hyperdiffusion, st, FT)

# The DG diffusion auxiliary variable: σ = D ∇ρ
vars_state(::Diffusion{1}, ::GradientFlux, FT) = @vars(σ::SVector{3, FT})
vars_state(::Diffusion{N}, ::GradientFlux, FT) where {N} =
    @vars(σ::SMatrix{3, N, FT, 3N})
vars_state(m::AdvectionDiffusion, st::GradientFlux, FT) =
    vars_state(m.diffusion, st, FT)

# The DG hyperdiffusion auxiliary variable: η = H ∇ Δρ
vars_state(::HyperDiffusion{1}, ::Hyperdiffusive, FT) = @vars(η::SVector{3, FT})
vars_state(::HyperDiffusion{N}, ::Hyperdiffusive, FT) where {N} =
    @vars(η::SMatrix{3, N, FT, 3N})
vars_state(m::AdvectionDiffusion, st::Hyperdiffusive, FT) =
    vars_state(m.hyperdiffusion, st, FT)

"""
    flux_first_order!(::Advection, flux::Grad, state::Vars, aux::Vars)

Computes non-diffusive flux `F_adv = u ρ` where

 - `u` is the advection velocity
 - `ρ` is the advected quantity
"""
function flux_first_order!(
    ::Advection{N},
    flux::Grad,
    state::Vars,
    aux::Vars,
) where {N}
    ρ = state.ρ
    u = aux.advection.u
    flux.ρ += u .* ρ'
end
flux_first_order!(::NoAdvection, flux::Grad, state::Vars, aux::Vars) = nothing
flux_first_order!(
    m::AdvectionDiffusion,
    flux::Grad,
    state::Vars,
    aux::Vars,
    t::Real,
    direction,
) = flux_first_order!(m.advection, flux, state, aux)

"""
    flux_second_order!(::Diffusion, flux::Grad, auxDG::Vars)

Computes diffusive flux `F_diff = -σ` where:

 - `σ` is DG diffusion auxiliary variable (`σ = D ∇ ρ`
    with `D` being the diffusion tensor)
"""
function flux_second_order!(::Diffusion, flux::Grad, auxDG::Vars)
    σ = auxDG.σ
    flux.ρ += -σ
end
flux_second_order!(::NoDiffusion, flux::Grad, auxDG::Vars) = nothing

"""
    flux_second_order!(::HyperDiffusion, flux::Grad, auxHDG::Vars)

Computes hyperdiffusive flux `F_hyperdiff = η` where:

 - `η` is DG hyperdiffusion auxiliary variable (`η = H ∇ Δρ`
    with `H` being the hyperdiffusion tensor)
"""
function flux_second_order!(::HyperDiffusion, flux::Grad, auxHDG::Vars)
    η = auxHDG.η
    flux.ρ += η
end
flux_second_order!(::NoHyperDiffusion, flux::Grad, auxHDG::Vars) = nothing

function flux_second_order!(
    m::AdvectionDiffusion,
    flux::Grad,
    state::Vars,
    auxDG::Vars,
    auxHDG::Vars,
    aux::Vars,
    t::Real,
)
    flux_second_order!(m.diffusion, flux, auxDG)
    flux_second_order!(m.hyperdiffusion, flux, auxHDG)
end

"""
    compute_gradient_argument!(m::AdvectionDiffusion, transform::Vars, state::Vars,
                   aux::Vars, t::Real)

Set the variable to take the gradient of (`ρ` in this case)
"""
function compute_gradient_argument!(
    m::AdvectionDiffusion,
    transform::Vars,
    state::Vars,
    aux::Vars,
    t::Real,
)
    transform.ρ = state.ρ
end

"""
    compute_gradient_flux!(::Diffusion, auxDG::Vars, gradvars::Grad, aux::Vars)

Computes the DG diffusion auxiliary variable `σ = D ∇ ρ` where `D` is
the diffusion tensor.
"""
function compute_gradient_flux!(
    ::Diffusion{N},
    auxDG::Vars,
    gradvars::Grad,
    aux::Vars,
) where {N}
    ∇ρ = gradvars.ρ
    D = aux.diffusion.D
    if N == 1
        auxDG.σ = D * ∇ρ
    else
        auxDG.σ = hcat(ntuple(n -> D[:, :, n] * ∇ρ[:, n], Val(N))...)
    end
end
compute_gradient_flux!(::NoDiffusion, auxDG::Vars, gradvars::Grad, aux::Vars) =
    nothing
compute_gradient_flux!(
    m::AdvectionDiffusion,
    auxDG::Vars,
    gradvars::Grad,
    state::Vars,
    aux::Vars,
    t::Real,
) = compute_gradient_flux!(m.diffusion, auxDG, gradvars, aux)

"""
    transform_post_gradient_laplacian!(::AdvectionDiffusion, auxHDG::Vars,
        gradvars::Grad, state::Vars, aux::Vars, t::Real)

Computes the DG hyperdiffusion auxiliary variable `η = H ∇ Δρ` where `H` is
the hyperdiffusion tensor.
"""
function transform_post_gradient_laplacian!(
    m::AdvectionDiffusion{N},
    auxHDG::Vars,
    gradvars::Grad,
    state::Vars,
    aux::Vars,
    t::Real,
) where {N}
    ∇Δρ = gradvars.ρ
    H = aux.hyperdiffusion.H
    if N == 1
        auxHDG.η = H * ∇Δρ
    else
        auxHDG.η = hcat(ntuple(n -> H[:, :, n] * ∇Δρ[:, n], Val(N))...)
    end
end

"""
    source!(m::AdvectionDiffusion, _...)

There is no source in the advection-diffusion model
"""
source!(m::AdvectionDiffusion, _...) = nothing

"""
    wavespeed(m::AdvectionDiffusion, nM, state::Vars, aux::Vars, t::Real)

Wavespeed with respect to vector `nM`
"""
wavespeed(
    m::AdvectionDiffusion,
    nM,
    state::Vars,
    aux::Vars,
    t::Real,
    direction,
) = wavespeed(m.advection, nM, aux)
function wavespeed(::Advection{N}, nM, aux::Vars) where {N}
    u = aux.advection.u
    if N == 1
        abs(nM' * u)
    else
        SVector(ntuple(n -> abs(nM' * u[:, n]), Val(N)))
    end
end
wavespeed(::NoAdvection, nM, aux::Vars) = 0

function nodal_init_state_auxiliary!(
    m::AdvectionDiffusion,
    aux::Vars,
    tmp::Vars,
    geom::LocalGeometry,
)
    aux.coord = geom.coord
    init_velocity_diffusion!(m.problem, aux, geom)
end

has_variable_coefficients(::AdvectionDiffusionProblem) = false
function update_auxiliary_state!(
    dg::DGModel,
    m::AdvectionDiffusion,
    Q::MPIStateArray,
    t::Real,
    elems::UnitRange,
)
    if has_variable_coefficients(m.problem)
        update_auxiliary_state!(dg, m, Q, t, elems) do m, state, aux, t
            update_velocity_diffusion!(m.problem, m, state, aux, t)
        end
        return true
    end
    return false
end

function init_state_prognostic!(
    m::AdvectionDiffusion,
    state::Vars,
    aux::Vars,
    localgeo,
    t::Real,
)
    initial_condition!(m.problem, state, aux, localgeo, t)
end

Neumann_data!(problem, ∇state, aux, x, t) = nothing
Dirichlet_data!(problem, state, aux, x, t) = nothing
boundary_conditions(::AdvectionDiffusion) = (1, 2, 3, 4)
function boundary_state!(
    nf,
    bctype,
    m::AdvectionDiffusion{N},
    stateP::Vars,
    auxP::Vars,
    nM,
    stateM::Vars,
    auxM::Vars,
    t,
    _...,
) where {N}
    if bctype == 1 # Dirichlet
        Dirichlet_data!(m.problem, stateP, auxP, (coord = auxP.coord,), t)
    elseif bctype ∈ (2, 4) # Neumann
        stateP.ρ = stateM.ρ
    elseif bctype == 3 # zero Dirichlet
        stateP.ρ = N == 1 ? 0 : zeros(typeof(stateP.ρ))
    end
end

function boundary_state!(
    nf::CentralNumericalFluxSecondOrder,
    bctype,
    m::AdvectionDiffusion,
    state⁺::Vars,
    diff⁺::Vars,
    aux⁺::Vars,
    n⁻::SVector,
    state⁻::Vars,
    diff⁻::Vars,
    aux⁻::Vars,
    t,
    _...,
)
    if m.diffusion isa NoDiffusion && m.hyperdiffusion isa NoHyperDiffusion
        return nothing
    end

    if bctype ∈ (1, 3) # Dirchlet
        # Just use the minus side values since Dirchlet
        diff⁺.σ = diff⁻.σ
    elseif bctype == 2 # Neumann with data
        FT = eltype(diff⁺)
        ngrad = number_states(m, Gradient())
        ∇state = Grad{vars_state(m, Gradient(), FT)}(similar(
            parent(diff⁺),
            Size(3, ngrad),
        ))
        # Get analytic gradient
        Neumann_data!(m.problem, ∇state, aux⁻, aux⁻.coord, t)
        compute_gradient_flux!(m.diffusion, diff⁺, ∇state, aux⁻)
        # compute the diffusive flux using the boundary state
    elseif bctype == 4 # zero Neumann
        FT = eltype(diff⁺)
        ngrad = number_states(m, Gradient())
        ∇state = Grad{vars_state(m, Gradient(), FT)}(similar(
            parent(diff⁺),
            Size(3, ngrad),
        ))
        # Get analytic gradient
        ∇state.ρ = zeros(typeof(∇state.ρ))
        # convert to auxDG variables
        compute_gradient_flux!(m.diffusion, diff⁺, ∇state, aux⁻)
    end
    nothing
end

function boundary_flux_second_order!(
    nf::CentralNumericalFluxSecondOrder,
    bctype,
    m::AdvectionDiffusion{N, dim, P, true},
    F,
    state⁺,
    diff⁺,
    hyperdiff⁺,
    aux⁺,
    n⁻,
    state⁻,
    diff⁻,
    hyperdiff⁻,
    aux⁻,
    t,
    _...,
) where {N, dim, P}
    if m.diffusion isa NoDiffusion && m.hyperdiffusion isa NoHyperDiffusion
        return nothing
    end

    # Default initialize flux to minus side
    if bctype ∈ (1, 3) # Dirchlet
        # Just use the minus side values since Dirchlet
        flux_second_order!(m, F, state⁻, diff⁻, hyperdiff⁻, aux⁻, t)
    elseif bctype == 2 # Neumann data
        FT = eltype(diff⁺)
        ngrad = number_states(m, Gradient())
        ∇state = Grad{vars_state(m, Gradient(), FT)}(similar(
            parent(diff⁺),
            Size(3, ngrad),
        ))
        # Get analytic gradient
        Neumann_data!(m.problem, ∇state, aux⁻, aux⁻.coord, t)
        # get the diffusion coefficient
        D = aux⁻.diffusion.D
        # exact the exact data
        ∇ρ = ∇state.ρ
        # set the flux
        if N == 1
            F.ρ = -D * ∇ρ
        else
            F.ρ = hcat(ntuple(n -> -D[:, :, n] * ∇ρ[:, n], Val(N))...)
        end
    elseif bctype == 4 # Zero Neumann
        F.ρ = zeros(typeof(F.ρ))
    end
    nothing
end

# Bcs for hyperdiffusion not implemented
boundary_state!(
    ::Union{GradNumericalFlux, DivNumericalPenalty},
    bctype,
    ::AdvectionDiffusion,
    _...,
) = nothing

max_forks_repo_name: jm-c/CLIMA


text:
"""
type `SimUnit` has fields
* x::Vehicle
* f::AbstractFilter
* p::Policy
* cm::CostModel

Constructor: `SimUnit(x, f, p, cm)`
"""
struct SimUnit
    x::Vehicle
    f::AbstractFilter
    p::Policy
    cm::CostModel
    tc::TerminationCondition
end
function SimUnit(x::Vehicle,
                 f::AbstractFilter,
                 p::Policy,
                 tc::TerminationCondition = StepThreshold(10),
                 cm::CostModel = ConstantCost(1.0)
                )
    return SimUnit(x, f, p, cm, tc)
end

get_cost(su::SimUnit, m) = get_cost(su.cm, m, su.x, su.f)
get_cost(su::SimUnit, m, a::Action) = get_cost(su.cm, m, su.x, su.f, a)

function reset!(m::SearchDomain, su::SimUnit)
    reset!(su.f)
    reset!(m, su.x)
    reset!(su.p)
end

# The following are just helper functions that simplify signatures
update!(su::SimUnit, o) = update!(su.f, su.x, o)
action(m::SearchDomain, su::SimUnit, o) = action(m, su.x, o, su.f, su.p)

function is_complete(vsu::Vector{SimUnit}, step_count::Int)
    ret_val = false
    for uav in vsu
        if is_complete(uav.f, uav.tc, step_count)
            ret_val = true
        end
    end
    return ret_val
end


text:
const GERG2008_consts = (
names =[
        "methane",
        "nitrogen",
        "carbon dioxide",
        "ethane",
        "propane",
        "butane",
        "isobutane",
        "pentane",
        "isopentane",
        "hexane",
        "heptane",
        "octane",
        "nonane",
        "decane",
        "hydrogen",
        "oxygen",
        "carbon monoxide",
        "water",
        "hydrogen sulfide",
        "helium",
        "argon"
        ],
        #calculated from Tmin and Pmax
lb_volumes = [
    2.952789230664869e-5
    2.595917396536852e-5
    3.0350520997504374e-5
    4.181265338895899e-5
    5.7773265228561896e-5
    7.778071972263543e-5
    7.529511394064539e-5
    9.03473213517436e-5
    9.138272716465976e-5
    0.00010746846161549492
    0.00012307769607878738
    0.00014149360905504386
    0.00013613745434844812
    0.00015142994339404644
    1.6005013101722185e-5
    2.3229544237311946e-5
    2.9381748144372254e-5
    1.732307538792257e-5
    3.154736660208043e-5
    1.2820394828598075e-5
    2.511877404958246e-5
],
acentricfactor = [
        0.0115478,
        0.0377215,
        0.223621,
        0.099493,
        0.152291,
        0.200164,
        0.183521,
        0.251506,
        0.227875,
        0.301261,
        0.349469,
        0.399552,
        0.44346,
        0.492328,
        -0.215993,
        0.0221798,
        0.0481621,
        0.344861,
        0.0941677,
        -0.390032,
        0.0,
    ],
    Mw = [
        16.04246,
        28.0134,
        44.0095,
        30.06904,
        44.09562,
        58.1222,
        58.1222,
        72.14878,
        72.14878,
        86.17536,
        100.20194,
        114.22852,
        128.2551,
        142.28168,
        2.01588,
        31.9988,
        28.0101,
        18.01528,
        34.08088,
        4.002602,
        39.948,
    ],
    Tc = [
        190.564,
        126.192,
        304.1282,
        305.322,
        369.825,
        407.817,
        425.125,
        460.35,
        469.7,
        507.82,
        540.13,
        569.32,
        594.55,
        617.7,
        33.19,
        154.595,
        132.86,
        647.096,
        373.1,
        5.1953,
        150.687,
    ],
    rhoc = [
        10.139342719,
        11.1839,
        10.624978698,
        6.87085454,
        5.000043088,
        3.86014294,
        3.920016792,
        3.271,
        3.215577588,
        2.705877875,
        2.315324434,
        2.056404127,
        1.81,
        1.64,
        14.94,
        13.63,
        10.85,
        17.87371609,
        10.19,
        17.399,
        13.407429659,
    ], #molar density, in mol/dm3
    pc = [
        4599000.,
        3400000.,
        7383000.,
        4872000.,
        4248000.,
        3796000.,
        3640000.,
        3370000.,
        3380000.,
        3025000.,
        2740000.,
        2490000.,
        2290000.,
        2110000.,
        1313000.,
        5043000.,
        3499000.,
        22064000.,
        8962910.,
        117000.,
        4898000.,
    ]#critical pressure, Pa
)

function GERG2008_xsel(names::Vector{String})

symvalues = GERG2008_consts.names
    xsel = Vector{Int}()
    for i in names
        i in symvalues || throw(ArgumentError("invalid name: $i"))
        push!(xsel, findall(isequal(i), symvalues)[1])

    end
    if iszero(length(xsel))
        return ArgumentError("no match found for input components. valid components are: $(GERG2008constsProps.names)")
    else
        return xsel
    end
end


text:
include("CFRcs.jl")

function rotn!(hist, n)
    if length(hist) > n
        popfirst!(hist)
    end
end

function makebids(m, n, r)
    vcat([100i + 10j for j = 0:(r-1), i = m:n]...)
end

function moves(bids, history)
    n = length(history)
    if n > 2 && history[end-2] % 10 == 5 && history[end] % 10 == 5
        return [1]
    elseif n > 1 && history[end] % 10 == 5
        b = findfirst(isequal(history[end-1]), bids) + 1
        return pushfirst!(bids[b:end], 1)
    elseif n > 0
        b = findfirst(isequal(history[end]), bids) + 1
        return pushfirst!(bids[b:end], history[end] + 5)
    else
        return bids
    end
end

function mkGraph(q, r, m)
    bids = makebids(1, 2q, r)
    initBids = map(x -> [x], bids)
    graph = Dict([] => initBids)
    parent = []
    children = Set(initBids)
    while !isempty(children)
        vertex = pop!(children)
        if vertex[end] == 1
            graph[vertex] = []
            continue
        end
        acts = moves(bids, vertex)
        for a in acts
            v = push!(copy(vertex), a)
           rotn!(v, m)
            if !haskey(graph, v)
                push!(children, v)
            end
            if haskey(graph, vertex)
                graph[vertex] = push!(copy(graph[vertex]), v)
            else
                graph[vertex] = [v]
            end
        end
    end
    graph
end

function hasincoming(graph, v)
    vals = values(graph)
    for xs in vals
        if v ∈ xs
            return true
        end
    end
    false
end

function kahn(graph)
    graph = deepcopy(graph)
    result = []
    s = Set([[]])
    while !isempty(s)
        # v = pop!(s)
        v = minimum(s)
        delete!(s, v)
        push!(result, v)
        while !isempty(graph[v])
            u = pop!(graph[v])
            if !hasincoming(graph, u)
                push!(s, u)
            end
        end
    end
    result
end

max_stars_repo_name: jeffreyrosenbluth/LP

max_issues_repo_name: jeffreyrosenbluth/LP

max_forks_repo_name: jeffreyrosenbluth/LP


text:
include("../src/multivector.jl")
@basis e1 e2 e3
# Vectors
a = e1
b = e2
c = e3

# Bivectors
A = e1*e2
B = e2*e3
C = e3*e1

# Pseudoscalar
I = e1*e2*e3

# Duals
i = -1*I*A
j = -1*I*B
k = -1*I*C