text:
### A Pluto.jl notebook ###
# v0.15.1

using Markdown
using InteractiveUtils

# ╔═╡ ab976d70-f59d-4dcd-ac89-771606a74108
begin
	using PlutoUI
	using PlutoTest
	using PlutoTeachingTools
	#eval(Meta.parse(code_for_check_type_funcs))
end

# ╔═╡ 5b826e17-ad51-43d0-9756-89a9fd0ea834
md"""
# Astro 528, Lab 1, Exercise 1
"""

# ╔═╡ acb10cb1-d03d-47c5-90dc-7216e5110ed2
md"""
## Access the lab
Before we delve into the meat of the course, we'll practice accessing the assignment via your private git repository ("repo"), cloning the repo, starting Julia, and launching the Pluto notebook server.
Once you open the Pluto notebook, it will begin executing the entire notebook.  (Depending on your system and how complex the notebook is, sometimes this can take a little time.)  You'll read through this notebook, make some small changes, see how that affects the results, save your changes to your local machine, commit your changes to your local git repository, and then push those to your repository on github.
For many of you, this is likely your first time using all of these tools, and it probably sounds like a lot.  The good news is that before long all of this will become second nature and you'll wonder how you got along without such tools.
"""

# ╔═╡ c15ce452-494f-407c-9c08-02d547929e6a
md"""
## Assigning Variables
In Julia, you assign variables values using the `=` sign.
Assigning integer or floating point numbers uses a natural notation.
Once a cell is executed, Pluto shows the output of each cell above the cell with the relevant code.
For example,"""

# ╔═╡ 25d5497c-3a67-4495-9f37-0bc8f2a64e68
current_year = 2021

# ╔═╡ 774b2e4a-6cd7-4c9a-920d-b4f38354d42d
md"""Strings are created by enclosing text between double quotes.
For example,"""

# ╔═╡ 5917df75-0dab-4778-917e-91bbfddbe03a
begin
	course_number = "Astro 528"
	semester = "Fall 2021"
	instructor = "Eric Ford"
end

# ╔═╡ 0d15ba33-f2c4-46d6-9b08-9181d40596cd
md"""Notice that, when working in a Pluto notebook, cells that assign multiple variables need to be wraped in a `begin...end` block.  (This is not necessary in general (e.g., julia scripts, Jupyter notebooks, or module code.)
Note that only the output from the last line is displayed.  The other variables were still assigned and can be accessed in other cells."""

# ╔═╡ dc0bf251-b3a6-42f7-84f6-796ac3761b2d
course_number

# ╔═╡ dbfceab2-994d-45a1-97ff-d978827bd129
tip(md"""Single quotes are used to indicate a single character (as opposed to a string of length 1).
[Triple quoted string literals](https://docs.julialang.org/en/v1/manual/strings/#Triple-Quoted-String-Literals) are useful for quotes spanning multiple lines (and when you want to include `"` symbols).
""")

# ╔═╡ e2bb48e5-a2b2-42d7-a8e5-5c7532d62e93
md"Assigning integer or floating point numbers uses a natural notation."

# ╔═╡ cda49166-b9fb-4e68-8b3e-b97845c7bf3b
instructor_year = current_year - 1999

# ╔═╡ 0c5d5b14-6c9a-4ac9-8126-0d6e1edba593
favorite_number = 3.1415926535897

# ╔═╡ a3450384-ccfc-4bda-af52-24ed8f429cf8
md"""Now, it's your turn to create variables named `student_name`, `student dept` and `student_year` with values based on yourself.  You can modify the code in the cells below."""

# ╔═╡ 33350a21-7971-406e-b316-6539eca72bf4
begin  # INSERT YOUR RESPONSES BELOW
	student_name = missing
	student_dept = missing
	student_year = missing
end;

# ╔═╡ f15643a4-39d5-4492-98d3-9db4f4cb0d3b
md"""
Run the cell with your code and make sure there are no error messages.
Sometimes, I may provide a few quick checks in the Pluto notebook itself.  For example,
"""

# ╔═╡ 1bc2edaf-a0a7-4c4e-96c2-25872dd30e2e
if !@isdefined(student_year)
	not_defined(:student_year)
else
	let
		if ismissing(student_name)
			still_missing()
		elseif isnothing(student_name)
			keep_working("Could there be a typeo?")
		elseif !(student_name isa String)
			keep_working(md"The value should be a String.")
		else
			correct(md"You set `student_name` to $student_name.")
		end
	end
end

# ╔═╡ c05d62b5-6063-41ca-9823-2499b2d4fbc4
if !@isdefined(student_dept)
	not_defined(:student_dept)
else
	let
		if ismissing(student_dept)
			still_missing()
		elseif isnothing(student_dept)
			keep_working("Could there be a typeo?")
		elseif !(student_dept isa String)
			keep_working(md"The value of `student_year` should be a String.")
		elseif occursin("Astro", student_dept)
			correct(md"Welcome back.")
		else
			correct(md"We're glad to have you join us.")
		end
	end
end

# ╔═╡ 8b666639-e3d0-4535-803d-40ddc4360a9c
if !@isdefined(student_year)
	not_defined(:student_year)
else
	let
		if ismissing(student_year)
			still_missing()
		elseif isnothing(student_year)
			keep_working("Could there be a typeo?")
		elseif !(student_year isa Integer)
			keep_working(md"The value should be an integer.")
		elseif ! (1<= student_year)
			keep_working(md"I'm expecting `student_year` to be one or larger.")
		elseif ! (student_year <= 10)
			keep_working(md"Are you sure?  Most grad students finish in less than 10 years.")
		else
			correct(md"I think that $student_year seems like a reasonable value for `student_year`.")
		end
	end
end

# ╔═╡ ec4ec6e7-9548-4bd6-970a-48fefa9861d7
md"""
## Pass a Test
Once we start doing real calculations, it will be important to test our work.  Therefore, you and I will both write tests that your code should pass. For example, in this case we can test that variables with the correct names exist, they're not empty, that they have the expected type, that their value falls within a reasonable range, etc.
"""

# ╔═╡ ac9e6787-5c42-44d6-907e-e7033e7ab4bc
tip(md"Before using a test, we need to tell Julia to import a module that will provide the `@test` macro.  To simplify the notebook, I've put that code at the bottom of the notebook, so it doesn't distract from our main task.")

# ╔═╡ 087b40c4-7d59-4f9a-b584-65c2aa9f9578
md"
The first few tests below should pass.  Later, I've written a couple of tests with bugs in them for you to fix.  Look at the results of the cells below, find the tests that doesn't pass, and modify the code so that it passes.
"

# ╔═╡ eefa92a1-1d4c-4536-abd3-1c1af98e146f
md"One of the first things to check is that variables have been defined.  If not, then later tests will fail and we might spend time trying to dive deeper into the problem than we need to."

# ╔═╡ 78f81544-fead-4649-b306-48fb1463d397
@test @isdefined(student_name)

# ╔═╡ 3d224633-0eb8-407c-b44c-0d9ec927dfc6
@test @isdefined(student_dept)

# ╔═╡ 2fe3f4e8-5288-4f61-9346-7a455c2abc85
@test @isdefined(student_year)

# ╔═╡ b49789b2-49c0-49e3-b81f-db04e60f9d23
md"Next, we'll test that the types of the variables are what we expect.  If not, then it's unlikely that the values stored will be what we expect."

# ╔═╡ 34f767c9-f95f-4968-9423-ad87bb456775
@test typeof(student_name) == String

# ╔═╡ d6417bdc-9e36-4d2b-8b3e-f5ddafd82437
@test typeof(student_dept) == String

# ╔═╡ 89ac58ea-b782-4dcc-802f-d6d063e7b855
@test typeof(student_year) <: Number

# ╔═╡ 2479fb6a-2e51-4e6f-8155-d47e9e3ae27d
@test typeof(student_year) <: Integer

# ╔═╡ 6d8e0985-3613-49bc-8050-91e2ee73de52
md"We often see the `==` operator when testing for equality of two numbers.  Here the `==` operator is testing for the equivalence of two types.
Alternatively, we can use the `<:` operator to ask if one type is a 'sub-type' of another 'abstract type'."

# ╔═╡ 7f031249-4413-4cd6-8ba5-b3369a1216a0
md"Now, we're ready to test that the values of the variables are consistent with what we expect."

# ╔═╡ 8ff88c72-7982-4025-9057-384a48186dc9
@test !ismissing(student_name) && length(student_name) > 3

# ╔═╡ c408ddb6-af5d-4b6a-ae81-9d338fd15f06
@test !ismissing(student_name) && occursin(' ',student_name)

# ╔═╡ 939a5680-c387-454d-a7f0-16e0d51254d8
md"In the above cell, we test that the `occursin` function retures true.  There should be a slider widget.  Try dragging it left or right.  This can be helpful when figuring out *why* a complicated test fails."

# ╔═╡ d732ab18-65b3-417a-8f8d-5f4502baa777
@test !ismissing(student_dept) && length(student_dept) >= 3

# ╔═╡ bc491991-1c10-43bb-8afe-3ed08c91d3d7
@test !ismissing(student_year) && 1 <= student_year <= 10

# ╔═╡ a328c8c6-6268-4b1f-ae1c-84ee295285c8
@test student_name == instructor

# ╔═╡ 852d7d9d-fc12-42f2-863f-af0b647979df
md"**TODO:** There are (at least) two bugs in the above tests.  Find them and fix them."

# ╔═╡ 6584e095-465d-4b3d-a383-9ac569e307cc
md"In the future, we'll separate out the tests for each exercise into a separate file (located in the `test` directory of the lab repository), so that they don't distract from the code we're writing."

# ╔═╡ 68341392-aa5a-41b7-a76a-76db1c6357f5
md"""Once you're happy with the results, it's time to save your notebook.  Look for the button in the upper right or `Ctrl`+'S'.  If it's not there, then your notebook has been autosaved.

Then read the [instructions for submitting assignments](https://psuastro528.github.io/tips/submitting) for how to commit your changes to your local git repository and push them to github.

When you're done you can look into the second notebook (ex2) in this repository to begin the next exercise."""

# ╔═╡ 51a1bf7c-6a1d-4d31-a735-de3414a97ab6
md"""## Helper code"""

# ╔═╡ 4d860e75-0631-4ef6-853c-818bafe080f4
md"""Students can ignore the code below was quietly used above."""

# ╔═╡ b9f89dc4-a8e8-4afe-9da6-9880cfcf2a37
md"""
The above cell both loads the code from packages and brings the functions (and macros) that it "exports" "into scope", meaing that we can "just call them" without having to include the package name as a prefix.
"""

# ╔═╡ 45efd82d-63e1-4016-838e-d0fbb21db99f
md"The PlutoUI package provides several cool features, but for now we just use it to add the Table of Contents in the upper right corent."

# ╔═╡ 60348e2d-9943-486c-9061-3c5c3bf88883
md"The PlutoTest package provides an alternative `@test` macro for use in Pluto notebooks."

# ╔═╡ ceff74da-4cd2-4091-af1a-8ecde8e13700
md"The PlutoTeachingTools package has several small functions that we use for things like tips."

# ╔═╡ 98ecc733-15f3-4263-821a-483dadfc8ab8
Foldable("Pro Tip!",tip(md"Normally, we'd use the `Test` module for the `@test` macro.  Julia has a large set of modules and packages, that range from very basic functionality to complex science codes.  The quality also varries widely.  Several modules (like `Test`) are included in Julia standard library, so they're already installed for us.

However, inside Pluto, it can be helpful to instead import `PlutoTest`, since it displays the results particularly nicely.  (It's an external package and it's still experimental, so if things break in the future, then we can revert to just using `Test`."))

# ╔═╡ 73a3ce7e-9967-4d11-8ab7-0e0409190975
TableOfContents()

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
PlutoTeachingTools = "661c6b06-c737-4d37-b85c-46df65de6f69"
PlutoTest = "cb4044da-4d16-4ffa-a6a3-8cad7f73ebdc"
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"

[compat]
PlutoTeachingTools = "~0.1.0"
PlutoTest = "~0.1.0"
PlutoUI = "~0.7.9"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

[[Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[HypertextLiteral]]
git-tree-sha1 = "1e3ccdc7a6f7b577623028e0095479f4727d8ec1"
uuid = "ac1192a8-f4b3-4bfe-ba22-af5b92cd3ab2"
version = "0.8.0"

[[InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[JSON]]
deps = ["Dates", "Mmap", "Parsers", "Unicode"]
git-tree-sha1 = "81690084b6198a2e1da36fcfda16eeca9f9f24e4"
uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
version = "0.21.1"

[[LaTeXStrings]]
git-tree-sha1 = "c7f1c695e06c01b95a67f0cd1d34994f3e7db104"
uuid = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
version = "1.2.1"

[[Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"

[[Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"

[[Mmap]]
uuid = "a63ad114-7e13-5084-954f-fe012c677804"

[[Parsers]]
deps = ["Dates"]
git-tree-sha1 = "94bf17e83a0e4b20c8d77f6af8ffe8cc3b386c0a"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "1.1.1"

[[PlutoTeachingTools]]
deps = ["LaTeXStrings", "Markdown", "PlutoUI", "Random"]
git-tree-sha1 = "265980831960aabe7e1f5ae47c898a8459588ee7"
uuid = "661c6b06-c737-4d37-b85c-46df65de6f69"
version = "0.1.3"

[[PlutoTest]]
deps = ["HypertextLiteral", "InteractiveUtils", "Markdown", "Test"]
git-tree-sha1 = "3479836b31a31c29a7bac1f09d95f9c843ce1ade"
uuid = "cb4044da-4d16-4ffa-a6a3-8cad7f73ebdc"
version = "0.1.0"

[[PlutoUI]]
deps = ["Base64", "Dates", "InteractiveUtils", "JSON", "Logging", "Markdown", "Random", "Reexport", "Suppressor"]
git-tree-sha1 = "44e225d5837e2a2345e69a1d1e01ac2443ff9fcb"
uuid = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
version = "0.7.9"

[[Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[Random]]
deps = ["Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[Reexport]]
git-tree-sha1 = "5f6c21241f0f655da3952fd60aa18477cf96c220"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "1.1.0"

[[Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[Suppressor]]
git-tree-sha1 = "a819d77f31f83e5792a76081eee1ea6342ab8787"
uuid = "fd094767-a336-5f1f-9728-57cf17d0bbfb"
version = "0.2.0"

[[Test]]
deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[[Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"
"""

# ╔═╡ Cell order:
# ╟─5b826e17-ad51-43d0-9756-89a9fd0ea834
# ╟─acb10cb1-d03d-47c5-90dc-7216e5110ed2
# ╟─c15ce452-494f-407c-9c08-02d547929e6a
# ╠═25d5497c-3a67-4495-9f37-0bc8f2a64e68
# ╟─774b2e4a-6cd7-4c9a-920d-b4f38354d42d
# ╠═5917df75-0dab-4778-917e-91bbfddbe03a
# ╟─0d15ba33-f2c4-46d6-9b08-9181d40596cd
# ╠═dc0bf251-b3a6-42f7-84f6-796ac3761b2d
# ╠═dbfceab2-994d-45a1-97ff-d978827bd129
# ╟─e2bb48e5-a2b2-42d7-a8e5-5c7532d62e93
# ╠═cda49166-b9fb-4e68-8b3e-b97845c7bf3b
# ╠═0c5d5b14-6c9a-4ac9-8126-0d6e1edba593
# ╟─a3450384-ccfc-4bda-af52-24ed8f429cf8
# ╠═33350a21-7971-406e-b316-6539eca72bf4
# ╟─f15643a4-39d5-4492-98d3-9db4f4cb0d3b
# ╟─1bc2edaf-a0a7-4c4e-96c2-25872dd30e2e
# ╟─c05d62b5-6063-41ca-9823-2499b2d4fbc4
# ╟─8b666639-e3d0-4535-803d-40ddc4360a9c
# ╟─ec4ec6e7-9548-4bd6-970a-48fefa9861d7
# ╟─ac9e6787-5c42-44d6-907e-e7033e7ab4bc
# ╟─087b40c4-7d59-4f9a-b584-65c2aa9f9578
# ╟─eefa92a1-1d4c-4536-abd3-1c1af98e146f
# ╠═78f81544-fead-4649-b306-48fb1463d397
# ╠═3d224633-0eb8-407c-b44c-0d9ec927dfc6
# ╠═2fe3f4e8-5288-4f61-9346-7a455c2abc85
# ╟─b49789b2-49c0-49e3-b81f-db04e60f9d23
# ╠═34f767c9-f95f-4968-9423-ad87bb456775
# ╠═d6417bdc-9e36-4d2b-8b3e-f5ddafd82437
# ╠═89ac58ea-b782-4dcc-802f-d6d063e7b855
# ╠═2479fb6a-2e51-4e6f-8155-d47e9e3ae27d
# ╟─6d8e0985-3613-49bc-8050-91e2ee73de52
# ╟─7f031249-4413-4cd6-8ba5-b3369a1216a0
# ╠═8ff88c72-7982-4025-9057-384a48186dc9
# ╠═c408ddb6-af5d-4b6a-ae81-9d338fd15f06
# ╟─939a5680-c387-454d-a7f0-16e0d51254d8
# ╠═d732ab18-65b3-417a-8f8d-5f4502baa777
# ╠═bc491991-1c10-43bb-8afe-3ed08c91d3d7
# ╠═a328c8c6-6268-4b1f-ae1c-84ee295285c8
# ╟─852d7d9d-fc12-42f2-863f-af0b647979df
# ╟─6584e095-465d-4b3d-a383-9ac569e307cc
# ╟─68341392-aa5a-41b7-a76a-76db1c6357f5
# ╟─51a1bf7c-6a1d-4d31-a735-de3414a97ab6
# ╟─4d860e75-0631-4ef6-853c-818bafe080f4
# ╠═ab976d70-f59d-4dcd-ac89-771606a74108
# ╟─b9f89dc4-a8e8-4afe-9da6-9880cfcf2a37
# ╠═45efd82d-63e1-4016-838e-d0fbb21db99f
# ╟─60348e2d-9943-486c-9061-3c5c3bf88883
# ╟─ceff74da-4cd2-4091-af1a-8ecde8e13700
# ╟─98ecc733-15f3-4263-821a-483dadfc8ab8
# ╠═73a3ce7e-9967-4d11-8ab7-0e0409190975
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
# Communicator

export UniqueID, Communicator, rank

import CUDAdrv: device

const NCCL_UNIQUE_ID_BYTES = 128
const ncclUniqueId_t = NTuple{NCCL_UNIQUE_ID_BYTES, Cchar}

struct UniqueID
    internal::ncclUniqueId_t

    function UniqueID()
        buf = zeros(Cchar, NCCL_UNIQUE_ID_BYTES)
        ncclGetUniqueId(buf)
        new(Tuple(buf))
    end
end

Base.convert(::Type{ncclUniqueId_t}, id::UniqueID) = id.internal

const ncclComm_t = Ptr{Cvoid}

struct Communicator
    handle::ncclComm_t
end

# creates a new communicator (multi thread/process version)
function Communicator(nranks, comm_id, rank)
    handle_ref = Ref{ncclComm_t}(C_NULL)
    ncclCommInitRank(handle_ref, nranks, comm_id.internal, rank)
    c = Communicator(handle_ref[])
    finalizer(c) do x
        ncclCommDestroy(x.handle)
    end
    return c
end

# creates a clique of communicators (single process version)
function Communicator(devices::Union{CUDAdrv.DeviceSet,AbstractVector})
    ndev    = length(devices)
    comms   = Vector{ncclComm_t}(undef, ndev)
    devlist = Cint.([i-1 for i in 1:ndev])
    ncclCommInitAll(comms, ndev, devlist)
    cs = Communicator.(comms)
    finalizer(cs) do xs
        for x in xs
            ncclCommDestroy(x.handle)
        end
    end
    return cs
end

function device(comm::Communicator)
    dev_ref = Ref{Cint}(C_NULL)
    ncclCommCuDevice(comm.handle, dev_ref)
    return dev_ref[]
end

function Base.size(comm::Communicator)
    size_ref = Ref{Cint}(C_NULL)
    ncclCommCount(comm.handle, size_ref)
    return size_ref[]
end

function rank(comm::Communicator)
    rank_ref = Ref{Cint}(C_NULL)
    ncclCommUserRank(comm.handle, rank_ref)
    return rank_ref[]
end

ext: jl


text:
module BrickMesh

export brickmesh, centroidtocode, connectmesh, partition, connectmeshfull

using MPI

"""
    linearpartition(n, p, np)

Partition the range `1:n` into `np` pieces and return the `p`th piece as a
range.

This will provide an equal partition when `n` is divisible by `np` and
otherwise the ranges will have lengths of either `floor(Int, n/np)` or
`ceil(Int, n/np)`.
"""
linearpartition(n, p, np) =
    range(div((p - 1) * n, np) + 1, stop = div(p * n, np))

"""
    hilbertcode(Y::AbstractArray{T}; bits=8sizeof(T)) where T

Given an array of axes coordinates `Y` stored as `bits`-bit integers
the function returns the Hilbert integer `H`.

The encoding of the Hilbert integer is best described by example.
If 5-bits are used from each of 3 coordinates then the function performs

     X[2]|                       H[0] = A B C D E
         | /X[1]       ------->  H[1] = F G H I J
    axes |/                      H[2] = K L M N O
         0------ X[0]                   high low

where the 15-bit Hilbert integer = `A B C D E F G H I J K L M N O` is stored
in `H`

This function is based on public domain code from John Skilling which can be
found in [Skilling2004](@cite).
"""
function hilbertcode(Y::AbstractArray{T}; bits = 8 * sizeof(T)) where {T}
    # Below is Skilling's AxestoTranspose
    X = deepcopy(Y)
    n = length(X)
    M = one(T) << (bits - 1)

    Q = M
    for j in 1:(bits - 1)
        P = Q - one(T)
        for i in 1:n
            if X[i] & Q != zero(T)
                X[1] ⊻= P
            else
                t = (X[1] ⊻ X[i]) & P
                X[1] ⊻= t
                X[i] ⊻= t
            end
        end
        Q >>>= one(T)
    end

    for i in 2:n
        X[i] ⊻= X[i - 1]
    end

    t = zero(T)
    Q = M
    for j in 1:(bits - 1)
        if X[n] & Q != zero(T)
            t ⊻= Q - one(T)
        end
        Q >>>= one(T)
    end

    for i in 1:n
        X[i] ⊻= t
    end

    # Below we transpose X and store it in H, i.e.:
    #
    #   X[0] = A D G J M               H[0] = A B C D E
    #   X[1] = B E H K N   <------->   H[1] = F G H I J
    #   X[2] = C F I L O               H[2] = K L M N O
    #
    # The 15-bit Hilbert integer is then = A B C D E F G H I J K L M N O
    H = zero(X)
    for i in 0:(n - 1), j in 0:(bits - 1)
        k = i * bits + j
        bit = (X[n - mod(k, n)] >>> div(k, n)) & one(T)
        H[n - i] |= (bit << j)
    end

    return H
end

"""
    centroidtocode(comm::MPI.Comm, elemtocorner; coortocode, CT)

Returns a code for each element based on its centroid.

These element codes can be used to determine a linear ordering for the
partition function.

The communicator `comm` is used to calculate the bounding box for representing
the centroids in coordinates of type `CT`, defaulting to `CT=UInt64`.  These
integer coordinates are converted to a code using the function `coortocode`,
which defaults to `hilbertcode`.

The array containing the element corner coordinates, `elemtocorner`, is used
to compute the centroids.  `elemtocorner` is a dimension by number of corners
by number of elements array.
"""
function centroidtocode(
    comm::MPI.Comm,
    elemtocorner;
    coortocode = hilbertcode,
    CT = UInt64,
)
    (d, nvert, nelem) = size(elemtocorner)

    centroids = sum(elemtocorner, dims = 2) ./ nvert
    T = eltype(centroids)

    centroidmin =
        (nelem > 0) ? minimum(centroids, dims = 3) : fill(typemax(T), d)
    centroidmax =
        (nelem > 0) ? maximum(centroids, dims = 3) : fill(typemin(T), d)

    centroidmin = MPI.Allreduce(centroidmin, min, comm)
    centroidmax = MPI.Allreduce(centroidmax, max, comm)
    centroidsize = centroidmax - centroidmin

    # Fix centroidsize to be nonzero.  It can be zero for a couple of reasons.
    # For example, it will be zero if we have just one element.
    if iszero(centroidsize)
        centroidsize = ones(T, d)
    else
        for i in 1:d
            if iszero(centroidsize[i])
                centroidsize[i] = maximum(centroidsize)
            end
        end
    end

    code = Array{CT}(undef, d, nelem)
    for e in 1:nelem
        c = (centroids[:, 1, e] .- centroidmin) ./ centroidsize
        X =
            CT.(floor.(
                typemax(CT) .* BigFloat.(c; precision = 16 * sizeof(CT)),
            ))
        code[:, e] = coortocode(X)
    end

    code
end

"""
    brickmesh(x, periodic; part=1, numparts=1; boundary)

Generate a brick mesh with coordinates given by the tuple `x` and the
periodic dimensions given by the `periodic` tuple.

The brick can optionally be partitioned into `numparts` and this returns
partition `part`.  This is a simple Cartesian partition and further
partitioning (e.g, based on a space-filling curve) should be done before the
mesh is used for computation.

By default boundary faces will be marked with a one and other faces with a
zero.  Specific boundary numbers can also be passed for each face of the brick
in `boundary`.  This will mark the nonperiodic brick faces with the given
boundary number.

# Examples

We can build a 3 by 2 element two-dimensional mesh that is periodic in the
\$x_2\$-direction with
```jldoctest brickmesh
julia> (elemtovert, elemtocoord, elemtobndy, faceconnections) =
        brickmesh((2:5,4:6), (false,true); boundary=((1,2), (3,4)));
```
This returns the mesh structure for

             x_2

              ^
              |
             6-  9----10----11----12
              |  |     |     |     |
              |  |  4  |  5  |  6  |
              |  |     |     |     |
             5-  5-----6-----7-----8
              |  |     |     |     |
              |  |  1  |  2  |  3  |
              |  |     |     |     |
             4-  1-----2-----3-----4
              |
              +--|-----|-----|-----|--> x_1
                 2     3     4     5

The (number of corners by number of elements) array `elemtovert` gives the
global vertex number for the corners of each element.
```jldoctest brickmesh
julia> elemtovert
4×6 Array{Int64,2}:
 1  2  3   5   6   7
 2  3  4   6   7   8
 5  6  7   9  10  11
 6  7  8  10  11  12
```
Note that the vertices are listed in Cartesian order.

The (dimension by number of corners by number of elements) array `elemtocoord`
gives the coordinates of the corners of each element.
```jldoctes brickmesh
julia> elemtocoord
2×4×6 Array{Int64,3}:
[:, :, 1] =
 2  3  2  3
 4  4  5  5

[:, :, 2] =
 3  4  3  4
 4  4  5  5

[:, :, 3] =
 4  5  4  5
 4  4  5  5

[:, :, 4] =
 2  3  2  3
 5  5  6  6

[:, :, 5] =
 3  4  3  4
 5  5  6  6

[:, :, 6] =
 4  5  4  5
 5  5  6  6
```

The (number of faces by number of elements) array `elemtobndy` gives the
boundary number for each face of each element.  A zero will be given for
connected faces.
```jldoctest brickmesh
julia> elemtobndy
4×6 Array{Int64,2}:
 1  0  0  1  0  0
 0  0  2  0  0  2
 0  0  0  0  0  0
 0  0  0  0  0  0
```
Note that the faces are listed in Cartesian order.

Finally, the periodic face connections are given in `faceconnections` which is a
list of arrays, one for each connection.
Each array in the list is given in the format `[e, f, vs...]` where
 - `e`  is the element number;
 - `f`  is the face number; and
 - `vs` is the global vertices that face associated with.
I the example
```jldoctest brickmesh
julia> faceconnections
3-element Array{Array{Int64,1},1}:
 [4, 4, 1, 2]
 [5, 4, 2, 3]
 [6, 4, 3, 4]
```
we see that face `4` of element `5` is associated with vertices `[2 3]` (the
vertices for face `1` of element `2`).
"""
function brickmesh(
    x,
    periodic;
    part = 1,
    numparts = 1,
    boundary = ntuple(j -> (1, 1), length(x)),
)
    if boundary isa Matrix
        boundary = tuple(mapslices(x -> tuple(x...), boundary, dims = 1)...)
    end

    @assert length(x) == length(periodic)
    @assert length(x) >= 1
    @assert 1 <= part <= numparts

    T = promote_type(eltype.(x)...)
    d = length(x)
    nvert = 2^d
    nface = 2d

    nelemdim = length.(x) .- 1
    elemlocal = linearpartition(prod(nelemdim), part, numparts)

    elemtovert = Array{Int}(undef, nvert, length(elemlocal))
    elemtocoord = Array{T}(undef, d, nvert, length(elemlocal))
    elemtobndy = zeros(Int, nface, length(elemlocal))
    faceconnections = Array{Array{Int, 1}}(undef, 0)

    verts = LinearIndices(ntuple(j -> 1:length(x[j]), d))
    elems = CartesianIndices(ntuple(j -> 1:(length(x[j]) - 1), d))

    p = reshape(1:nvert, ntuple(j -> 2, d))
    fmask = hcat((
        p[ntuple(
            j ->
                (j == div(f - 1, 2) + 1) ?
                ((mod(f - 1, 2) + 1):(mod(f - 1, 2) + 1)) : (:),
            d,
        )...,][:] for f in 1:nface
    )...)

    for (e, ec) in enumerate(elems[elemlocal])
        corners = CartesianIndices(ntuple(j -> ec[j]:(ec[j] + 1), d))
        for (v, vc) in enumerate(corners)
            elemtovert[v, e] = verts[vc]

            for j in 1:d
                elemtocoord[j, v, e] = x[j][vc[j]]
            end
        end

        for i in 1:d
            if !periodic[i] && ec[i] == 1
                elemtobndy[2 * (i - 1) + 1, e] = boundary[i][1]
            end
            if !periodic[i] && ec[i] == nelemdim[i]
                elemtobndy[2 * (i - 1) + 2, e] = boundary[i][2]
            end
        end

        for i in 1:d
            if periodic[i] && ec[i] == nelemdim[i]
                neighcorners = CartesianIndices(ntuple(
                    j -> (i == j) ? (1:2) : (ec[j]:(ec[j] + 1)),
                    d,
                ))
                push!(
                    faceconnections,
                    vcat(e, 2i, verts[neighcorners[fmask[:, 2i - 1]]]),
                )
            end
        end
    end

    (elemtovert, elemtocoord, elemtobndy, faceconnections)
end

"""
    parallelsortcolumns(comm::MPI.Comm, A;
                        alg::Base.Sort.Algorithm=Base.Sort.DEFAULT_UNSTABLE,
                        lt=isless,
                        by=identity,
                        rev::Union{Bool,Nothing}=nothing)

Sorts the columns of the distributed matrix `A`.

See the documentation of `sort!` for a description of the keyword arguments.

This function assumes `A` has the same number of rows on each MPI rank but can
have a different number of columns.
"""
function parallelsortcolumns(
    comm::MPI.Comm,
    A;
    alg::Base.Sort.Algorithm = Base.Sort.DEFAULT_UNSTABLE,
    lt = isless,
    by = identity,
    rev::Union{Bool, Nothing} = nothing,
)

    m, n = size(A)
    T = eltype(A)

    csize = MPI.Comm_size(comm)
    crank = MPI.Comm_rank(comm)
    croot = 0

    A = sortslices(A, dims = 2, alg = alg, lt = lt, by = by, rev = rev)

    npivots = clamp(n, 0, csize)
    pivots = T[A[i, div(n * p, npivots) + 1] for i in 1:m, p in 0:(npivots - 1)]
    pivotcounts = MPI.Allgather(Cint(length(pivots)), comm)
    pivots = MPI.Allgatherv!(
        pivots,
        VBuffer(similar(pivots, sum(pivotcounts)), pivotcounts),
        comm,
    )
    pivots = reshape(pivots, m, div(length(pivots), m))
    pivots =
        sortslices(pivots, dims = 2, alg = alg, lt = lt, by = by, rev = rev)

    # if we don't have any pivots then we must have zero columns
    if size(pivots) == (m, 0)
        return A
    end

    pivots = [
        pivots[i, div(div(length(pivots), m) * r, csize) + 1]
        for i in 1:m, r in 0:(csize - 1)
    ]

    cols = map(i -> view(A, :, i), 1:n)
    senddispls = Cint[
        (
            searchsortedfirst(cols, pivots[:, i], lt = lt, by = by, rev = rev) - 1
        ) * m for i in 1:csize
    ]
    sendcounts = Cint[
        (i == csize ? n * m : senddispls[i + 1]) - senddispls[i]
        for i in 1:csize
    ]

    recvcounts = similar(sendcounts)
    MPI.Alltoall!(UBuffer(sendcounts, 1), MPI.UBuffer(recvcounts, 1), comm)
    B = similar(A, sum(recvcounts))
    MPI.Alltoallv!(
        VBuffer(A, sendcounts, senddispls),
        VBuffer(B, recvcounts),
        comm,
    )
    B = reshape(B, m, div(length(B), m))
    sortslices(B, dims = 2, alg = alg, lt = lt, by = by, rev = rev)
end

"""
    getpartition(comm::MPI.Comm, elemtocode)

Returns an equally weighted partition of a distributed set of elements by
sorting their codes given in `elemtocode`.

The codes for each element, `elemtocode`, are given as an array with a single
entry per local element or as a matrix with a column for each local element.

The partition is returned as a tuple three parts:

 - `partsendorder`: permutation of elements into sending order
 - `partsendstarts`: start entries in the send array for each rank
 - `partrecvstarts`: start entries in the receive array for each rank

Note that both `partsendstarts` and `partrecvstarts` are of length
`MPI.Comm_size(comm)+1` where the last entry has the total number of elements
to send or receive, respectively.
"""
getpartition(comm::MPI.Comm, elemtocode::AbstractVector) =
    getpartition(comm, reshape(elemtocode, 1, length(elemtocode)))

function getpartition(comm::MPI.Comm, elemtocode::AbstractMatrix)
    (ncode, nelem) = size(elemtocode)

    csize = MPI.Comm_size(comm)
    crank = MPI.Comm_rank(comm)

    CT = eltype(elemtocode)

    A = CT[
        elemtocode                                 # code
        collect(CT, 1:nelem)'                      # original element number
        fill(CT(MPI.Comm_rank(comm)), (1, nelem))  # original rank
        fill(typemax(CT), (1, nelem))
    ]              # new rank
    m, n = size(A)

    # sort by just code
    A = parallelsortcolumns(comm, A)

    # count the distribution of A
    counts = MPI.Allgather(last(size(A)), comm)
    starts = ones(Int, csize + 1)
    for i in 1:csize
        starts[i + 1] = counts[i] + starts[i]
    end

    # loop to determine new rank
    j = range(starts[crank + 1], stop = starts[crank + 2] - 1)
    for r in 0:(csize - 1)
        k = linearpartition(starts[end] - 1, r + 1, csize)
        o = intersect(k, j) .- (starts[crank + 1] - 1)
        A[ncode + 3, o] .= r
    end

    # sort by original rank and code
    A = sortslices(A, dims = 2, by = x -> x[[ncode + 2, (1:ncode)...]])

    # count number of elements that are going to be sent
    sendcounts = zeros(Cint, csize)
    for i in 1:last(size(A))
        sendcounts[A[ncode + 2, i] + 1] += m
    end

    # communicate columns of A to original rank
    recvcounts = similar(sendcounts)
    MPI.Alltoall!(UBuffer(sendcounts, 1), UBuffer(recvcounts, 1), comm)
    B = similar(A, sum(recvcounts))
    MPI.Alltoallv!(VBuffer(A, sendcounts), VBuffer(B, recvcounts), comm)
    B = reshape(B, m, div(length(B), m))

    # check to make sure we didn't drop any elements
    @assert nelem == n == size(B)[2]

    partsendcounts = zeros(Cint, csize)
    for i in 1:last(size(B))
        partsendcounts[B[ncode + 3, i] + 1] += 1
    end
    partsendstarts = ones(Int, csize + 1)
    for i in 1:csize
        partsendstarts[i + 1] = partsendcounts[i] + partsendstarts[i]
    end

    partsendorder = Int.(B[ncode + 1, :])

    partrecvcounts = similar(partsendcounts)
    MPI.Alltoall!(UBuffer(partsendcounts, 1), UBuffer(partrecvcounts, 1), comm)

    partrecvstarts = ones(Int, csize + 1)
    for i in 1:csize
        partrecvstarts[i + 1] = partrecvcounts[i] + partrecvstarts[i]
    end

    partsendorder, partsendstarts, partrecvstarts
end

"""
    partition(comm::MPI.Comm, elemtovert, elemtocoord, elemtobndy,
              faceconnections)

This function takes in a mesh (as returned for example by `brickmesh`) and
returns a Hilbert curve based partitioned mesh.
"""
function partition(
    comm::MPI.Comm,
    elemtovert,
    elemtocoord,
    elemtobndy,
    faceconnections,
    globord = [],
)
    (d, nvert, nelem) = size(elemtocoord)

    csize = MPI.Comm_size(comm)
    crank = MPI.Comm_rank(comm)

    nface = 2d
    nfacevert = 2^(d - 1)

    # Here we expand the list of face connections into a structure that is easy
    # to partition.  The cost is extra memory transfer.  If this becomes a
    # bottleneck something more efficient may be implemented.
    #
    elemtofaceconnect =
        zeros(eltype(eltype(faceconnections)), nfacevert, nface, nelem)
    for fc in faceconnections
        elemtofaceconnect[:, fc[2], fc[1]] = fc[3:end]
    end

    elemtocode = centroidtocode(comm, elemtocoord; CT = UInt64)
    sendorder, sendstarts, recvstarts = getpartition(comm, elemtocode)

    elemtovert = elemtovert[:, sendorder]
    elemtocoord = elemtocoord[:, :, sendorder]
    elemtobndy = elemtobndy[:, sendorder]
    elemtofaceconnect = elemtofaceconnect[:, :, sendorder]

    if !isempty(globord)
        globord = globord[sendorder]
    end

    sendcounts = diff(sendstarts)
    recvcounts = similar(sendcounts)
    MPI.Alltoall!(UBuffer(sendcounts, 1), UBuffer(recvcounts, 1), comm)

    newelemtovert = similar(elemtovert, nvert * sum(recvcounts))
    MPI.Alltoallv!(
        VBuffer(elemtovert, Cint(nvert) .* sendcounts),
        VBuffer(newelemtovert, Cint(nvert) .* recvcounts),
        comm,
    )

    newelemtocoord = similar(elemtocoord, d * nvert * sum(recvcounts))
    MPI.Alltoallv!(
        VBuffer(elemtocoord, Cint(d * nvert) .* sendcounts),
        VBuffer(newelemtocoord, Cint(d * nvert) .* recvcounts),
        comm,
    )

    newelemtobndy = similar(elemtobndy, nface * sum(recvcounts))
    MPI.Alltoallv!(
        VBuffer(elemtobndy, Cint(nface) .* sendcounts),
        VBuffer(newelemtobndy, Cint(nface) .* recvcounts),
        comm,
    )

    newelemtofaceconnect =
        similar(elemtofaceconnect, nfacevert * nface * sum(recvcounts))
    MPI.Alltoallv!(
        VBuffer(elemtofaceconnect, Cint(nfacevert * nface) .* sendcounts),
        VBuffer(newelemtofaceconnect, Cint(nfacevert * nface) .* recvcounts),
        comm,
    )

    if !isempty(globord)
        newglobord = similar(globord, sum(recvcounts))
        MPI.Alltoallv!(
            VBuffer(globord, sendcounts),
            VBuffer(newglobord, recvcounts),
            comm,
        )
    else
        newglobord = similar(globord)
    end

    newnelem = recvstarts[end] - 1
    newelemtovert = reshape(newelemtovert, nvert, newnelem)
    newelemtocoord = reshape(newelemtocoord, d, nvert, newnelem)
    newelemtobndy = reshape(newelemtobndy, nface, newnelem)
    newelemtofaceconnect =
        reshape(newelemtofaceconnect, nfacevert, nface, newnelem)

    # reorder local elements based on code of new elements
    A = UInt64[
        centroidtocode(comm, newelemtocoord; CT = UInt64)
        collect(1:newnelem)'
    ]
    A = sortslices(A, dims = 2)
    newsortorder = view(A, d + 1, :)

    newelemtovert = newelemtovert[:, newsortorder]
    newelemtocoord = newelemtocoord[:, :, newsortorder]
    newelemtobndy = newelemtobndy[:, newsortorder]
    newelemtofaceconnect = newelemtofaceconnect[:, :, newsortorder]

    newfaceconnections = similar(faceconnections, 0)
    for e in 1:newnelem, f in 1:nface
        if newelemtofaceconnect[1, f, e] > 0
            push!(newfaceconnections, vcat(e, f, newelemtofaceconnect[:, f, e]))
        end
    end

    if !isempty(globord)
        newglobord = newglobord[newsortorder]
    end

    (
        newelemtovert,
        newelemtocoord,
        newelemtobndy,
        newfaceconnections,
        newglobord,
    )#sendorder)
end

"""
    minmaxflip(x, y)

Returns `x, y` sorted lowest to highest and a bool that indicates if a swap
was needed.
"""
minmaxflip(x, y) = y < x ? (y, x, true) : (x, y, false)

"""
    vertsortandorder(a)

Returns `(a)` and an ordering `o==0`.
"""
vertsortandorder(a) = ((a,), 1)

"""
    vertsortandorder(a, b)

Returns sorted vertex numbers `(a,b)` and an ordering `o` depending on the
order needed to sort the elements.  This ordering is given below including the
vetex ordering for faces.

    o=    0      1

        (a,b)  (b,a)

          a      b
          |      |
          |      |
          b      a
"""
function vertsortandorder(a, b)
    a, b, s1 = minmaxflip(a, b)
    o = s1 ? 2 : 1
    ((a, b), o)
end

"""
    vertsortandorder(a, b, c)

Returns sorted vertex numbers `(a,b,c)` and an ordering `o` depending on the
order needed to sort the elements.  This ordering is given below including the
vetex ordering for faces.

    o=     1         2         3         4         5         6

        (a,b,c)   (c,a,b)   (b,c,a)   (b,a,c)   (c,b,a)   (a,c,b)

          /c\\      /b\\      /a\\      /c\\      /a\\      /b\\
         /   \\    /   \\    /   \\    /   \\    /   \\    /   \\
        /a___b\\  /c___a\\  /b___c\\  /b___a\\  /c___b\\  /a___c\\
"""
function vertsortandorder(a, b, c)
    # Use a (Bose-Nelson Algorithm based) sorting network from
    # <http://pages.ripco.net/~jgamble/nw.html> to sort the vertices.
    b, c, s1 = minmaxflip(b, c)
    a, c, s2 = minmaxflip(a, c)
    a, b, s3 = minmaxflip(a, b)

    if !s1 && !s2 && !s3
        o = 1
    elseif !s1 && s2 && s3
        o = 2
    elseif s1 && !s2 && s3
        o = 3
    elseif !s1 && !s2 && s3
        o = 4
    elseif s1 && s2 && s3
        o = 5
    elseif s1 && !s2 && !s3
        o = 6
    else
        error("Problem finding vertex ordering $((a,b,c)) with flips
              $((s1,s2,s3))")
    end

    ((a, b, c), o)
end

"""
    vertsortandorder(a, b, c, d)

Returns sorted vertex numbers `(a,b,c,d)` and an ordering `o` depending on the
order needed to sort the elements.  This ordering is given below including the
vetex ordering for faces.

    o=   1      2      3      4      5      6      7      8

       (a,b,  (a,c,  (b,a,  (b,d,  (c,a,  (c,d,  (d,b,  (d,c,
        c,d)   b,d)   c,d)   a,c)   d,b)   a,b)   c,a)   b,a)

       c---d  b---d  c---d  a---c  d---b  a---b  c---a  b---a
       |   |  |   |  |   |  |   |  |   |  |   |  |   |  |   |
       a---b  a---c  b---a  b---d  c---a  c---d  d---b  d---c
"""
function vertsortandorder(a, b, c, d)
    # Use a (Bose-Nelson Algorithm based) sorting network from
    # <http://pages.ripco.net/~jgamble/nw.html> to sort the vertices.
    a, b, s1 = minmaxflip(a, b)
    c, d, s2 = minmaxflip(c, d)
    a, c, s3 = minmaxflip(a, c)
    b, d, s4 = minmaxflip(b, d)
    b, c, s5 = minmaxflip(b, c)

    if !s1 && !s2 && !s3 && !s4 && !s5
        o = 1
    elseif !s1 && !s2 && !s3 && !s4 && s5
        o = 2
    elseif s1 && !s2 && !s3 && !s4 && !s5
        o = 3
    elseif !s1 && !s2 && s3 && s4 && s5
        o = 4
    elseif s1 && s2 && !s3 && !s4 && s5
        o = 5
    elseif !s1 && !s2 && s3 && s4 && !s5
        o = 6
    elseif s1 && s2 && s3 && s4 && s5
        o = 7
    elseif s1 && s2 && s3 && s4 && !s5
        o = 8
    else
        # FIXME: some possible orientations are missing since there are a total of
        # 24. Missing orientations:
        #=
           d---c  d---c  b---c  a---d  c---b  d---a  a---b  b---a
           |   |  |   |  |   |  |   |  |   |  |   |  |   |  |   |
           a---b  b---a  a---d  b---c  d---a  c---b  d---c  c---d

           c---b  d---b  b---d  b---c  c---a  d---a  a---d  a---c
           |   |  |   |  |   |  |   |  |   |  |   |  |   |  |   |
           a---d  a---c  c---a  d---a  b---d  b---c  c---b  d---b
        =#
        error("Problem finding vertex ordering $((a,b,c,d))
                with flips $((s1,s2,s3,s4,s5))")
    end

    ((a, b, c, d), o)
end

"""
    connectmesh(comm::MPI.Comm, elemtovert, elemtocoord, elemtobndy,
                faceconnections)

This function takes in a mesh (as returned for example by `brickmesh`) and
returns a connected mesh.  This returns a `NamedTuple` of:

 - `elems` the range of element indices
 - `realelems` the range of real (aka nonghost) element indices
 - `ghostelems` the range of ghost element indices
 - `ghostfaces` ghost element to face is received;
   `ghostfaces[f,ge] == true` if face `f` of ghost element `ge` is received.
 - `sendelems` an array of send element indices
 - `sendfaces` send element to face is sent;
   `sendfaces[f,se] == true` if face `f` of send element `se` is sent.
 - `elemtocoord` element to vertex coordinates; `elemtocoord[d,i,e]` is the
    `d`th coordinate of corner `i` of element `e`
 - `elemtoelem` element to neighboring element; `elemtoelem[f,e]` is the
   number of the element neighboring element `e` across face `f`.  If there is
   no neighboring element then `elemtoelem[f,e] == e`.
 - `elemtoface` element to neighboring element face; `elemtoface[f,e]` is the
   face number of the element neighboring element `e` across face `f`.  If
   there is no neighboring element then `elemtoface[f,e] == f`.
 - `elemtoordr` element to neighboring element order; `elemtoordr[f,e]` is the
   ordering number of the element neighboring element `e` across face `f`.  If
   there is no neighboring element then `elemtoordr[f,e] == 1`.
 - `elemtobndy` element to bounday number; `elemtobndy[f,e]` is the
   boundary number of face `f` of element `e`.  If there is a neighboring
   element then `elemtobndy[f,e] == 0`.
 - `nabrtorank` a list of the MPI ranks for the neighboring processes
 - `nabrtorecv` a range in ghost elements to receive for each neighbor
 - `nabrtosend` a range in `sendelems` to send for each neighbor

"""
function connectmesh(
    comm::MPI.Comm,
    elemtovert,
    elemtocoord,
    elemtobndy,
    faceconnections;
    dim = size(elemtocoord, 1),
)
    d = dim
    (coorddim, nvert, nelem) = size(elemtocoord)
    nface, nfacevert = 2d, 2^(d - 1)

    p = reshape(1:nvert, ntuple(j -> 2, d))
    fmask = hcat((
        p[ntuple(
            j ->
                (j == div(f - 1, 2) + 1) ?
                ((mod(f - 1, 2) + 1):(mod(f - 1, 2) + 1)) : (:),
            d,
        )...][:] for f in 1:nface
    )...)

    csize = MPI.Comm_size(comm)
    crank = MPI.Comm_rank(comm)

    VT = eltype(elemtovert)
    A = Array{VT}(undef, nfacevert + 8, nface * nelem)

    MR, ME, MF, MO, NR, NE, NF, NO = nfacevert .+ (1:8)
    for e in 1:nelem
        v = reshape(elemtovert[:, e], ntuple(j -> 2, d))
        for f in 1:nface
            j = (e - 1) * nface + f
            fv, o = vertsortandorder(v[fmask[:, f]]...)
            A[1:nfacevert, j] .= fv
            A[MR, j] = crank
            A[ME, j] = e
            A[MF, j] = f
            A[MO, j] = o
            A[NR, j] = typemax(VT)
            A[NE, j] = typemax(VT)
            A[NF, j] = typemax(VT)
            A[NO, j] = typemax(VT)
        end
    end

    # use neighboring vertices for connected faces
    for fc in faceconnections
        e = fc[1]
        f = fc[2]
        v = fc[3:end]
        j = (e - 1) * nface + f
        fv, o = vertsortandorder(v...)
        A[1:nfacevert, j] .= fv
        A[MO, j] = o
    end

    A = parallelsortcolumns(comm, A, by = x -> x[1:nfacevert])
    m, n = size(A)

    # match faces
    j = 1
    while j <= n
        if j + 1 <= n && A[1:nfacevert, j] == A[1:nfacevert, j + 1]
            # found connected face
            A[NR:NO, j] = A[MR:MO, j + 1]
            A[NR:NO, j + 1] = A[MR:MO, j]
            j += 2
        else
            # found unconnect face
            A[NR:NO, j] = A[MR:MO, j]
            j += 1
        end
    end

    A = sortslices(A, dims = 2, by = x -> (x[MR], x[NR], x[ME], x[MF]))

    # count number of elements that are going to be sent
    sendcounts = zeros(Cint, csize)
    for i in 1:last(size(A))
        sendcounts[A[MR, i] + 1] += m
    end
    sendstarts = ones(Int, csize + 1)
    for i in 1:csize
        sendstarts[i + 1] = sendcounts[i] + sendstarts[i]
    end

    # communicate columns of A to original rank
    recvcounts = similar(sendcounts)
    MPI.Alltoall!(UBuffer(sendcounts, 1), UBuffer(recvcounts, 1), comm)

    B = similar(A, sum(recvcounts))
    MPI.Alltoallv!(VBuffer(A, sendcounts), VBuffer(B, recvcounts), comm)
    B = reshape(B, m, nface * nelem)

    # get element sending information
    B = sortslices(B, dims = 2, by = x -> (x[NR], x[ME]))
    sendelems = Int[]
    counts = zeros(Int, csize + 1)
    counts[1] = (last(size(B)) > 0) ? 1 : 0
    sr, se = -1, 0
    for i in 1:last(size(B))
        r, e = B[NR, i], B[ME, i]
        # See if we need to send element `e` to rank `r` and make sure that we
        # didn't already mark it for sending.
        if r != crank && !(sr == r && se == e)
            counts[r + 2] += 1
            append!(sendelems, e)
            sr, se = r, e
        end
    end

    # Mark which faces need to be sent
    sendfaces = BitArray(undef, nface, length(sendelems))
    sendfaces .= false
    sr, se, n = -1, 0, 0
    for i in 1:last(size(B))
        r, e, f = B[NR, i], B[ME, i], B[MF, i]
        if r != crank
            if !(sr == r && se == e)
                n += 1
                sr, se = r, e
            end
            sendfaces[f, n] = true
        end
    end

    sendstarts = cumsum(counts)
    nabrtosendrank = Int[
        r for r = 0:(csize - 1) if sendstarts[r + 2] - sendstarts[r + 1] > 0
    ]
    nabrtosend = UnitRange{Int}[
        (sendstarts[r + 1]:(sendstarts[r + 2] - 1))
        for r = 0:(csize - 1) if sendstarts[r + 2] - sendstarts[r + 1] > 0
    ]

    # get element receiving information
    B = sortslices(B, dims = 2, by = x -> (x[NR], x[NE]))
    counts = zeros(Int, csize + 1)
    counts[1] = (last(size(B)) > 0) ? 1 : 0
    sr, se = -1, 0
    nghost = 0
    for i in 1:last(size(B))
        r, e = B[NR, i], B[NE, i]
        if r != crank
            # Check to make sure we have not already marked the element for
            # receiving since we could be connected to the receiving element across
            # multiple faces.
            if !(sr == r && se == e)
                nghost += 1
                counts[r + 2] += 1
                sr, se = r, e
            end
            B[NE, i] = nelem + nghost
        end
    end

    # Mark which faces will be received
    ghostfaces = BitArray(undef, nface, nghost)
    ghostfaces .= false
    sr, se, ge = -1, 0, 0
    for i in 1:last(size(B))
        r, e, f = B[NR, i], B[NE, i], B[NF, i]
        if r != crank
            if !(sr == r && se == e)
                ge += 1
                sr, se = r, e
            end
            B[NR, i] = crank
            ghostfaces[f, ge] = true
        end
    end

    recvstarts = cumsum(counts)
    nabrtorecvrank = Int[
        r for r = 0:(csize - 1) if recvstarts[r + 2] - recvstarts[r + 1] > 0
    ]
    nabrtorecv = UnitRange{Int}[
        (recvstarts[r + 1]:(recvstarts[r + 2] - 1))
        for r = 0:(csize - 1) if recvstarts[r + 2] - recvstarts[r + 1] > 0
    ]

    @assert nabrtorecvrank == nabrtosendrank
    nabrtorank = nabrtorecvrank

    elemtoelem = repeat((1:(nelem + nghost))', nface, 1)
    elemtoface = repeat(1:nface, 1, nelem + nghost)
    elemtoordr = ones(Int, nface, nelem + nghost)

    if d == 2
        for i in 1:last(size(B))
            me, mf, mo = B[ME, i], B[MF, i], B[MO, i]
            ne, nf, no = B[NE, i], B[NF, i], B[NO, i]

            elemtoelem[mf, me] = ne
            elemtoface[mf, me] = nf
            elemtoordr[mf, me] = (no == mo ? 1 : 2)
        end
    else
        for i in 1:last(size(B))
            me, mf, mo = B[ME, i], B[MF, i], B[MO, i]
            ne, nf, no = B[NE, i], B[NF, i], B[NO, i]

            elemtoelem[mf, me] = ne
            elemtoface[mf, me] = nf
            if no != 1 || mo != 1
                error("TODO add support for other orientations")
            end
            elemtoordr[mf, me] = 1
        end
    end

    # fill the ghost values in elemtocoord
    newelemtocoord = similar(elemtocoord, coorddim, nvert, nelem + nghost)
    newelemtobndy = similar(elemtobndy, nface, nelem + nghost)

    sendelemtocoord = elemtocoord[:, :, sendelems]
    sendelemtobndy = elemtobndy[:, sendelems]

    crreq = [
        MPI.Irecv!(view(newelemtocoord, :, :, nelem .+ er), r, 666, comm)
        for (r, er) in zip(nabrtorank, nabrtorecv)
    ]

    brreq = [
        MPI.Irecv!(view(newelemtobndy, :, nelem .+ er), r, 666, comm)
        for (r, er) in zip(nabrtorank, nabrtorecv)
    ]

    csreq = [
        MPI.Isend(view(sendelemtocoord, :, :, es), r, 666, comm)
        for (r, es) in zip(nabrtorank, nabrtosend)
    ]

    bsreq = [
        MPI.Isend(view(sendelemtobndy, :, es), r, 666, comm)
        for (r, es) in zip(nabrtorank, nabrtosend)
    ]

    newelemtocoord[:, :, 1:nelem] .= elemtocoord
    newelemtobndy[:, 1:nelem] .= elemtobndy

    MPI.Waitall!([csreq; crreq; bsreq; brreq])

    (
        elems = 1:(nelem + nghost),       # range of          element indices
        realelems = 1:nelem,            # range of real     element indices
        ghostelems = nelem .+ (1:nghost), # range of ghost    element indices
        ghostfaces = ghostfaces,
        sendelems = sendelems,          # array of send     element indices
        sendfaces = sendfaces,
        elemtocoord = newelemtocoord,   # element to vertex coordinates
        elemtovert = nothing,           # element to locally unique global vertex number (for direct stiffness summation)
        elemtoelem = elemtoelem,        # element to neighboring element
        elemtoface = elemtoface,        # element to neighboring element face
        elemtoordr = elemtoordr,        # element to neighboring element order
        elemtobndy = newelemtobndy,     # element to boundary number
        nabrtorank = nabrtorank,        # list of neighboring processes MPI ranks
        nabrtorecv = nabrtorecv,        # neighbor receive ranges into `ghostelems`
        nabrtosend = nabrtosend,
    )        # neighbor send ranges into `sendelems`
end

"""
    (bndytoelem, bndytoface) = enumerateboundaryfaces!(elemtoelem, elemtobndy, periodicity, boundary)

Update the `elemtoelem` array based on the boundary faces specified in
`elemtobndy`. Builds the `bndytoelem` and `bndytoface` tuples.
"""
function enumerateboundaryfaces!(elemtoelem, elemtobndy, periodicity, boundary)
    nb = 0
    for i in 1:length(periodicity)
        if !periodicity[i]
            nb = max(nb, boundary[i]...)
        end
    end
    # Limitation of the boundary condition unrolling in the DG kernels
    # (should never be violated unless more general unstructured meshes are used
    # since cube meshes only have 6 faces in 3D, and only 1 bcs is currently allowed
    # per face)

    @assert nb <= 6

    bndytoelem = ntuple(b -> Vector{Int64}(), nb)
    bndytoface = ntuple(b -> Vector{Int64}(), nb)

    nface, nelem = size(elemtoelem)

    N = zeros(Int, nb)
    for e in 1:nelem
        for f in 1:nface
            d = elemtobndy[f, e]
            @assert 0 <= d <= nb
            if d != 0
                elemtoelem[f, e] = N[d] += 1
                push!(bndytoelem[d], e)
                push!(bndytoface[d], f)
            end
        end
    end
    return (bndytoelem, bndytoface)
end

"""
    connectmeshfull(comm::MPI.Comm, elemtovert, elemtocoord, elemtobndy,
                faceconnections)

This function takes in a mesh (as returned for example by `brickmesh`) and
returns a corner connected mesh. It is similar to [`connectmesh`](@ref) but returns
a corner connected mesh, i.e. `ghostelems` will also include remote elements
that are only connected by vertices. This returns a `NamedTuple` of:

 - `elems` the range of element indices
 - `realelems` the range of real (aka nonghost) element indices
 - `ghostelems` the range of ghost element indices
 - `ghostfaces` ghost element to face is received;
   `ghostfaces[f,ge] == true` if face `f` of ghost element `ge` is received.
 - `sendelems` an array of send element indices
 - `sendfaces` send element to face is sent;
   `sendfaces[f,se] == true` if face `f` of send element `se` is sent.
 - `elemtocoord` element to vertex coordinates; `elemtocoord[d,i,e]` is the
    `d`th coordinate of corner `i` of element `e`
 - `elemtoelem` element to neighboring element; `elemtoelem[f,e]` is the
   number of the element neighboring element `e` across face `f`.  If there is
   no neighboring element then `elemtoelem[f,e] == e`.
 - `elemtoface` element to neighboring element face; `elemtoface[f,e]` is the
   face number of the element neighboring element `e` across face `f`.  If
   there is no neighboring element then `elemtoface[f,e] == f`.
 - `elemtoordr` element to neighboring element order; `elemtoordr[f,e]` is the
   ordering number of the element neighboring element `e` across face `f`.  If
   there is no neighboring element then `elemtoordr[f,e] == 1`.
 - `elemtobndy` element to bounday number; `elemtobndy[f,e]` is the
   boundary number of face `f` of element `e`.  If there is a neighboring
   element then `elemtobndy[f,e] == 0`.
 - `nabrtorank` a list of the MPI ranks for the neighboring processes
 - `nabrtorecv` a range in ghost elements to receive for each neighbor
 - `nabrtosend` a range in `sendelems` to send for each neighbor

The algorithm assumes that the 2-D mesh can be gathered on single process,
which is reasonable for the 2-D mesh.
"""
function connectmeshfull(
    comm::MPI.Comm,
    elemtovert,
    elemtocoord,
    elemtobndy,
    faceconnections;
    dim = size(elemtocoord, 1),
)
    @assert dim == 2
    dim_coord = size(elemtocoord, 1)
    csize = MPI.Comm_size(comm)
    crank = MPI.Comm_rank(comm)
    root = 0
    I = eltype(elemtovert)
    FT = eltype(elemtocoord)
    nvert, nelem = size(elemtovert)
    nfaces = 2 * dim
    nelemv = zeros(I, csize)
    fmask = build_fmask(dim)
    nfvert = size(fmask, 1)
    # collecting # of realelems on each process
    MPI.Allgather!(nelem, UBuffer(nelemv, Cint(1)), comm)
    offset = cumsum(nelemv) .+ 1
    pushfirst!(offset, 1)
    nelemg = sum(nelemv)
    # collecting elemtovert on each process
    elemtovertg = zeros(I, nvert + 2, nelemg)
    MPI.Allgatherv!(
        [elemtovert; ones(I, 1, nelem) .* crank; reshape(Array(1:nelem), 1, :)],
        VBuffer(elemtovertg, nelemv .* (nvert + 2)),
        comm,
    )

    nvertg = maximum(elemtovertg[1:nvert, :])
    # collecting elemtocoordg on each process
    elemtocoordg = Array{FT}(undef, dim_coord, nvert, nelemg)
    MPI.Allgatherv!(
        elemtocoord,
        VBuffer(elemtocoordg, nelemv .* (dim_coord * nvert)),
        comm,
    )
    # collecting elemtobndy on each process
    elemtobndyg = zeros(I, nfaces, nelemg)
    MPI.Allgatherv!(elemtobndy, VBuffer(elemtobndyg, nelemv .* nfaces), comm)
    # accounting for vertices on periodic faces
    nper = length(faceconnections) * nfvert
    vconn = Array{Int}(undef, 2, nper)          # stores the periodic and corresponding shadow vertex
    ctr = 1
    for fc in faceconnections
        e, f, v = fc[1], fc[2], fc[3:end]
        fv = elemtovert[fmask[:, f], e]
        fv, o = vertsortandorder(fv)
        v, o = vertsortandorder(v)
        for i in 1:nfvert
            vconn[1, ctr] = fv[1][i]
            vconn[2, ctr] = v[1][i]
            ctr += 1
        end
    end
    vconn = unique(vconn, dims = 2)
    nper = size(vconn, 2)
    nperv = zeros(I, csize)
    MPI.Allgather!(nper, UBuffer(nperv, Cint(1)), comm)
    nperg = sum(nperv)
    vconng = Array{Int}(undef, 2, nperg)
    MPI.Allgatherv!(vconn, VBuffer(vconng, nperv .* 2), comm)

    gldofv = -ones(Int, nvertg)
    pmarker = -ones(Int, nperg)
    # conflict-free periodic nodes
    for i in 1:nperg
        v1, v2 = vconng[1, i], vconng[2, i]
        if gldofv[v1] == -1 && gldofv[v2] == -1
            id = min(v1, v2)
            gldofv[v1] = id
            gldofv[v2] = id
            pmarker[i] = 1
        end
    end
    # dealing with doubly periodic nodes (whenever applicable)
    for i in 1:nperg
        if pmarker[i] == -1
            v1, v2 = vconng[1, i], vconng[2, i]
            id = min(gldofv[v1], gldofv[v2])
            gldofv[v1] = id
            gldofv[v2] = id
        end
    end
    # labeling non-periodic vertices
    for i in 1:nvertg
        if gldofv[i] == -1
            gldofv[i] = i
        end
    end
    #-------------------------------------------------------------
    elemtovertg_orig = deepcopy(elemtovertg)
    for i in 1:nelemg, j in 1:nvert
        elemtovertg[j, i] = gldofv[elemtovertg[j, i]]
    end
    #-------------------------------------------------
    nsend, nghost = 0, 0
    interiorelems, exteriorelems = Int[], Int[]
    vertgtoprocs = map(i -> zeros(Int, i), zeros(Int, nvertg)) # procs for each vertex
    vertgtolelem = map(i -> zeros(Int, i), zeros(Int, nvertg)) # local cell # for each vertex
    vertgtolvert = map(i -> zeros(Int, i), zeros(Int, nvertg)) # local vertex # in local cell for each vertex
    sendelems = map(i -> zeros(Int, i), zeros(Int, csize))
    recvelems = map(i -> zeros(Int, i), zeros(Int, csize))

    for icls in 1:nelemg # gathering connectivity info for each global vertex
        for ivt in 1:nvert
            gvt = elemtovertg[ivt, icls]
            prc = elemtovertg[nvert + 1, icls]
            lcl = elemtovertg[nvert + 2, icls]
            push!(vertgtoprocs[gvt], prc)
            push!(vertgtolelem[gvt], lcl)
            push!(vertgtolvert[gvt], ivt)
        end
    end

    for icls in 1:nelem # building sendelems and recvelems
        interior = true
        for ivt in 1:nvert
            vt = elemtovert[ivt, icls]
            gvt = gldofv[vt]
            for ip in 1:length(vertgtoprocs[gvt])
                proc = vertgtoprocs[gvt][ip]
                if proc ≠ crank
                    lcell = vertgtolelem[gvt][ip]
                    if findfirst(recvelems[proc + 1] .== lcell) == nothing
                        push!(recvelems[proc + 1], lcell)
                        nghost += 1
                    end
                    if findfirst(sendelems[proc + 1] .== icls) == nothing
                        push!(sendelems[proc + 1], icls)
                        nsend += 1
                    end
                    interior = false
                end
            end
        end
        interior ? push!(interiorelems, icls) : push!(exteriorelems, icls)
    end

    nabrtorank = Int[]
    newsendelems = Array{Int}(undef, nsend)
    nabrtosend = Array{UnitRange{Int64}}(undef, 0)
    nabrtorecv = Array{UnitRange{Int64}}(undef, 0)
    st_recv, en_recv = 1, 1
    st_send, en_send = 1, 1

    for ipr in 0:(csize - 1)
        l_send = length(sendelems[ipr + 1])
        l_recv = length(recvelems[ipr + 1])
        if l_send > 0
            en_send = st_send + l_send - 1
            sort!(sendelems[ipr + 1])
            newsendelems[st_send:en_send] .= sendelems[ipr + 1][:]
            push!(nabrtosend, st_send:en_send)
            st_send = en_send + 1
            push!(nabrtorank, ipr)

        end
        if l_recv > 0
            en_recv = st_recv + l_recv - 1
            sort!(recvelems[ipr + 1])
            push!(nabrtorecv, st_recv:en_recv)
            st_recv = en_recv + 1
        end
    end

    sendfaces = BitArray(zeros(nfaces, nsend))
    ghostfaces = BitArray(zeros(nfaces, nghost))

    newelemtovert = similar(elemtovert, nvert, (nelem + nghost))
    newelemtocoord = similar(elemtocoord, dim_coord, nvert, (nelem + nghost))
    newelemtobndy = similar(elemtobndy, nfaces, (nelem + nghost))
    newelemtovert[1:nvert, 1:nelem] .= elemtovert
    newelemtocoord[:, :, 1:nelem] .= elemtocoord
    newelemtobndy[:, 1:nelem] .= elemtobndy
    vmarker = BitArray(undef, nvert)
    ctrg, ctrs = 1, 1

    for ipr in 0:(csize - 1)
        # building ghost faces
        for icls in recvelems[ipr + 1]
            vmarker .= 0
            off = offset[ipr + 1]
            newelemtovert[1:nvert, nelem + ctrg] .=
                elemtovertg_orig[1:nvert, off + icls - 1]
            newelemtocoord[:, :, nelem + ctrg] .=
                elemtocoordg[:, :, off + icls - 1]
            newelemtobndy[:, nelem + ctrg] .= elemtobndyg[:, off + icls - 1]
            for ivt in 1:nvert
                gvt = elemtovertg[ivt, icls + off - 1]
                if findfirst(vertgtoprocs[gvt] .== crank) ≠ nothing
                    vmarker[ivt] = 1
                end
            end
            for fc in 1:nfaces
                if findfirst(vmarker[fmask[:, fc]]) ≠ nothing
                    ghostfaces[fc, ctrg] = 1
                end
            end
            ctrg += 1
        end
        # building send faces
        for icls in sendelems[ipr + 1]
            vmarker .= 0
            for ivt in 1:nvert
                vt = elemtovert[ivt, icls]
                gvt = gldofv[vt]
                if findfirst(vertgtoprocs[gvt] .== ipr) ≠ nothing
                    vmarker[ivt] = 1
                end
            end
            for fc in 1:nfaces
                if findfirst(vmarker[fmask[:, fc]]) ≠ nothing
                    sendfaces[fc, ctrs] = 1
                end
            end
            ctrs += 1
        end
    end

    A = zeros(Int, (nfvert + 3), nfaces * (nelem + nghost))
    for e in 1:(nelem + nghost)
        v = reshape(newelemtovert[:, e], ntuple(j -> 2, dim))
        for f in 1:nfaces
            j = (e - 1) * nfaces + f
            fv, o = vertsortandorder(v[fmask[:, f]]...) # faces vertices, B-N orientation
            for fvt in 1:nfvert
                A[fvt, j] = gldofv[fv[fvt]]
            end
            A[nfvert + 1, j] = o # orientation
            A[nfvert + 2, j] = e # local element number
            A[nfvert + 3, j] = f # local face number for element e
        end
    end

    A = sortslices(A, dims = 2, by = x -> x[1:nfvert])
    elemtoelem = Array{Int}(undef, nfaces, (nelem + nghost))
    elemtoface = Array{Int}(undef, nfaces, (nelem + nghost))
    elemtoordr = Array{Int}(undef, nfaces, (nelem + nghost))
    # match faces
    n = size(A, 2)
    j = 1
    while j ≤ n
        lel = A[nfvert + 2, j]
        lfc = A[nfvert + 3, j]
        if j + 1 ≤ n && A[1:nfvert, j] == A[1:nfvert, j + 1]
            # found connected face
            nel = A[nfvert + 2, j + 1] # local neighboring element
            nfc = A[nfvert + 3, j + 1] # local face # of local neighboring element
            elemtoelem[lfc, lel] = nel
            elemtoface[lfc, lel] = nfc
            elemtoelem[nfc, nel] = lel
            elemtoface[nfc, nel] = lfc
            if A[nfvert + 1, j] == A[nfvert + 1, j + 1]
                elemtoordr[lfc, lel] = 1
                elemtoordr[nfc, nel] = 1
            else
                elemtoordr[lfc, lel] = 2
                elemtoordr[nfc, nel] = 2
            end
            j += 2
        else
            # found unconnected face
            elemtoelem[lfc, lel] = lel
            elemtoface[lfc, lel] = lfc
            elemtoordr[lfc, lel] = 1
            j += 1
        end
    end
    # provide locally unique elemtovert for DSS
    uvert = -ones(Int, nvertg)
    for el in 1:last(size(newelemtovert)), i in 1:nvert
        newelemtovert[i, el] = gldofv[newelemtovert[i, el]]
        uvert[newelemtovert[i, el]] = 0
    end
    ctr = 1
    for i in 1:length(uvert)
        if uvert[i] == 0
            uvert[i] = ctr
            ctr += 1
        end
    end
    for el in 1:last(size(newelemtovert)), i in 1:nvert
        newelemtovert[i, el] = uvert[newelemtovert[i, el]]
    end
    #---------------------------------------------------------------
    (
        elems = 1:(nelem + nghost),       # range of element indices
        realelems = 1:nelem,              # range of real element indices
        ghostelems = nelem .+ (1:nghost), # range of ghost element indices
        ghostfaces = ghostfaces,          # bit array of recv faces for ghost elems
        sendelems = newsendelems,         # array of send element indices
        sendfaces = sendfaces,            # bit array of send faces for send elems
        elemtocoord = newelemtocoord,     # element to vertex coordinates
        elemtovert = newelemtovert,       # element to locally unique global vertex number (for direct stiffness summation)
        elemtoelem = elemtoelem,          # element to neighboring element
        elemtoface = elemtoface,          # element to neighboring element face
        elemtoordr = elemtoordr,          # element to neighboring element order
        elemtobndy = newelemtobndy,       # element to boundary number
        nabrtorank = nabrtorank,          # list of neighboring processes MPI ranks
        nabrtorecv = nabrtorecv,          # neighbor receive ranges into `ghostelems`
        nabrtosend = nabrtosend,          # neighbor send ranges into `sendelems`
    )
end

"""
    build_fmask(dim)

Returns the face mask for mapping element vertices to face vertices.

ex: for 2D element with vertices (1, 2, 3, 4)

       3---4
       |   |
       1---2

the function returns the face mask
f1 | f2 | f3 | f4
=================
 1 |  2 |  1 | 3
 3 |  4 |  2 | 4
=================
"""
function build_fmask(dim)
    nvert = 2^dim
    nfaces = 2 * dim
    p = reshape(1:nvert, ntuple(j -> 2, dim))
    fmask = Array{Int64}(undef, 2^(dim - 1), nfaces)
    f = 0
    for d in 1:dim
        for slice in eachslice(p, dims = d)
            fmask[:, f += 1] = vec(slice)
        end
    end
    return fmask
end

end # module


text:
@static if VERSION < v"0.6"
    ⊻ = $
end

const global IS_LITTLE_ENDIAN = (Base.ENDIAN_BOM == 0x04030201)
const global IS_64_BIT = (Sys.WORD_SIZE == 64)

if IS_LITTLE_ENDIAN && IS_64_BIT
    include("fastmemory.jl")
else
    include("safememory.jl")
end

# byte patterns for each of the tag types in the compression format.
const global SNAPPY_LITERAL = 0x00
const global SNAPPY_COPY_1_BYTE_OFFSET = 0x01
const global SNAPPY_COPY_2_BYTE_OFFSET = 0x02
const global SNAPPY_COPY_4_BYTE_OFFSET = 0x03

# : snappy.h
# the size of a compression block. note that many parts of the compression
# code assumes that K_BLOCK_SIZE <= 65536; in particular, the hash table
# can only store 16-bit offsets, and emit_copy! also assumes the offset
# is 65535 bytes or less.
#
# note that there might be older data around that is compressed with larger
# block sizes, so the decompression code should not rely on the
# non-existence of long backreferences.
const global K_BLOCK_SIZE = 1 << 16
const global K_INPUT_MARGIN_BYTES = 15
const global K_MAX_HASH_TABLE_SIZE = 1 << 14

# : snappy-internal.h
# data stored per entry in lookup table:
#      range   bits-used       description
#      ------------------------------------
#      1..64   0..7            literal/copy length encoded in opcode byte
#      0..7    8..10           copy offset encoded in opcode byte / 256
#      0..4    11..13          extra bytes after opcode
#
# we use eight bits for the length even though 7 would have sufficed
# because of efficiency reasons:
#      (1) extracting a byte is faster than a bit-field
#      (2) it properly aligns copy offset so we do not need a <<8
const global CHAR_TABLE = UInt16[
    0x0001, 0x0804, 0x1001, 0x2001, 0x0002, 0x0805, 0x1002, 0x2002,
    0x0003, 0x0806, 0x1003, 0x2003, 0x0004, 0x0807, 0x1004, 0x2004,
    0x0005, 0x0808, 0x1005, 0x2005, 0x0006, 0x0809, 0x1006, 0x2006,
    0x0007, 0x080a, 0x1007, 0x2007, 0x0008, 0x080b, 0x1008, 0x2008,
    0x0009, 0x0904, 0x1009, 0x2009, 0x000a, 0x0905, 0x100a, 0x200a,
    0x000b, 0x0906, 0x100b, 0x200b, 0x000c, 0x0907, 0x100c, 0x200c,
    0x000d, 0x0908, 0x100d, 0x200d, 0x000e, 0x0909, 0x100e, 0x200e,
    0x000f, 0x090a, 0x100f, 0x200f, 0x0010, 0x090b, 0x1010, 0x2010,
    0x0011, 0x0a04, 0x1011, 0x2011, 0x0012, 0x0a05, 0x1012, 0x2012,
    0x0013, 0x0a06, 0x1013, 0x2013, 0x0014, 0x0a07, 0x1014, 0x2014,
    0x0015, 0x0a08, 0x1015, 0x2015, 0x0016, 0x0a09, 0x1016, 0x2016,
    0x0017, 0x0a0a, 0x1017, 0x2017, 0x0018, 0x0a0b, 0x1018, 0x2018,
    0x0019, 0x0b04, 0x1019, 0x2019, 0x001a, 0x0b05, 0x101a, 0x201a,
    0x001b, 0x0b06, 0x101b, 0x201b, 0x001c, 0x0b07, 0x101c, 0x201c,
    0x001d, 0x0b08, 0x101d, 0x201d, 0x001e, 0x0b09, 0x101e, 0x201e,
    0x001f, 0x0b0a, 0x101f, 0x201f, 0x0020, 0x0b0b, 0x1020, 0x2020,
    0x0021, 0x0c04, 0x1021, 0x2021, 0x0022, 0x0c05, 0x1022, 0x2022,
    0x0023, 0x0c06, 0x1023, 0x2023, 0x0024, 0x0c07, 0x1024, 0x2024,
    0x0025, 0x0c08, 0x1025, 0x2025, 0x0026, 0x0c09, 0x1026, 0x2026,
    0x0027, 0x0c0a, 0x1027, 0x2027, 0x0028, 0x0c0b, 0x1028, 0x2028,
    0x0029, 0x0d04, 0x1029, 0x2029, 0x002a, 0x0d05, 0x102a, 0x202a,
    0x002b, 0x0d06, 0x102b, 0x202b, 0x002c, 0x0d07, 0x102c, 0x202c,
    0x002d, 0x0d08, 0x102d, 0x202d, 0x002e, 0x0d09, 0x102e, 0x202e,
    0x002f, 0x0d0a, 0x102f, 0x202f, 0x0030, 0x0d0b, 0x1030, 0x2030,
    0x0031, 0x0e04, 0x1031, 0x2031, 0x0032, 0x0e05, 0x1032, 0x2032,
    0x0033, 0x0e06, 0x1033, 0x2033, 0x0034, 0x0e07, 0x1034, 0x2034,
    0x0035, 0x0e08, 0x1035, 0x2035, 0x0036, 0x0e09, 0x1036, 0x2036,
    0x0037, 0x0e0a, 0x1037, 0x2037, 0x0038, 0x0e0b, 0x1038, 0x2038,
    0x0039, 0x0f04, 0x1039, 0x2039, 0x003a, 0x0f05, 0x103a, 0x203a,
    0x003b, 0x0f06, 0x103b, 0x203b, 0x003c, 0x0f07, 0x103c, 0x203c,
    0x0801, 0x0f08, 0x103d, 0x203d, 0x1001, 0x0f09, 0x103e, 0x203e,
    0x1801, 0x0f0a, 0x103f, 0x203f, 0x2001, 0x0f0b, 0x1040, 0x2040,
]

# mapping from i in range [0,4] to a mask to extract the bottom 8*i bits.
const global WORDMASK = UInt32[
    0x00000000, 0x000000ff, 0x0000ffff, 0x00ffffff, 0xffffffff
]

# : snappy.cc
# any hash function will produce a valid compressed bitstream, but a good
# hash function reduces the number of collisions and thus yields better
# compression for compressible input, and more speed for incompressible
# input. Of course, it doesn't hurt if the hash function is reasonably fast
# either, as it gets called a lot.
@inline hashdword(bytes::UInt32, shift::Unsigned) = (bytes * 0x1e35a7bd) >> shift

# : snappy.cc
# leturn floor(log2(n)) for positive integer n. returns -1 iff n == 0.
@inline log2floor(n::UInt32) = n == 0 ? -1 : 31 ⊻ leading_zeros(n)

# : snappy.cc
# use smaller hash table when length(input) is smaller, since we
# fill the table, incurring O(hash table size) overhead for
# compression, and if the input is short, we won't need that
# many hash table entries anyway.
function alloc_hashtable(n::Unsigned)
    htsize = 256
    while htsize < K_MAX_HASH_TABLE_SIZE && htsize < n
        htsize <<= 0x01
    end
    return zeros(UInt16, htsize)
end

# : snappy.cc
# flat array compression that does not emit the "uncompressed length"
# prefix. compresses "input" bytes into "output".
#
# REQUIRES: "input" is at most "K_BLOCK_SIZE" bytes long.
# REQUIRES: "output" array has been allocated with at least
# "maxlength_compressed(length(input))" bytes of memory.
# REQUIRES: All elements in "table" are initialized to 0xffff.
# REQUIRES: "length(table)" is a power of two.
#
# returns an index into output one past the end of the compressed
# "input" bytes. The compressed length of "input" is outputindex-1.
function compress_fragment!(output::Vector{UInt8}, outputindex::Integer, input::Vector{UInt8}, ip::Integer, ip_end::Integer, table::Vector{UInt16})
    local shift::UInt = 32 - log2floor(convert(UInt32,length(table)))
    candidate = next_emit = base_ip = ip
    input_size = ip_end-ip+1
    ip_limit = ip_end - K_INPUT_MARGIN_BYTES

    @inbounds if input_size >= K_INPUT_MARGIN_BYTES
        while true

            # : snappy.cc
            # the body of this loop calls emit_literal! once and then emit_copy! one or
            # more times (the exception is that when we're close to exhausting
            # the input we @goto emit_remainder).
            #
            # in the first iteration of this loop we're just starting, so
            # there's nothing to copy, so calling emit_literal! once is
            # necessary. and we only start a new iteration when the
            # current iteration has determined that a call to emit_literal! will
            # precede the next call to emit_copy! (if any).
            #
            # step 1: scan forward in the input looking for a 4-byte-long match.
            # if we get close to exhausting the input then @goto emit_remainder.
            #
            # heuristic match skipping: if 32 bytes are scanned with no matches
            # found, start looking only at every other byte. if 32 more bytes are
            # scanned (or skipped), look at every third byte, etc.. When a match is
            # found, immediately go back to looking at every byte. this is a small
            # loss (~5% performance, ~0.1% density) for compressible data due to more
            # bookkeeping, but for non-compressible data (such as JPEG) it's a huge
            # win since the compressor quickly "realizes" the data is incompressible
            # and doesn't bother looking for matches everywhere.
            #
            # the "skip" variable keeps track of how many bytes there are since the
            # last match; dividing it by 32 (ie. right-shifting by five) gives the
            # number of bytes to move ahead for each iteration.
            skip = 32
            next_hash = hashdword(load32u(input, ip+=1), shift)
            next_ip = ip

            # this loop is a do {...} while (load32u(input, candidate) == load32u(input,ip))
            while true
                ip = next_ip
                cur_hash = next_hash
                bytes_between_hash_lookups = skip >> 0x05
                skip += bytes_between_hash_lookups
                next_ip = ip + bytes_between_hash_lookups

                # early exit if no match is found and we're nearing the end of the input
                (next_ip > ip_limit) && @goto emit_remainder

                # the table needs to hold indices up to 65536, but is constrained to a UInt16.
                # so we subtract 1 when storing, and add 1 when retrieving. the trick is that
                # +1 should only happen when retrieving values that have already been stored,
                # that is to say, the default value needs to be 0. so rather than filling the
                # array with 0's, fill it with 0xffff instead, and force the overflow back to 0.
                # in this case the addition is not associative, the +1 needs to evaluate on the
                # table first to force the overflow, and only then then add to base_ip.
                #
                # it is.. unfortunate that array indices are 1-based. the rest of the code can likely
                # be refactored and made more better, but as most of it is a ported from C, it is
                # easier to handle it this way for now. TODO (krm01): consider using the neat new
                # feature that allows arrays with arbitrary indexes and start it from 0.
                next_hash = hashdword(load32u(input, next_ip), shift)
                candidate = base_ip + (table[cur_hash+1] + 0x0001)
                table[cur_hash+1] = (ip - base_ip) - 1

                (load32u(input, candidate) == load32u(input,ip)) && break
            end

            # : snappy.cc
            # step 2: a 4-byte match has been found. we'll later see if more
            # than 4 bytes match. but, prior to the match, input
            # bytes [next_emit, ip) are unmatched. emit them as "literal bytes."
            outputindex = emit_literal!(output, outputindex, input, next_emit, ip-next_emit, true)

            # : snappy.cc
            # step 3: call emit_copy!, and then see if another emit_copy! could
            # be our next move. repeat until we find no match for the
            # input immediately after what was consumed by the last emit_copy! call.
            #
            # if we exit this loop normally then we need to call emitLiteral next,
            # though we don't yet know how big the literal will be. we handle that
            # by proceeding to the next iteration of the main loop. we also can exit
            # this loop via goto if we get close to exhausting the input.
            while true

                # : snappy.cc
                # we have a 4-byte match at ip, and no need to emit any
                # "literal bytes" prior to ip.
                matched = 4 + find_match_length(input, candidate+4, ip+4, ip_end)
                outputindex = emit_copy!(output, outputindex, ip - candidate, matched)
                ip += matched
                next_emit = ip

                # early exit if no match is found and we're nearing the end of the input
                (ip >= ip_limit) && @goto emit_remainder

                # : snappy.cc
                # we are now looking for a 4-byte match again. we read
                # table[hashdword(load32u(input, ip), shift)] for that.
                # to improve compression, we also update the table at ip-1 and ip.
                prev_hash = hashdword(load32u(input, ip-1), shift)
                input_bytes = load32u(input, ip)
                cur_hash = hashdword(input_bytes, shift)

                # we need to do the same awkward +0x0001, -1 trick here.
                table[prev_hash+1] = ip - base_ip - 1 - 1
                candidate = base_ip + (table[cur_hash+1] + 0x0001)
                table[cur_hash+1] = ip - base_ip - 1

                # another do {...} while () break
                (input_bytes != load32u(input, candidate)) && break
            end
        end
    end
    @label emit_remainder
    # emit the remaining bytes as a literal
    if next_emit <= ip_end
        # we are near the end of the input, so can't allow emit_literal! to
        # write beyond its bounds. calling with false indicates to use the slow/safe path.
        outputindex = emit_literal!(output, outputindex, input, next_emit, ip_end-next_emit+1, false)
    end
    return outputindex
end

@inline function emit_literal!(output::Vector{UInt8}, outputindex::Integer, input::Vector{UInt8}, inputindex::Integer, len::Integer, allow_fast_path::Bool)
    # : snappy.cc
    # the vast majority of copies are below 16 bytes, for which a
    # call to memcpy is overkill. this fast path can sometimes
    # copy up to 15 bytes too much, but that is okay in the
    # main loop, since we have a bit to go on for both sides:
    #
    #   - the input will always have K_INPUT_MARGIN_BYTES = 15 extra
    #     available bytes, as long as we're in the main loop, and
    #     if not, allow_fast_path = false.
    #   - the output will always have 32 spare bytes (see
    #     maxlength_compressed).
    local n::UInt32 = (len - 1) % UInt32
    if (allow_fast_path && len <= 16)
        # fits in tag byte
        output[outputindex] = SNAPPY_LITERAL | ((n << 2) % UInt8)
        unaligned_copy_128u!(output, outputindex+=1, input, inputindex)
        return outputindex + len
    end
    if len < 60
        # fits in tag byte
        fb = SNAPPY_LITERAL | ((n << 2) % UInt8)
        output[outputindex] = fb
    else
        # encode in upcoming bytes
        count = 0
        base = outputindex
        while n > 0
            output[outputindex+=1] = n % UInt8
            n >>= 8; count += 1;
        end
        output[base] = SNAPPY_LITERAL | (((59+count) << 2) % UInt8)
    end
    copy!(output, outputindex+=1, input, inputindex, len)
    return outputindex + len
end

@inline function emit_copy_upto_64!(output::Vector{UInt8}, outputindex::Integer, offset::Integer, len::Integer)
    if len < 12 && offset < 2048
        # offset fits in 11 bits. the 3 highest go in the top of the first byte,
        # and the rest go in the second byte.
        output[outputindex] = (SNAPPY_COPY_1_BYTE_OFFSET + ((len - 4) << 2) + ((offset >> 3) & 0xe0)) % UInt8
        output[outputindex+1] = (offset & 0xff) % UInt8
        outputindex += 2
    else
        # write 4 bytes, though we only care about 3 of them. The output buffer
        # is required to have some slack, so the extra byte won't overrun it.
        local u::UInt32 = SNAPPY_COPY_2_BYTE_OFFSET + ((len - 1) << 2) + (offset << 8)
        store32u!(output, outputindex, u)
        outputindex += 3
    end
    return outputindex
end

@inline function emit_copy!(output::Vector{UInt8}, outputindex::Integer, offset::Integer, len::Integer)
    if len < 12
        outputindex = emit_copy_upto_64!(output, outputindex, offset, len)
        return outputindex
    else
        # : snappy.cc
        # a special case for len <= 64 might help, but so far measurements suggest
        # it's in the noise. TODO (krm01): confirm this is still true in Julia

        # emit 64 byte copies but make sure to keep at least four bytes reserved.
        while len >= 68
            outputindex = emit_copy_upto_64!(output, outputindex, offset, 64)
            len -= 64
        end
        # one or two copies will now finish the job.
        if len > 64
            outputindex = emit_copy_upto_64!(output, outputindex, offset, 60)
            len -= 60
        end
        # emit remainder.
        outputindex = emit_copy_upto_64!(output, outputindex, offset, len)
    end
    return outputindex
end

# : snappy-internal.h
# find the largest n such that
#
#   a[i1, n-1] == a[i2, n-1]
#   and n <= (limit - i2).
#
# does not read a[limit] or beyond.
# does not read a[(i1 + (limit - i2))] or beyond.
# requires that limit >= i2.
#
# separate implementation for 64-bit, little-endian cpus.
@static if IS_64_BIT && IS_LITTLE_ENDIAN
    @inline function find_match_length(a::Vector{UInt8}, i1::Integer, i2::Integer, limit::Integer)
        matched = 0

        # : snappy-internal.h
        # this block isn't necessary for correctness; we could just start looping
        # immediately. as an optimization though, it is useful. it creates some not
        # uncommon code paths that determine, without extra effort, whether the match
        # length is less than 8. in short, we are hoping to avoid a conditional
        # branch, and perhaps get better code layout from the C++ compiler.

        # it turns out, this is also true for Julia. also we check (limit - 7)
        # instead of (limit - 8) because 1-based arrays.
        if i2 <= limit - 7
            a1 = load64u(a, i1)
            a2 = load64u(a, i2)
            if a1 != a2
                return trailing_zeros(a1 ⊻ a2) >> 3
            else
                i2 += 8
                matched = 8
            end
        end

        # find out how long the match is. we loop over the data 64 bits at a
        # time until we find a 64-bit block that doesn't match; then we find
        # the first non-matching bit and use that to calculate the total
        # length of the match.
        while i2 <= limit - 7
            if load64u(a, i2) == load64u(a, i1+matched)
                i2 += 8
                matched += 8
            else
                x = load64u(a, i2) ⊻ load64u(a, i1+matched)
                matching_bits = trailing_zeros(x)
                matched += (matching_bits >> 3)
                return matched
            end
        end
        while i2 <= limit && a[i1+matched] == a[i2]
            i2 += 1
            matched += 1
        end
        return matched
    end
else
# 32bit version, based on above routine.
    @inline function find_match_length(a::Vector{UInt8}, i1::Integer, i2::Integer, limit::Integer)
        matched = 0
        while i2 <= limit - 4 && load32u(a, i2) == load32u(a, i1+matched)
            i2 += 4
            matched += 4
        end
        if IS_LITTLE_ENDIAN && i2 <= limit - 4
            x = load32u(a, i2) ⊻ load32u(a, i1+matched)
            matching_bits = trailing_zeros(x)
            matched += (matching_bits >> 3)
        else
            while i2 <= limit && a[i1+matched] == a[i2]
                i2 += 1
                matched += 1
            end
        end
        return matched
    end
end

function decompress_all_tags!(output::Vector{UInt8}, input::Vector{UInt8}, ip::Integer)
    ip_limit = endof(input)
    op = start(output)
    op_limit = endof(output)

    @inbounds while ip < ip_limit

        c = input[ip]
        ip += 1

        # load32u reads 4 bytes blindly from input, so when we're nearing the
        # end of the array it needs to be handled differently. the Java version
        # built up the trailer byte-by-byte, but my testing shows it is faster
        # in Julia to create a slice padded with 0 bytes and let it read over
        # into that buffer. it seems kind of silly, but...
        if ip+4 <= ip_limit
            tag = load32u(input, ip)
        else
            tag = load32u(vcat(input[ip:ip_limit], zeros(UInt8, 4)), 1)
        end

        # here we have read c, the tag byte marker, and need to extract the
        # upper six bits encoding the length. the CHAR_TABLE is a precomputed
        # lookup table. then advance ip 1 beyond the end of the tag.
        entry = CHAR_TABLE[c+1]
        len = entry & 0xff
        taglen = entry >> 11
        trailer = tag & WORDMASK[taglen+1]
        ip += taglen

        # : snappy.cc
        # ratio of iterations that have LITERAL vs non-LITERAL for different
        # inputs.
        #
        # input          LITERAL  NON_LITERAL
        # -----------------------------------
        # html|html4|cp   23%        77%
        # urls            36%        64%
        # jpg             47%        53%
        # pdf             19%        81%
        # txt[1-4]        25%        75%
        # pb              24%        76%
        # bin             24%        76%
        #
        # in the C++ version, this table was used to justify a compiler directive
        # to help with branch prediction, but no such feature exists in Julia that
        # i know of. still, copies are more common, so that case is checked first.
        if ((c & 0x3) % UInt8) != SNAPPY_LITERAL
            copy_offset = (entry & 0x700)
            op = incremental_copy!(output, op, copy_offset + trailer, len)
        else
            op, ip = copy_literal!(output, op, input, ip, len + trailer)
        end
    end
    return op
end

# snappy supports run-length encoding (RLE), so handling like this is necessary
# in cases where the input and output regions may overlap.
# for example:
#       ip ---v
#            [ab                      ]    copy length: 20
#       op -----^
# after the copy the buffer will look like:
#            [ababababababababababab  ]
@inline function incremental_copy_slow!(output::Vector{UInt8}, op::Integer, ip::Integer, len::Integer)
    for i in 0:len-1
        output[op+i] = output[ip+i]
    end
end

# : snappy.cc
# copy [src, src+(op_limit-op)) to [op, (op_limit-op)) but faster than
# IncrementalCopySlow. buf_limit is the address past the end of the writable
# region of the buffer.
#
# copy output[offset:len] to output[op:len], optimized for x64 little endian
# systems to be faster than a call to copy!() - also catches the end of the copy
# to be handled by incremental_copy_slow!
@inline function incremental_copy!(output::Vector{UInt8}, op::Integer, offset::Integer, len::Integer)
    avail_out = endof(output) - op + 1
    # this "offset - 0x01" trick is from the original C++, we could just
    # check op - start < offset, but this version is equivalent in all cases
    # except when offset == 0, a case we do not want to allow. subtracting by
    # 0x01 causes an offset of 0 to be cast to unsigned and wrap around, and
    # is then caught by this check.
    (op - start(output) <= offset - 0x01) && error("Invalid input: corrupt copy offset")
    if len <= 16 && offset >= 8 && avail_out >= 16
        # use two 64s instead of a 128, because op - offset can be < 16. A 128 copy might overwrite into op.
        unaligned_copy_64u!(output, op, output, op - offset)
        unaligned_copy_64u!(output, op + 8, output, op - offset + 8)
    else
        (avail_out < len) && error("Invalid input: corrupt copy length")
        incremental_copy_slow!(output, op, op - offset, len)
    end
    return op + len
end

# when copying literal bytes there is no concern for overlapping input
# and output regions, so take the fast path when possible and otherwise
# use copy!() , which is still faster than incremental_copy_slow!()
@inline function copy_literal!(output::Vector{UInt8}, op::Integer, input::Vector{UInt8}, ip::Integer, len::Integer)
    avail_out = endof(output) - op + 1
    avail_in = endof(input) - ip + 1
    (avail_out < len || avail_in < len) && error("Invalid input: corrupt literal")

    # fast path, can copy more than the literal but that's OK here. most literals are short.
    if len <= 16 && avail_out >= 16 && avail_in >= 16
        unaligned_copy_128u!(output, op, input, ip)
    else
        copy!(output, op, input, ip, len)
    end
    return op+len, ip+len
end


text:
#=
Chr support

Copyright 2017-2018 Gandalf Software, Inc., Scott P. Jones
Licensed under MIT License, see LICENSE.md
In part based on code for Char in Julia
=#

"""Default value for Str /Chr types"""
codeunit(::Type{<:CSE}) = UInt8
codeunit(::Type{<:Word_CSEs}) = UInt16
codeunit(::Type{<:Quad_CSEs}) = UInt32

"""Type of codeunits"""
codeunit(::Type{UniStr})    = UInt32 # Note, the real type could be UInt8, UInt16, or UInt32
codeunit(::Type{<:Str{C}}) where {C<:CSE} = codeunit(C)

codeunit(::MaybeSub{S}) where {S<:Str} = codeunit(S)
codeunit(::Type{<:MaybeSub{S}}) where {S<:Str} = codeunit(S)

bytoff(::Type{UInt8},  off) = off
bytoff(::Type{UInt16}, off) = off << 1
bytoff(::Type{UInt32}, off) = off << 2
chroff(::Type{UInt8},  off) = off
chroff(::Type{UInt16}, off) = off >>> 1
chroff(::Type{UInt32}, off) = off >>> 2

chrdiff(pnt::Ptr{T}, beg::Ptr{T}) where {T<:CodeUnitTypes} = Int(chroff(T, pnt - beg))

bytoff(pnt::Ptr{T}, off) where {T<:CodeUnitTypes} = pnt + bytoff(T, off)

typemax(::Type{ASCIIChr}) = ASCIIChr(0x7f)
typemax(::Type{UTF32Chr}) = UTF32Chr(0x10ffff)

basetype(::Type{Char})        = UInt32
basetype(::Type{T}) where {T<:CodeUnitTypes} = T

eltype(::Type{<:Str{BinaryCSE}}) = UInt8

eltype(::Type{<:Str{Text1CSE}})       = Text1Chr
eltype(::Type{<:Str{Text2CSE}})       = Text2Chr
eltype(::Type{<:Str{Text4CSE}})       = Text4Chr
eltype(::Type{<:Str{ASCIICSE}})       = ASCIIChr
eltype(::Type{<:Str{LatinCSE}})       = LatinChr
eltype(::Type{<:Str{_LatinCSE}})      = _LatinChr
eltype(::Type{<:Str{<:Union{UCS2CSE, _UCS2CSE}}}) = UCS2Chr
eltype(::Type{<:Str{<:Union{UTF8CSE, UTF16CSE, UTF32_CSEs}}}) = UTF32Chr

codepoint_size(::Type{T}) where {T<:Union{String,Str}} = sizeof(eltype(T))

get_codeunit(dat, pos) = codeunit(dat, pos)
get_codeunit(pnt::Ptr{<:CodeUnitTypes}, pos) = unsafe_load(pnt, pos)
get_codeunit(dat::AbstractVector{<:CodeUnitTypes}, pos) = dat[pos]
get_codeunit(str::MaybeSub{<:Str}, pos::Integer) = @preserve str get_codeunit(pointer(str), pos)

get_codeunit(dat) = get_codeunit(dat, 1)
get_codeunit(pnt::Ptr{<:CodeUnitTypes}) = unsafe_load(pnt)

codeunit(str::Str, pos::Integer) = get_codeunit(str, pos)

set_codeunit!(pnt::Ptr{<:CodeUnitTypes}, pos, ch) = unsafe_store!(pnt, ch, pos)
set_codeunit!(dat::AbstractVector{<:CodeUnitTypes}, pos, ch) = (dat[pos] = ch)
set_codeunit!(dat::String, pos, ch) = unsafe_store!(pointer(dat), pos, ch)

set_codeunit!(pnt::Ptr{<:CodeUnitTypes}, ch) = unsafe_store!(pnt, ch)
set_codeunit!(dat::AbstractVector{<:CodeUnitTypes}, ch) = (dat[1] = ch)
set_codeunit!(dat::String, ch) = set_codeunit!(dat, 1, ch)

convert(::Type{T}, v::S) where {T<:Integer, S<:Chr} = convert(T, codepoint(v))::T
convert(::Type{T}, v::Signed) where {T<:Chr} =
    (v >= 0 && is_valid(T, v%Unsigned)) ? convert(T, v%Unsigned) : codepoint_error(T, v)
convert(::Type{T}, v::Unsigned) where {CS,B,T<:Chr{CS,B}} =
    is_valid(T, v) ? Chr(CS, v%B) : codepoint_error(T, v)
convert(::Type{Char}, v::T) where {T<:Chr} = convert(Char, codepoint(v))
convert(::Type{T}, v::Char) where {T<:Chr} = convert(T, codepoint(v))::T

rem(x::Number, ::Type{<:Chr{CS,B}}) where {CS,B} = Chr(CS, x%B)
rem(x::Char, ::Type{T}) where {T<:Chr}   = x%UInt32%T
rem(x::Chr, ::Type{T}) where {T<:Chr}    = (x.v)%T
rem(x::Chr, ::Type{T}) where {T<:Char}   = (x.v)%T
rem(x::Chr, ::Type{T}) where {T<:Number} = (x.v)%T

(::Type{S})(v::T) where {S<:Union{UInt32, Int, UInt}, T<:Chr} = codepoint(v)%S
(::Type{Char})(v::Chr) = Char(codepoint(v))
(::Type{T})(v::Char) where {T<:Chr} = T(codepoint(v))

eltype(::Type{T}) where {T<:Chr} = T
size(cp::Chr, dim) = convert(Int, dim) < 1 ? boundserr(cp, dim) : 1
getindex(cp::Chr, i::Integer) = i == 1 ? cp : boundserr(cp, i)
getindex(cp::Chr, I::Integer...) = all(x -> x == 1, I) ? cp : boundserr(cp, I)

_uni_rng(m) = 0x00000:ifelse(m < 0xd800, m, m-0x800)
codepoint_rng(::Type{T}) where {T<:Chr} = _uni_rng(typemax(T)%UInt32)
codepoint_rng(::Type{Char}) = _uni_rng(0x10ffff)
codepoint_rng(::Type{Text2Chr}) = 0%UInt16:typemax(UInt16)
codepoint_rng(::Type{Text4Chr}) = 0%UInt32:typemax(UInt32)

codepoint_adj(::Type{T}, ch) where {T} = ifelse(ch < 0xd800, ch, ch+0x800)%T
codepoint_adj(::Type{T}, ch) where {T<:Union{Text2Chr,Text4Chr}} = ch%T

# returns a random valid Unicode scalar value in the correct range for the type of character
@static if V6_COMPAT
import Base.Random: rand!, rand, AbstractRNG
rand(r::AbstractRNG, ::Type{T}) where {T<:Chr} =
    codepoint_adj(T, rand(r, codepoint_rng(T)))
rand!(rng::AbstractRNG, A::AbstractArray, r::UnitRange{<:Chr}) =
    rand!(rng, A, Base.Random.RangeGenerator(r))
else
Random.rand(r::Random.AbstractRNG, ::Random.SamplerType{T}) where {T<:Chr} =
    codepoint_adj(T, rand(r, codepoint_rng(T)))
end

==(x::Chr, y::AbsChar) = codepoint(x) == codepoint(y)
==(x::AbsChar, y::Chr) = codepoint(x) == codepoint(y)
==(x::Chr, y::Chr)     = codepoint(x) == codepoint(y)

isless(x::Chr, y::AbsChar) = codepoint(x) < codepoint(y)
isless(x::AbsChar, y::Chr) = codepoint(x) < codepoint(y)
isless(x::Chr, y::Chr)     = codepoint(x) < codepoint(y)

# This is so that the hash is compatible with isless, but it's very inefficient
Base.hash(x::Chr, h::UInt) = hash(Char(x), h)


text:
# ---
# title: 77. Combinations
# id: problem77
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Backtracking
# link: <https://leetcode.com/problems/combinations/description/>
# hidden: true
# ---
#
# Given two integers _n_ and _k_ , return all possible combinations of _k_
# numbers out of 1 ... _n_.
#
# You may return the answer in **any order**.
#
#
#
# **Example 1:**
#
#
#
#     Input: n = 4, k = 2
#     Output:
#     [
#       [2,4],
#       [3,4],
#       [2,3],
#       [1,2],
#       [1,3],
#       [1,4],
#     ]
#
#
# **Example 2:**
#
#
#
#     Input: n = 1, k = 1
#     Output: [[1]]
#
#
#
#
# **Constraints:**
#
#   * `1 <= n <= 20`
#   * `1 <= k <= n`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
#################################################################################################################
# AIMS 1-3:
# * Computes E[X] and E[XX] from model file used in aims 1-3
#################################################################################################################
using CpelNano
using LinearAlgebra

#################################################################################################################
# FUNCTIONS
#################################################################################################################
function process_model_file(chr::String,fasta::String,mod_file::String)::Nothing

    # IO
    outdir = dirname(mod_file)
    prefix = split(basename(mod_file),".")[1]
    ex_file = "$(outdir)/$(prefix)_ex.txt"
    exx_file = "$(outdir)/$(prefix)_exx.txt"

    # FASTA record
    fa_rec = CpelNano.get_chr_fa_rec(chr,fasta)

    # Go over all lines in file
    open(mod_file) do f
        for line in enumerate(eachline(f))

            # Get line data
            line_vec = split(line[2],"\t")

            # Get chromosome info
            reg_chr = line_vec[1]

            # Init struct
            rs = CpelNano.RegStruct()
            rs.chr = reg_chr

            # Get start or end
            rs.chrst = parse(Int64,line_vec[2])
            rs.chrend = parse(Int64,line_vec[3])
            CpelNano.print_log("st: $(rs.chrst) - end: $(rs.chrend)")

            # Get genome info
            CpelNano.get_grp_info!(rs,fa_rec,1)

            # Get αs & βs
            length(String(line_vec[5]))>0 || continue
            αs = parse.(Float64,String.(split(line_vec[5],",")))
            βs = parse.(Float64,String.(split(line_vec[6],",")))

            # Compute E[X] & E[XX]
            u1 = CpelNano.get_log_u(αs[1])
            uN = CpelNano.get_log_u(αs[end])
            Ws = [CpelNano.get_log_W(αs[n],αs[n+1],βs[n]) for n=1:length(βs)]
            Z = CpelNano.get_log_Z(u1,uN,Ws)
            ex = CpelNano.get_E_X_log(u1,uN,Ws,Z)
            exx = CpelNano.get_E_XX_log(u1,uN,Ws,Z)

            # Write E[X]
            open(ex_file,"a") do io
                for n=1:rs.N
                    write(io,"$(rs.chr)\t$(rs.cpg_pos[n])\t$(rs.cpg_pos[n])\t$(ex[n])\n")
                end
            end

            # Write E[XX]
            open(exx_file,"a") do io
                for n=1:(rs.N-1)
                    write(io,"$(rs.chr)\t$(rs.cpg_pos[n])\t$(rs.cpg_pos[n+1])\t$(exx[n])\n")
                end
            end

        end
    end

    # Return nothing
    return nothing

end
#################################################################################################################
# Calls
#################################################################################################################

# Model file aims 1-2
chr = "chr22"
fasta = "/Users/jordiabante/OneDrive - Johns Hopkins/CpelNano/Data/Simulations/reference/chr22.fa"
mod_file = "/Users/jordiabante/OneDrive - Johns Hopkins/CpelNano/Data/Simulations/ground_truth/GM12878_wgbs_cpelnano_theta_aug.txt"
process_model_file(chr,fasta,mod_file)

# Model file aim 3
chr = "chr22"
fasta = "/Users/jordiabante/OneDrive - Johns Hopkins/CpelNano/Data/Simulations/reference/chr22.fa"
mod_file = "/Users/jordiabante/OneDrive - Johns Hopkins/CpelNano/Data/Simulations/ground_truth/GM12878_wgbs_cpelnano_theta.txt"
process_model_file(chr,fasta,mod_file)


text:
using Drawing

ignore = false  # when re-generating

with(PNG("square-align.png", 100, 100), Axes(negative=true)) do
    draw() do
        square(0.6, align=9)
        square(0.8, align=5)
        square(1.0, align=1)
    end
end
ignore || compare("square-align.png")

with(PNG("rectangle-align.png", 100, 100), Axes(negative=true)) do
    draw() do
        rectangle(0.6, 0.3, align=9)
        rectangle(0.8, 0.4, align=5)
        rectangle(1.0, 0.5, align=1)
    end
end
ignore || compare("rectangle-align.png")