text:
module sphere_dipole_examples
using FinEtools
using LinearAlgebra
using Gaston

function sphere_dipole_transient()
    println("The interior sphere accelerates in the alternately in the positive
    and negative x-direction, generating positive pressure ahead of it, negative
    pressure behind. Time-dependent simulation.
    ")
    rho = 1.21*phun("kg/m^3");# mass density
    c  = 343.0*phun("m/s");# sound speed
    bulk =  c^2*rho;
    a_amplitude=1.0*phun("mm/s^2");# amplitude of the  acceleration of the sphere
    # omega = 2000*phun("rev/s");      # frequency of the incident wave
    R = 50.0*phun("mm"); # radius of the interior sphere
    Ro = 8*R # radius of the external sphere
    P_amplitude = R*rho*a_amplitude; # pressure amplitude
    nref=2;
    nlayers=40;
    tolerance = Ro/(nlayers)/100
    frequency=1200.; # Hz
    omega=2*pi*frequency;
    dt=1.0/frequency/20;
    tfinal=90*dt;
    nsteps=round(tfinal/dt)+1;

    println("""
            Spherical domain filled with acoustic medium, no scatterer.
    Hexahedral mesh.
            """)

    t0  =  time()

    # Hexahedral mesh
    fens,fes  =  H8sphere(R,nref);
    fens1,fes1  =  mirrormesh(fens, fes, [-1.0, 0.0, 0.0], [0.0, 0.0, 0.0],
    renumb =  r(c) = c[[1, 4, 3, 2, 5, 8, 7, 6]]);
    fens,newfes1,fes2 =  mergemeshes(fens1, fes1, fens, fes, tolerance)
    fes = cat(newfes1,fes2)

    # Derive the finite element sets for the boundary
    bfes  =  meshboundary(fes)
    # Outer spherical boundary
    function dout(xyz)
        return xyz/norm(xyz)
    end

    louter = selectelem(fens, bfes, facing = true, direction = dout)
    outer_fes = subset(bfes, louter);

    fens,fes  = H8extrudeQ4(fens, outer_fes, nlayers, (xyz, layer)->(R+layer/nlayers*(Ro-R))*xyz/norm(xyz));
    connected = findunconnnodes(fens, fes); fens, new_numbering = compactnodes(fens, connected);
    fess = renumberconn!(fes, new_numbering);

    geom  =  NodalField(fens.xyz)
    P  =  NodalField(zeros(FCplxFlt,size(fens.xyz,1),1))

    bfes  =  meshboundary(fes)

    # File  =   "Sphere.vtk"
    # vtkexportmesh(File, bfes.conn, geom.values, FinEtools.MeshExportModule.Q4)
    # @async run(`"paraview.exe" $File`)

    # File  =   "Sphere.vtk"
    # vtkexportmesh (File, fes.conn, fens.xyz, MeshExportModule.H8)
    # @async run(`"C:/Program Files (x86)/ParaView 4.2.0/bin/paraview.exe" $File`)

    linner = selectelem(fens, bfes, distance = R, from = [0.0 0.0 0.0], inflate = tolerance)
    louter = selectelem(fens, bfes, facing = true, direction = dout)

    println("Pre-processing time elapsed  =  ",time() - t0,"s")

    t1  =  time()

    numberdofs!(P)

    material = MatAcoustFluid(bulk,rho)
    femm  =  FEMMAcoust(IntegData(fes, GaussRule(3, 2)), material)

    S  =  acousticstiffness(femm, geom, P);
    C  =  acousticmass(femm, geom, P);

    abcfemm  =  FEMMAcoustSurf(IntegData(subset(bfes, louter), GaussRule(2, 2)), material)
    D  =  acousticABC(abcfemm, geom, P);

    # Inner sphere pressure loading
    function dipole(dpdn, xyz, J, label, t)
        n = cross(J[:,1],J[:,2]);
        n = vec(n/norm(n));
        dpdn[1] = -rho*a_amplitude*sin(omega*t)*n[1]
    end

    dipfemm  =  FEMMAcoustSurf(IntegData(subset(bfes, linner), GaussRule(2, 2)), material)
    # F  = distribloads(dipfemm, geom, P, fi, 2);
    #
    # K = lu((1.0+0.0im)*(-omega^2*S + omega*1.0im*D + C)) # We fake a complex matrix here
    # p = K\F  #
    #
    # scattersysvec!(P, p[:])

    # println(" Minimum/maximum pressure= $(minimum(real(p)))/$(maximum(real(p)))")

    println("Computing time elapsed  =  ",time() - t1,"s")
    println("Total time elapsed  =  ",time() - t0,"s")

    # Solve
    P0 = deepcopy(P)
    P0.values[:] = 0.0; # initially all pressure is zero
    vP0 = gathersysvec(P0);
    vP1 = zeros(eltype(vP0), size(vP0));
    vQ0 = zeros(eltype(vP0), size(vP0));
    vQ1 = zeros(eltype(vP0), size(vP0));
    t = 0.0;
    P1 = deepcopy(P0);

    fi  =  ForceIntensity(FCplxFlt, 1, (dpdn, xyz, J, label)->dipole(dpdn, xyz, J, label, t));
    La0 = distribloads(dipfemm, geom, P1, fi, 2);

    A = (2.0/dt)*S + D + (dt/2.)*C;

    step =0;
    while t <= tfinal
      step = step  + 1;
      println("Time $t ($(step)/$(round(tfinal/dt)+1))")
      # if graphics
      #   Camera =1.0e+003 * [-1.4303   -3.0332   -2.0483    0.1553    0.2912    0.2533         0         0    0.0010    0.0066];
      #   wv = get(P1,'values');
      #   dcm=data_colormap(struct ('range',P_amplitude*[-1,1], 'colormap',bwr));
      #   colors=map_data(dcm, wv);
      #   colorfield = field(struct ('name',['cf'], 'dim', 3, 'data',colors));
      #   gv=reset (gv,[]);
      #   camset(gv,Camera);
      #   draw(bfeb,gv, struct ('x', geom,'u',0*geom,'colorfield',colorfield));
      #   draw_axes (gv,struct('length',R));
      #   pause(0.55);
      # end
      t = t+dt;
      fi  =  ForceIntensity(FCplxFlt, 1, (dpdn, xyz, J, label)->dipole(dpdn, xyz, J, label, t));
      La1 = distribloads(dipfemm, geom, P1, fi, 2);
      vQ1 = A\((2/dt)*(S*vQ0)-D*vQ0-C*(2*vP0+(dt/2)*vQ0)+La0+La1);
      vP1 = vP0 + (dt/2)*(vQ0+vQ1);
      vP0 = vP1;
      vQ0 = vQ1;
      P1 = scattersysvec!(P1, vec(vP1));
      P0 = P1;
      La0 = La1;
    end

    println("$(vP1[1])")

    File  =   "sphere_dipole_1.vtk"
    vtkexportmesh(File, fes.conn, geom.values, FinEtools.MeshExportModule.H8;
     scalars = [( "realP", real(P1.values))])
    @async run(`"paraview.exe" $File`)

    true
end # sphere_dipole_transient

function allrun()
    println("#####################################################")
    println("# sphere_dipole_transient ")
    sphere_dipole_transient()
    return true
end # function allrun

end # module baffled_piston_examples


text:
struct Sherlog32{T<:AbstractFloat,i} <: AbstractSherlog     # T is the bitpattern type for logging, i is the logbook #
    val::Float32                                            # the value is always Float32
end

# conversions back from Sherlog32
Base.Float64(x::Sherlog32) = Float64(x.val)
Base.Float32(x::Sherlog32) = x.val
Base.Float16(x::Sherlog32) = Float16(x.val)
Base.Int64(x::Sherlog32) = Int64(x.val)
Base.Int32(x::Sherlog32) = Int32(x.val)
Base.Int16(x::Sherlog32) = Int16(x.val)

# otherwise conversion to itself is ambigious
Sherlog32{T,i}(x::Sherlog32) where {T<:AbstractFloat,i} = Sherlog32{T,i}(x.val)
Sherlog32(x::Sherlog32{T,i}) where {T<:AbstractFloat,i} = Sherlog32{T,i}(x.val)

# generator functions from Float, Int and Bool
Sherlog32(x::AbstractFloat) = Sherlog32{Float16,1}(x)   # use Float16 and logbook #1 as standard
Sherlog32{T}(x::AbstractFloat) where T = Sherlog32{T,1}(x)      # use logbook #1 if no logbook provided
Sherlog32(x::Integer) = Sherlog32{Float16,1}(x)           # same for integers
Sherlog32{T}(x::Integer) where T = Sherlog32{T,1}(x)
Sherlog32{T}(x::Bool) where T = if x Sherlog32{T,1}(1) else Sherlog32{T,1}(0) end

Base.promote_rule(::Type{Int64},::Type{Sherlog32{T,i}}) where {T,i} = Sherlog32{T,i}
Base.promote_rule(::Type{Int32},::Type{Sherlog32{T,i}}) where {T,i} = Sherlog32{T,i}
Base.promote_rule(::Type{Int16},::Type{Sherlog32{T,i}}) where {T,i} = Sherlog32{T,i}

Base.promote_rule(::Type{Float64},::Type{Sherlog32{T,i}}) where {T,i} = Sherlog32{T,i}
Base.promote_rule(::Type{Float32},::Type{Sherlog32{T,i}}) where {T,i} = Sherlog32{T,i}
Base.promote_rule(::Type{Float16},::Type{Sherlog32{T,i}}) where {T,i} = Sherlog32{T,i}

Base.bitstring(x::Sherlog32) = bitstring(x.val)
Base.show(io::IO,x::Sherlog32) = print(io,"Sherlog32(",string(x.val),")")

Base.eps(::Type{Sherlog32{T,i}}) where {T,i} = eps(Float32)
Base.eps(x::Sherlog32) = eps(x.val)

Base.typemin(::Sherlog32{T,i}) where {T,i} = Sherlog32{T,i}(typemin(Float32))
Base.typemax(::Sherlog32{T,i}) where {T,i} = Sherlog32{T,i}(typemax(Float32))
Base.floatmin(::Sherlog32{T,i}) where {T,i} = Sherlog32{T,i}(floatmin(Float32))
Base.floatmax(::Sherlog32{T,i}) where {T,i} = Sherlog32{T,i}(floatmax(Float32))
Base.precision(::Sherlog32{T,i}) where {T,i} = Sherlog32{T,i}(precision(Float32))

function +(x::Sherlog32{T,i},y::Sherlog32{T,i}) where {T,i}
    r = x.val + y.val
    log_it(T,r,i)
    return Sherlog32{T,i}(r)
end

function -(x::Sherlog32{T,i},y::Sherlog32{T,i}) where {T,i}
    r = x.val - y.val
    log_it(T,r,i)
    return Sherlog32{T,i}(r)
end

function *(x::Sherlog32{T,i},y::Sherlog32{T,i}) where {T,i}
    r = x.val * y.val
    log_it(T,r,i)
    return Sherlog32{T,i}(r)
end

function /(x::Sherlog32{T,i},y::Sherlog32{T,i}) where {T,i}
    r = x.val / y.val
    log_it(T,r,i)
    return Sherlog32{T,i}(r)
end

function ^(x::Sherlog32{T,i},y::Sherlog32{T,i}) where {T,i}
    r = x.val ^ y.val
    log_it(T,r,i)
    return Sherlog32{T,i}(r)
end

for O in ( :(-), :(+),
           :sign,
           :prevfloat, :nextfloat,
           :round, :trunc, :ceil, :floor,
           :inv, :abs, :sqrt, :cbrt,
           :exp, :expm1, :exp2, :exp10,
           :log, :log1p, :log2, :log10,
           :rad2deg, :deg2rad, :mod2pi, :rem2pi,
           :sin, :cos, :tan, :csc, :sec, :cot,
           :asin, :acos, :atan, :acsc, :asec, :acot,
           :sinh, :cosh, :tanh, :csch, :sech, :coth,
           :asinh, :acosh, :atanh, :acsch, :asech, :acoth,
           :sinc, :sinpi, :cospi,
           :sind, :cosd, :tand, :cscd, :secd, :cotd,
           :asind, :acosd, :atand, :acscd, :asecd, :acotd
          )
    @eval begin
        function Base.$O(x::Sherlog32{T,i}) where {T,i}
            r = $O(x.val)
            log_it(T,r,i)
            return Sherlog32{T,i}(r)
        end
    end
end

for O in ( :(<), :(<=))
    @eval begin
        Base.$O(x::Sherlog32, y::Sherlog32) = $O(x.val, y.val)
    end
end

ext: jl


text:
export particle_type

const EXT_HEADER = ".IAEAheader"
const EXT_PHSP   = ".IAEAphsp"
struct IAEAPath
    header::String
    phsp::String
end
IAEAPath(path::IAEAPath) = path

function IAEAPath(path::AbstractString)
    stem, ext = splitext(path)
    if !(ext in (EXT_HEADER, EXT_PHSP))
        stem = path
    end
    header_path = stem * EXT_HEADER
    phsp_path = stem * EXT_PHSP
    IAEAPath(header_path, phsp_path)
end

function Base.rm(path::IAEAPath;kw...)
    rm(path.header;kw...)
    rm(path.phsp;kw...)
end

struct IAEAParticle{Nf, Ni} <: AbstractIAEAParticle
    typ::ParticleType
    E::Float32
    weight::Float32
    x::Float32
    y::Float32
    z::Float32
    u::Float32
    v::Float32
    w::Float32
    new_history::Bool
    extra_floats::NTuple{Nf, Float32}
    extra_ints::NTuple{Ni, Int32}
end

function equal_properties(p1::AbstractIAEAParticle, p2::AbstractIAEAParticle)
    p1.typ    == p2.typ || return false
    p1.E      == p2.E || return false
    p1.weight == p2.weight || return false
    p1.x      == p2.x || return false
    p1.y      == p2.y || return false
    p1.z      == p2.z || return false
    p1.u      == p2.u || return false
    p1.v      == p2.v || return false
    p1.w      == p2.w || return false
    p1.new_history  == p2.new_history || return false
    p1.extra_floats == p2.extra_floats || return false
    p1.extra_ints   == p2.extra_ints || return false
    return true
end

function Base.:(==)(p1::AbstractIAEAParticle, p2::AbstractIAEAParticle)
    equal_properties(p1, p2)
end

function IAEAParticle(typ, E, weight,
                          x,y,z,
                          u,v,w,
                          new_history,
                          extra_floats::NTuple{Nf,Float32},
                          extra_ints::NTuple{Ni,Int32}) where {Nf, Ni}

    IAEAParticle{Nf, Ni}(typ, E, weight,
                     x,y,z, u,v,w,
                     new_history,
                     extra_floats, extra_ints)
end

function IAEAParticle(;typ,E,weight=1,x,y,z,
                  u,v,w,
                  new_history=true,
                  extra_floats=(),
                  extra_ints=()
                 )
    IAEAParticle(typ, E, weight,
                     x,y,z, u,v,w,
                     new_history,
                     extra_floats, extra_ints)
end

particle_type(p::IAEAParticle) = p.typ
function Base.show(io::IO, p::IAEAParticle)
    kwshow(io, p)
end

function Base.isapprox(p1::AbstractIAEAParticle, p2::AbstractIAEAParticle;kw...)
    p1.typ == p2.typ  &&
    p1.new_history   == p2.new_history    &&
    p1.extra_floats  == p2.extra_floats   &&
    p1.extra_ints    == p2.extra_ints     &&
    isapprox(p1.E,      p2.E;      kw...) &&
    isapprox(p1.weight, p2.weight; kw...) &&
    isapprox(position(p1), position(p2); kw...) &&
    isapprox(direction(p1), direction(p2); kw...)
end

struct IAEAHeader{Nf, Ni, NT<:NamedTuple}
    record_contents::NT
    attributes::OrderedDict{Symbol, String}
    function IAEAHeader{Nf, Ni, NT}(
            rc::NT=NamedTuple(),
            raw=OrderedDict{Symbol,String}()) where {Nf, Ni, NT <: NamedTuple}
        @argcheck Ni isa Int
        @argcheck Nf isa Int
        for (keyword,val) in pairs(rc)
            @argcheck keyword in [:x, :y, :z, :u, :v, :w, :weight]
            @argcheck val isa Float32
        end
        new{Nf, Ni,NT}(rc, raw)
    end
end

function IAEAHeader{Nf, Ni}(
        rc::NamedTuple=NamedTuple(),
        raw=OrderedDict{Symbol,String}()) where {Nf, Ni}

    for (keyword,val) in pairs(rc)
        @argcheck keyword in [:x, :y, :z, :u, :v, :w, :weight]
    end
    rc32 = map(Float32, rc)
    NT = typeof(rc32)
    IAEAHeader{Nf, Ni,NT}(rc32, raw)
end

function get_originalcount(o::IAEAHeader)::Float64
    parse(Float64, o.attributes[:ORIG_HISTORIES])
end

struct IAEAPhspIterator{H,I<:IO} <: AbstractPhspIterator
    io::I
    header::H
    length::Int64
end

function IAEAPhspIterator(io::IO,h::IAEAHeader)
    H = typeof(h)
    I = typeof(io)
    bl = bytelength(io)
    length = Int64(bl / ptype_disksize(h))
    IAEAPhspIterator{H,I}(io, h, length)
end

function Base.eltype(::Type{<:IAEAPhspIterator{H}}) where {H}
    ptype(H)
end


text:
module DOM

using Hyperscript

const global_unique_dom_id_counter = Ref(0)

"""
    get_unique_dom_id()

We could use a unique ID like uuid4, but since every dom element gets
such an id, I prefer to keep the id as short as possible, so we just use a counter.
"""
function get_unique_dom_id()
    global_unique_dom_id_counter[] += 1
    return string(global_unique_dom_id_counter[])
end

"""
Dome node with unique ID, to make it easier to interpolate it.
"""
function um(tag, args...; kw...)
    m(tag, args..., dataJscallId = get_unique_dom_id(); kw...)
end

function m_unesc(tag, args...; kw...)
    m(Hyperscript.NOESCAPE_HTMLSVG_CONTEXT, tag, args..., dataJscallId = get_unique_dom_id(); kw...)
end

p(args...; kw...) = um("p", args...; kw...)
div(args...; kw...) = um("div", args...; kw...)
div_unesc(args...; kw...) = m_unesc("div", args...; kw...)
input(args...; kw...) = um("input", args...; kw...)
font(args...; kw...) = um("font", args...; kw...)

style(args...; kw...) = m(Hyperscript.NOESCAPE_HTMLSVG_CONTEXT, "style", args...; kw...)
script(args...; kw...) = m(Hyperscript.NOESCAPE_HTMLSVG_CONTEXT, "script", args...; kw...)

end

using .DOM

# default turn attributes into strings
attribute_render(session, parent, attribute, x) = string(x)
function attribute_render(session, parent, attribute, obs::Observable)
    onjs(session, obs, js"""
    function (value){
        var node = $(parent);
        if(node){
            if(node[$attribute] != value){
                node[$attribute] = value;
            }
            return true;
        }else{
            return false; //deregister
        }
    }
    """)
    return attribute_render(session, parent, attribute, obs[])
end

function attribute_render(session, parent, attribute, jss::JSCode)
    return serialize_string(jss)
end

render_node(session::Session, x) = x

function render_node(session::Session, node::Node{Hyperscript.CSS})
# do nothing
    return node
end
function render_node(session::Session, node::Node)
    # give each node a unique id inside the dom
    new_attributes = Dict{String, Any}()
    newchildren = map(children(node)) do elem
        childnode = jsrender(session, elem)
        # if a transform elem::Any -> ::Node happens, we need to
        # render the resulting node again, since the attr/children won't be
        # lowered yet!
        if !(elem isa Node)
            childnode = render_node(session, childnode)
        end
        return childnode
    end
    for (k, v) in Hyperscript.attrs(node)
        new_attributes[k] = attribute_render(session, node, k, v)
    end
    return Node(
        Hyperscript.context(node),
        Hyperscript.tag(node),
        newchildren,
        new_attributes
    )
end

# jsrender(session, x) will be called anywhere...
# if there is nothing sessions specific in the dom, fallback to jsrender without session
function jsrender(session::Session, node::Node)
    render_node(session, node)
end

function uuid(node::Node)
    get(Hyperscript.attrs(node), "data-jscall-id") do
        error("Node $(node) doesn't have a unique id. Make sure to use DOM.$(Hyperscript.tag(node))")
    end
end

# Handle interpolating into Javascript
function serialize_string(io::IO, node::Node)
    # This relies on jsrender to give each node a unique id under the
    # attribute data-jscall-id. This is a bit brittle
    # improving this would be nice
    print(io, "(document.querySelector('[data-jscall-id=$(repr(uuid(node)))]'))")
end

"""
    jsrender([::Session], x::Any)
Internal render method to create a valid dom. Registers used observables with a session
And makes sure the dom only contains valid elements. Overload jsrender(::YourType)
To enable putting YourType into a dom element/div.
You can also overload it to take a session as first argument, to register
messages with the current web session (e.g. via onjs).
"""
jsrender(::Session, x::Any) = jsrender(x)
jsrender(::Session, x::Union{Symbol, String}) = DOM.p(string(x))

jsrender(::Session, x::String) = x
jsrender(::Session, x::Hyperscript.Styled) = x


text:
include("run_scf_and_compare.jl")
include("testcases.jl")

function run_oxygen_pbe(T; kwargs...)
    # These values were computed using ABINIT at Ecut=13. They are not yet converged and
    # thus require the same discretisation parameters to be obtained.
    ref_evals = [
        [ -1.214034187039492, -0.617266880915023, -0.370358278719988, -0.370358278712184,
          -0.352666293942291,  0.008939230950524,  0.008939230958503,  0.098061783836035],
        [ -1.177293790864051, -0.551559952680091, -0.318494822087423, -0.308762934378086,
          -0.308762934370473,  0.090371516203330,  0.090371516211105,  0.115096531941182],
    ]
    ref_etot = -29.505154342215235
    ref_magn = 1.99231275

    # Produce reference data and guess for this configuration
    Ecut = 13
    O = ElementPsp(o2molecule.atnum, psp=load_psp("hgh/pbe/O-q6.hgh"))
    magnetic_moments = [O => [1., 1.]]
    model = model_PBE(Array{T}(o2molecule.lattice), [O => o2molecule.positions],
                      temperature=0.02, smearing=smearing=Smearing.Gaussian(),
                      magnetic_moments=magnetic_moments)
    basis = PlaneWaveBasis(model, Ecut; fft_size=[24, 24, 30], kgrid=[1, 1, 1])

    ρspin = guess_spin_density(basis, magnetic_moments)
    scfres = run_scf_and_compare(T, basis, ref_evals, ref_etot; ρspin=ρspin,
                                 test_etot=false, kwargs...)
    @test scfres.energies.total ≈ ref_etot atol=1e-4  # A little large a difference ...

    dVol = scfres.basis.model.unit_cell_volume / prod(scfres.basis.fft_size)
    magnetisation = sum(scfres.ρspin.real) * dVol
    @test magnetisation ≈ ref_magn atol=5e-6
end

@testset "Oxygen PBE (Float64)" begin
    run_oxygen_pbe(Float64, test_tol=5e-6, scf_tol=1e-8, n_ignored=1)
end


text:
using Convex
using Convex.ProblemDepot: run_tests
using Test
using SCS, ECOS, GLPKMathProgInterface

# Seed random number stream to improve test reliability
using Random
Random.seed!(2)

@testset "ProblemDepot" begin
    @testset "Problems can run without `solve!`ing if `test==false`" begin
        Convex.ProblemDepot.foreach_problem() do name, func
            @testset "$name" begin
                # We want to check to make sure this does not throw
                func(Convex.conic_problem, Val(false), 0.0, 0.0, Float64)
                @test true
            end
        end
    end
end

@testset "Convex" begin
    include("test_utilities.jl")

    @testset "SCS" begin
        run_tests(; exclude=[r"mip"]) do p
            solve!(p, SCSSolver(verbose=0, eps=1e-6))
        end
    end

    @testset "ECOS" begin
        run_tests(; exclude=[r"mip", r"sdp"]) do p
            solve!(p, ECOSSolver(verbose=0))
        end
    end

    @testset "GLPK MIP" begin
        run_tests(; exclude=[r"socp", r"sdp", r"exp", r"dual"]) do p
            solve!(p, GLPKSolverMIP())
        end
    end
end


text:
using ADCME
using PyCall
using LinearAlgebra
using PyPlot
using Random
Random.seed!(233)

function halo_exchange_two_d(u,fill_value,m,n)
    halo_exchange_two_d_ = load_op_and_grad("./build/libHaloExchangeTwoD","halo_exchange_two_d")
    u,fill_value,m,n = convert_to_tensor(Any[u,fill_value,m,n], [Float64,Float64,Int64,Int64])
    halo_exchange_two_d_(u,fill_value,m,n)
end

mpi_init()
U = reshape(1:24, 4, 6)'|>Array

m = 3
n = 2
fill_value = 1.0

M = mpi_rank()÷n+1
N = mpi_rank()%n+1
ulocal = U[(M-1)*2 + 1: M * 2, (N-1)*2+1:N*2]

sess = Session(); init(sess)

# TODO: change your test parameter to `m`
#       in the case of `multiple=true`, you also need to specify which component you are testings
# gradient check -- v
function scalar_function(x)
    x = mpi_bcast(x)
    u = x .* ulocal
    return mpi_sum(sum(halo_exchange_two_d(u,fill_value,m,n)^2))
end

# TODO: change `m_` and `v_` to appropriate values
m_ = constant(rand(2,2))
v_ = rand(2,2)
y_ = scalar_function(m_)
dy_ = gradients(y_, m_)
ms_ = Array{Any}(undef, 5)
ys_ = Array{Any}(undef, 5)
s_ = Array{Any}(undef, 5)
w_ = Array{Any}(undef, 5)
gs_ =  @. 1 / 10^(1:5)

for i = 1:5
    g_ = gs_[i]
    ms_[i] = m_ + g_*v_
    ys_[i] = scalar_function(ms_[i])
    s_[i] = ys_[i] - y_
    w_[i] = s_[i] - g_*sum(v_.*dy_)
end

sess = Session(); init(sess)
sval_ = run(sess, s_)
wval_ = run(sess, w_)

if mpi_rank()==0
    close("all")
    loglog(gs_, abs.(sval_), "*-", label="finite difference")
    loglog(gs_, abs.(wval_), "+-", label="automatic differentiation")
    loglog(gs_, gs_.^2 * 0.5*abs(wval_[1])/gs_[1]^2, "--",label="\$\\mathcal{O}(\\gamma^2)\$")
    loglog(gs_, gs_ * 0.5*abs(sval_[1])/gs_[1], "--",label="\$\\mathcal{O}(\\gamma)\$")

    plt.gca().invert_xaxis()
    legend()
    xlabel("\$\\gamma\$")
    ylabel("Error")
    savefig("gradtest.png")
end

mpi_finalize()


text:
"""
    is_valid_executable_document(str::String; throw_on_error=false)

Returns a `Bool` indicating whether the document described by `str` is valid.

The document is parsed and some validation performed.
For further information see [`validate_executable_document`](@ref) and package documentation.

To throw an exception if a validation error is found, set `throw_on_error` to `true`.
(Note, parsing errors will always throw an exception).

To retrieve a list of validation errors in the document, use [`validate_executable_document`](@ref) instead.
"""
function is_valid_executable_document(str::String; throw_on_error=false)
    errors = validate_executable_document(str)
    if throw_on_error && !isempty(errors)
        throw(ValidationException(errors))
    end
    return isempty(errors)
end

"""
    validate_executable_document(str::String)

Return a list of validation errors in the GraphQL executable document described by `str`.

Firstly the document will be parsed with any parsing errors being immediately thrown.

Secondly, the parsed document will be validated against **some** of the specification with all validation errors being returned.
See the package documentation for a full description of what validation is performed.
"""
validate_executable_document(str::String) = validate(parse(str))

function validate(doc::Document)
    errors = ValidationError[]
    validate_operations!(errors::Vector{ValidationError}, doc::Document)
    validate_fragments!(errors::Vector{ValidationError}, doc::Document)
    return errors
end

validate_variable_definitions!(errors, ::Nothing, used_variables, op) = errors
function validate_variable_definitions!(errors, var_defs::Vector{VariableDefinition}, used_variables, op::Operation)
    variable_names = (var_def.name for var_def in var_defs)
    for var_def in var_defs
        if count(==(var_def.name), variable_names) > 1
            push!(
                errors,
                RepeatedVariableDefinition(
                    "There can only be one variable named \"$(var_def.name)\".",
                    [var_def.loc]
                )
            )
        end
        if var_def.name ∉ used_variables
            push!(
                errors,
                UnusedVariable(
                    "Variable \"\$(var_def.name\" is never used in operation \"$(op.name)\".",
                    [var_def.loc, op.loc]
                )
            )
        end
    end
    return errors
end

function validate_operation!(errors, op::Operation, defined_fragments)
    used_variables = String[]
    validate_directive_uniqueness!(errors, op.directives)
    validate_arguments!(errors, used_variables, op.directives, op)
    validate_selection_set!(errors, used_variables, op.selection_set, op, defined_fragments)
    validate_variable_definitions!(errors, op.variable_definitions, used_variables, op)
    return errors
end

validate_selection_set!(errors, used_variables, ::Nothing, op::Operation, defined_fragments) = errors
function validate_selection_set!(errors, used_variables, ss::SelectionSet, op::Operation, defined_fragments)
    for selection in ss.selections
        selection isa FragmentSpread && validate_fragment_spread!(errors, used_variables, selection, op, defined_fragments)
        selection isa InlineFragment && validate_inline_fragment!(errors, used_variables, selection, op, defined_fragments)
        selection isa Field && validate_field!(errors, used_variables, selection, op, defined_fragments)
    end
    # validate merging of fields
    return errors
end

function validate_fragment_spread!(errors, used_variables, spread::FragmentSpread, op::Operation, defined_fragments)
    if isnothing(defined_fragments) || isnothing(findfirst(def -> (def.name==spread.name), defined_fragments))
        push!(
            errors,
            UnknownFragment(
                "Unknown fragment \"$(spread.name)\"",
                [spread.loc]
            )
        )
        return errors
    end
    fragment_definition = defined_fragments[findfirst(def -> (def.name==spread.name), defined_fragments)]
    validate_arguments!(errors, used_variables, fragment_definition, op)
    validate_selection_set!(errors, used_variables, fragment_definition.selection_set, op, defined_fragments)
    # Don't need to check directive Uniqueness here as that is done at fragment definition
    return errors
end

function validate_inline_fragment!(errors, used_variables, inline_frag::InlineFragment, op::Operation, defined_fragments)
    validate_arguments!(errors, used_variables, inline_frag, op)
    validate_directive_uniqueness!(errors, inline_frag.directives)
    validate_selection_set!(errors, used_variables, inline_frag.selection_set, op, defined_fragments)
    return errors
end

validate_field!(errors, used_variables, ::Nothing, op::Operation, defined_fragments) = errors
function validate_field!(errors, used_variables, field::Field, op::Operation, defined_fragments)
    # Validate arguments
    validate_arguments!(errors, used_variables, field.arguments, op)
    # Validate directives
    validate_arguments!(errors, used_variables, field.directives, op)
    validate_directive_uniqueness!(errors, field.directives)

    # validate selection set
    validate_selection_set!(errors, used_variables, field.selection_set, op, defined_fragments)
    return errors
end

validate_directive_uniqueness!(errors, ::Nothing) = errors
function validate_directive_uniqueness!(errors, directives::Vector{Directive})
    directive_names = (directive.name for directive in directives)
    for directive in unique(x -> x.name, directives)
        if count(==(directive.name), directive_names) > 1
            push!(
                errors,
                RepeatedDirectiveName(
                    "The directive \"$(directive.name)\" can only be used once at this location.",
                    [dir.loc for dir in directives if dir.name == directive.name]
                )
            )
        end
    end
    return errors
end

validate_arguments!(errors, used_variables, ::Nothing, op::Operation) = errors
function validate_arguments!(errors, used_variables, fragment::Union{FragmentDefinition, InlineFragment}, op::Operation)
    validate_arguments!(errors, used_variables, fragment.directives, op)
    return errors
end
function validate_arguments!(errors, used_variables, directives::Vector{Directive}, op::Operation)
    for directive in directives
        validate_arguments!(errors, used_variables, directive.arguments, op)
    end
    return errors
end
function validate_arguments!(errors, used_variables, arguments::Vector{Argument}, op::Operation)
    argument_names = (argument.name for argument in arguments)
    for argument in arguments
        if count(==(argument.name), argument_names) > 1
            push!(
                errors,
                RepeatedArgumentName(
                    "There can only be one argument named \"$(op.name)\".",
                    [op.loc]
                )
            )
        end
        if argument.value isa Variable
            var = argument.value
            if isnothing(op.variable_definitions) || !any(==(var.name), (def.name for def in op.variable_definitions))
                push!(
                    errors,
                    UnknownVariable(
                        "Variable \"\$$(var.name)\" is not defined by operation \"$(op.name)\".",
                        [var.loc, op.loc]
                    )
                )
            end
            push!(used_variables, var.name)
        elseif argument.value isa InputObject
            input_object = argument.value
            field_names = (field.name for field in input_object.object_fields)
            for field in unique(x -> x.name, input_object.object_fields)
                if count(==(field.name), field_names) > 1
                    push!(
                        errors,
                        RepeatedInputObjectField(
                            "There can only be one input object field named \"$(field.name)\".",
                            [f.loc for f in input_object.object_fields if f.name == field.name]
                        )
                    )
                end
            end
        end
    end
    return errors
end

"""
    get_defined_operations(doc::Document)

Returns a generator over all operations that are defined in the document.
"""
get_defined_operations(doc::Document) = (def for def in doc.definitions if isa(def, Operation))

function validate_operations!(errors::Vector{ValidationError}, doc::Document)
    defined_operations = get_defined_operations(doc)
    defined_names = [def.name for def in defined_operations]
    n_defined_names = length(defined_names)

    defined_fragments = get_defined_fragments(doc)
    for op in defined_operations
        # Check anonymouse operation is alone - 5.2.2.1
        if isnothing(op.name) && n_defined_names > 1
            push!(
                errors,
                AnonymousOperationNotAlone(
                    # TODO: message
                    "This anonymous operation must be the only defined operation.",
                    [op.loc]
                )
            )
        end
        # Check operation names are unique - 5.2.1.1
        if count(==(op.name), defined_names) > 1
            push!(
                errors,
                RepeatedOperationDefinition(
                    "There can only be one Operation named \"$(op.name)\".",
                    [op.loc]
                )
            )
        end
        validate_operation!(errors, op, defined_fragments)
    end

    return errors
end

"""
    get_defined_fragments(doc::Document)

Returns a vector of all fragments that are defined in the document.
"""
get_defined_fragments(doc::Document) = [def for def in doc.definitions if isa(def, FragmentDefinition)]

function validate_fragments!(errors::Vector{ValidationError}, doc::Document)
    defined_fragments = get_defined_fragments(doc)
    defined_names = [def.name for def in defined_fragments]

    # Check for multiple definitions of same fragment - 5.5.1.1
    if length(defined_names) > length(unique(defined_names))
        for fragment in defined_fragments
            if count(==(fragment.name), defined_names) > 1
                push!(
                    errors,
                    RepeatedFragmentDefinition(
                        "There can only be one fragment named \"$(fragment.name))\".",
                        [fragment.loc]
                    )
                )
            end
        end
    end

    # Check all defined fragments are used - 5.5.2.1
    used_fragments = find_fragments(doc)
    used_names = keys(used_fragments)
    for fragment in defined_fragments
        if fragment.name ∉ used_names
            push!(
                errors,
                UnusedFragment(
                    "Fragment \"$(fragment.name)\" is never used.",
                    [fragment.loc]
                )
            )
        end
        validate_directive_uniqueness!(errors, fragment.directives)
    end

    return errors
end

find_fragments(doc) = find_fragments!(Dict{String, Vector{Loc}}(), doc)
find_fragments!(used_fragments, doc::Document) = find_fragments!(used_fragments, doc.definitions)
find_fragments!(used_fragments, frag_def::FragmentDefinition) = find_fragments!(used_fragments, frag_def.selection_set)
find_fragments!(used_fragments, op::Operation) = find_fragments!(used_fragments, op.selection_set)
find_fragments!(used_fragments, ss::SelectionSet) = find_fragments!(used_fragments, ss.selections)
find_fragments!(used_fragments, field::Field) = find_fragments!(used_fragments, field.selection_set)
find_fragments!(used_fragments, inline_frag::InlineFragment) = find_fragments!(used_fragments, inline_frag.selection_set)
find_fragments!(used_fragments, objs::Vector) = (foreach((obj) -> find_fragments!(used_fragments, obj), objs); used_fragments)
find_fragments!(used_fragments, ::Nothing) = used_fragments
function find_fragments!(used_fragments, frag_spread::FragmentSpread)
    if haskey(used_fragments, frag_spread.name)
        push!(used_fragments[frag_spread.name], frag_spread.loc)
    else
        used_fragments[frag_spread.name] = [frag_spread.loc]
    end
    return used_fragments
end


text:
# Taken from https://github.com/simonster/Reexport.jl/blob/master/src/Reexport.jl by https://github.com/simonster
# As no release has been made yet: https://github.com/simonster/Reexport.jl/pull/23
macro reexport(ex)
    isa(ex, Expr) && (ex.head == :module ||
                      ex.head == :using ||
                      (ex.head == :toplevel &&
                       all(e->isa(e, Expr) && e.head == :using, ex.args))) ||
        error("@reexport: syntax error")

    if ex.head == :module
        modules = Any[ex.args[2]]
        ex = Expr(:toplevel, ex, :(using .$(ex.args[2])))
    elseif ex.head == :using && all(e->isa(e, Symbol), ex.args)
        modules = Any[ex.args[end]]
    elseif ex.head == :using && ex.args[1].head == :(:)
        symbols = [e.args[end] for e in ex.args[1].args[2:end]]
        return esc(Expr(:toplevel, ex, :(eval(Expr(:export, $symbols...)))))
    else
        modules = Any[e.args[end] for e in ex.args]
    end

    esc(Expr(:toplevel, ex,
             [:(eval(Expr(:export, names($mod)...))) for mod in modules]...))
end

@reexport using ArchGDAL: intersects, equals, disjoint, touches, crosses, within, contains, overlaps
@reexport using ArchGDAL: boundary, convexhull, buffer
@reexport using ArchGDAL: intersection, union, difference, symdifference, distance
@reexport using ArchGDAL: geomlength, geomarea, centroid
@reexport using ArchGDAL: isempty, isvalid, issimple, isring, geomarea, centroid
@reexport using ArchGDAL: createpoint, createlinestring, createlinearring, createpolygon, createmultilinestring, createmultipolygon
@reexport using ArchGDAL: reproject

AG.intersects(a::Vector{AG.IGeometry}, b::Vector{AG.IGeometry}) = AG.intersects.(a, b)
AG.equals(a::Vector{AG.IGeometry}, b::Vector{AG.IGeometry}) = AG.equals.(a, b)
AG.disjoint(a::Vector{AG.IGeometry}, b::Vector{AG.IGeometry}) = AG.disjoint.(a, b)
AG.touches(a::Vector{AG.IGeometry}, b::Vector{AG.IGeometry}) = AG.touches.(a, b)
AG.crosses(a::Vector{AG.IGeometry}, b::Vector{AG.IGeometry}) = AG.crosses.(a, b)
AG.within(a::Vector{AG.IGeometry}, b::Vector{AG.IGeometry}) = AG.within.(a, b)
AG.contains(a::Vector{AG.IGeometry}, b::Vector{AG.IGeometry}) = AG.contains.(a, b)
AG.overlaps(a::Vector{AG.IGeometry}, b::Vector{AG.IGeometry}) = AG.overlaps.(a, b)
AG.intersection(a::Vector{AG.IGeometry}, b::Vector{AG.IGeometry}) = AG.intersection.(a, b)
AG.union(a::Vector{AG.IGeometry}, b::Vector{AG.IGeometry}) = AG.union.(a, b)
AG.difference(a::Vector{AG.IGeometry}, b::Vector{AG.IGeometry}) = AG.difference.(a, b)
AG.symdifference(a::Vector{AG.IGeometry}, b::Vector{AG.IGeometry}) = AG.symdifference.(a, b)
AG.distance(a::Vector{AG.IGeometry}, b::Vector{AG.IGeometry}) = AG.distance.(a, b)

AG.boundary(v::Vector{AG.IGeometry}) = AG.boundary.(v)
AG.convexhull(v::Vector{AG.IGeometry}) = AG.convexhull.(v)
AG.buffer(v::Vector{AG.IGeometry}, d) = AG.buffer.(v, d)
AG.transform!(v::Vector{AG.IGeometry}, d) = AG.buffer.(v, d)
AG.geomlength(v::Vector{AG.IGeometry}) = AG.geomlength.(v)
AG.geomarea(v::Vector{AG.IGeometry}) = AG.geomarea.(v)
AG.centroid(v::Vector{AG.IGeometry}) = AG.centroid.(v)
AG.isempty(v::Vector{AG.IGeometry}) = AG.isempty.(v)
AG.isvalid(v::Vector{AG.IGeometry}) = AG.isvalid.(v)
AG.issimple(v::Vector{AG.IGeometry}) = AG.issimple.(v)
AG.isring(v::Vector{AG.IGeometry}) = AG.isring.(v)


text:
using Crystallography:
    CrystalSystem,
    Cubic,
    Hexagonal,
    Trigonal,
    Tetragonal,
    Orthorhombic,
    Monoclinic,
    Triclinic

export solve_elastic_constants, solve_stiffnesses, solve_compliances

function recombine_strains(::Cubic, strain::EngineeringStrain)
    ϵ₁, ϵ₂, ϵ₃, ϵ₄, ϵ₅, ϵ₆ = strain
    return [  # 6×3 matrix
        ϵ₁ ϵ₂+ϵ₃ 0
        ϵ₂ ϵ₁+ϵ₃ 0
        ϵ₃ ϵ₁+ϵ₂ 0
        0 0 ϵ₄
        0 0 ϵ₅
        0 0 ϵ₆
    ]
end
function recombine_strains(::Tetragonal, strain::EngineeringStrain)
    ϵ₁, ϵ₂, ϵ₃, ϵ₄, ϵ₅, ϵ₆ = strain
    # Tetragonal (I) class (c₁₆ = 0) is a special case of tetragonal (II) class
    return [  # 6×7 matrix
        ϵ₁ 0 ϵ₂ ϵ₃ ϵ₆ 0 0
        ϵ₂ 0 ϵ₁ ϵ₃ -ϵ₆ 0 0
        0 ϵ₃ 0 ϵ₁+ϵ₂ 0 0 0
        0 0 0 0 0 ϵ₄ 0
        0 0 0 0 0 ϵ₅ 0
        0 0 0 0 ϵ₁-ϵ₂ 0 ϵ₆
    ]
end
function recombine_strains(::Orthorhombic, strain::EngineeringStrain)
    ϵ₁, ϵ₂, ϵ₃, ϵ₄, ϵ₅, ϵ₆ = strain
    return [  # 6×9 matrix
        ϵ₁ 0 0 ϵ₂ ϵ₃ 0 0 0 0
        0 ϵ₂ 0 ϵ₁ 0 ϵ₃ 0 0 0
        0 0 ϵ₃ 0 ϵ₁ ϵ₂ 0 0 0
        0 0 0 0 0 0 ϵ₄ 0 0
        0 0 0 0 0 0 0 ϵ₅ 0
        0 0 0 0 0 0 0 0 ϵ₆
    ]
end
function recombine_strains(::Hexagonal, strain::EngineeringStrain)
    ϵ₁, ϵ₂, ϵ₃, ϵ₄, ϵ₅, ϵ₆ = strain
    return [  # 6×5 matrix
        ϵ₁ 0 ϵ₂ ϵ₃ 0
        ϵ₂ 0 ϵ₁ ϵ₃ 0
        0 ϵ₃ 0 ϵ₁+ϵ₂ 0
        0 0 0 0 ϵ₄
        0 0 0 0 ϵ₅
        ϵ₆/2 0 -ϵ₆/2 0 0
    ]
end
function recombine_strains(::Trigonal, strain::EngineeringStrain)
    ϵ₁, ϵ₂, ϵ₃, ϵ₄, ϵ₅, ϵ₆ = strain
    # Rhombohedral (I) class (c₁₅ = 0) is a special case of rhombohedral (II) class
    return [  # 6×7 matrix
        ϵ₁ 0 ϵ₂ ϵ₃ 0 ϵ₄ ϵ₅
        ϵ₂ 0 ϵ₁ ϵ₃ 0 -ϵ₄ -ϵ₅
        0 ϵ₃ 0 ϵ₁+ϵ₂ 0 0 0
        0 0 0 0 ϵ₄ ϵ₁-ϵ₂ -ϵ₆
        0 0 0 0 ϵ₅ ϵ₆ ϵ₁-ϵ₂
        ϵ₆/2 0 -ϵ₆/2 0 0 ϵ₅ -ϵ₄
    ]
end
function recombine_strains(::Monoclinic, strain::EngineeringStrain)  # Only standard orientation (diad ∥ x₂) is implemented
    ϵ₁, ϵ₂, ϵ₃, ϵ₄, ϵ₅, ϵ₆ = strain
    return [  # 6×13 matrix
        ϵ₁ 0 0 ϵ₂ ϵ₃ 0 0 0 0 ϵ₅ 0 0 0
        0 ϵ₂ 0 ϵ₁ 0 ϵ₃ 0 0 0 0 ϵ₅ 0 0
        0 0 ϵ₃ 0 ϵ₁ ϵ₂ 0 0 0 0 0 ϵ₅ 0
        0 0 0 0 0 0 ϵ₄ 0 0 0 0 0 ϵ₆
        0 0 0 0 0 0 0 ϵ₅ 0 ϵ₁ ϵ₂ ϵ₃ 0
        0 0 0 0 0 0 0 0 ϵ₆ 0 0 0 ϵ₄
    ]
end
function recombine_strains(::Triclinic, strain::EngineeringStrain)
    ϵ₁, ϵ₂, ϵ₃, ϵ₄, ϵ₅, ϵ₆ = strain
    return [  # 6×21 matrix
        ϵ₁ ϵ₂ ϵ₃ ϵ₄ ϵ₅ ϵ₆ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
        0 ϵ₁ 0 0 0 0 ϵ₂ ϵ₃ ϵ₄ ϵ₅ ϵ₆ 0 0 0 0 0 0 0 0 0 0
        0 0 ϵ₁ 0 0 0 0 ϵ₂ 0 0 0 ϵ₃ ϵ₄ ϵ₅ ϵ₆ 0 0 0 0 0 0
        0 0 0 ϵ₁ 0 0 0 0 ϵ₂ 0 0 0 ϵ₃ 0 0 ϵ₄ ϵ₅ ϵ₆ 0 0 0
        0 0 0 0 ϵ₁ 0 0 0 0 ϵ₂ 0 0 0 ϵ₃ 0 0 ϵ₄ 0 ϵ₅ ϵ₆ 0
        0 0 0 0 0 ϵ₁ 0 0 0 0 ϵ₂ 0 0 0 ϵ₃ 0 0 ϵ₄ 0 ϵ₅ ϵ₆
    ]
end
recombine_strains(system::CrystalSystem, strains::AbstractVector{<:EngineeringStrain}) =
    vcat((recombine_strains(system, strain) for strain in strains)...)

function reconstruct_cᵢⱼ(::Cubic, coefficients)
    c₁₁, c₁₂, c₄₄ = coefficients
    𝟎 = zero(c₁₁)
    return StiffnessMatrix(
        c₁₁,
        c₁₂,
        c₁₂,
        𝟎,
        𝟎,
        𝟎,
        c₁₁,
        c₁₂,
        𝟎,
        𝟎,
        𝟎,
        c₁₁,
        𝟎,
        𝟎,
        𝟎,
        c₄₄,
        𝟎,
        𝟎,
        c₄₄,
        𝟎,
        c₄₄,
    )
end
function reconstruct_cᵢⱼ(::Tetragonal, coefficients)
    c₁₁, c₃₃, c₁₂, c₁₃, c₁₆, c₄₄, c₆₆ = coefficients
    𝟎 = zero(c₁₁)
    return StiffnessMatrix(
        c₁₁,
        c₁₂,
        c₁₃,
        𝟎,
        𝟎,
        c₁₆,
        c₁₁,
        c₁₃,
        𝟎,
        𝟎,
        -c₁₆,
        c₃₃,
        𝟎,
        𝟎,
        𝟎,
        c₄₄,
        𝟎,
        𝟎,
        c₄₄,
        𝟎,
        c₆₆,
    )
end
function reconstruct_cᵢⱼ(::Orthorhombic, coefficients)
    c₁₁, c₂₂, c₃₃, c₁₂, c₁₃, c₂₃, c₄₄, c₅₅, c₆₆ = coefficients
    𝟎 = zero(c₁₁)
    return StiffnessMatrix(
        c₁₁,
        c₁₂,
        c₁₃,
        𝟎,
        𝟎,
        𝟎,
        c₂₂,
        c₂₃,
        𝟎,
        𝟎,
        𝟎,
        c₃₃,
        𝟎,
        𝟎,
        𝟎,
        c₄₄,
        𝟎,
        𝟎,
        c₅₅,
        𝟎,
        c₆₆,
    )
end
function reconstruct_cᵢⱼ(::Hexagonal, coefficients)
    c₁₁, c₃₃, c₁₂, c₁₃, c₄₄ = coefficients
    𝟎 = zero(c₁₁)
    return StiffnessMatrix(
        c₁₁,
        c₁₂,
        c₁₃,
        𝟎,
        𝟎,
        𝟎,
        c₁₁,
        c₁₃,
        𝟎,
        𝟎,
        𝟎,
        c₃₃,
        𝟎,
        𝟎,
        𝟎,
        c₄₄,
        𝟎,
        𝟎,
        c₄₄,
        𝟎,
        (c₁₁ - c₁₂) / 2,
    )
end
function reconstruct_cᵢⱼ(::Trigonal, coefficients)
    c₁₁, c₃₃, c₁₂, c₁₃, c₄₄, c₁₄, c₁₅ = coefficients
    𝟎 = zero(c₁₁)
    return StiffnessMatrix(
        c₁₁,
        c₁₂,
        c₁₃,
        c₁₄,
        c₁₅,
        𝟎,
        c₁₁,
        c₁₃,
        -c₁₄,
        -c₁₅,
        𝟎,
        c₃₃,
        𝟎,
        𝟎,
        𝟎,
        c₄₄,
        𝟎,
        -c₁₅,
        c₄₄,
        c₁₄,
        (c₁₁ - c₁₂) / 2,
    )
end
function reconstruct_cᵢⱼ(::Monoclinic, coefficients)
    c₁₁, c₂₂, c₃₃, c₁₂, c₁₃, c₂₃, c₄₄, c₅₅, c₆₆, c₁₅, c₂₅, c₃₅, c₄₆ = coefficients
    𝟎 = zero(c₁₁)
    return StiffnessMatrix(
        c₁₁,
        c₁₂,
        c₁₃,
        𝟎,
        c₁₅,
        𝟎,
        c₂₂,
        c₂₃,
        𝟎,
        c₂₅,
        𝟎,
        c₃₃,
        𝟎,
        c₃₅,
        𝟎,
        c₄₄,
        𝟎,
        c₄₆,
        c₅₅,
        𝟎,
        c₆₆,
    )
end
reconstruct_cᵢⱼ(::Triclinic, coefficients) = StiffnessMatrix(coefficients...)

function recombine_stresses(::Cubic, stress::EngineeringStress)
    σ₁, σ₂, σ₃, σ₄, σ₅, σ₆ = stress
    return [  # 6×3 matrix
        σ₁ σ₂+σ₃ 0
        σ₂ σ₁+σ₃ 0
        σ₃ σ₁+σ₂ 0
        0 0 σ₄
        0 0 σ₅
        0 0 σ₆
    ]
end
recombine_stresses(system::CrystalSystem, stresses::AbstractVector{<:EngineeringStress}) =
    vcat((recombine_stresses(system, stress) for stress in stresses)...)

function reconstruct_sᵢⱼ(::Cubic, coefficients)
    s₁₁, s₁₂, s₄₄ = coefficients
    𝟎 = zero(s₁₁)
    return StiffnessMatrix(
        s₁₁,
        s₁₂,
        s₁₂,
        𝟎,
        𝟎,
        𝟎,
        s₁₁,
        s₁₂,
        𝟎,
        𝟎,
        𝟎,
        s₁₁,
        𝟎,
        𝟎,
        𝟎,
        s₄₄,
        𝟎,
        𝟎,
        s₄₄,
        𝟎,
        s₄₄,
    )
end

function solve_elastic_constants(
    system::CrystalSystem,
    strains::AbstractVector{<:EngineeringStrain},
    stresses::AbstractVector{<:EngineeringStress},
)
    if length(strains) != length(stresses)
        throw(
            DimensionMismatch(
                "the number of strains and the number of stresses must match!",
            ),
        )
    end
    n = minimal_npairs(system)
    if length(strains) < n
        throw(ArgumentError("the number of strains/stresses must be at least $n."))
    end
    σ = vcat(stresses...)  # Length 6n vector, n = length(strains) = length(stresses)
    ε = recombine_strains(system, strains)  # Size 6n×N matrix, N = # independent coefficients
    coefficients = ε \ σ  # Length N vector
    return reconstruct_cᵢⱼ(system, coefficients)
end
function solve_elastic_constants(
    system::CrystalSystem,
    stresses::AbstractVector{<:EngineeringStress},
    strains::AbstractVector{<:EngineeringStrain},
)
    if length(strains) != length(stresses)
        throw(
            DimensionMismatch(
                "the number of strains and the number of stresses must match!",
            ),
        )
    end
    n = minimal_npairs(system)
    if length(strains) < n
        throw(ArgumentError("the number of strains/stresses must be at least $n."))
    end
    ε = vcat(strains...)
    σ = recombine_stresses(system, stresses)
    coefficients = σ \ ε
    return reconstruct_sᵢⱼ(system, coefficients)
end
function solve_elastic_constants(
    system::CrystalSystem,
    strains::AbstractVector{<:TensorStrain},
    stresses::AbstractVector{<:TensorStress},
)
    cᵢⱼ = solve_elastic_constants(
        system,
        EngineeringStrain.(strains),
        EngineeringStress.(stresses),
    )
    return StiffnessTensor(cᵢⱼ)
end
function solve_elastic_constants(
    system::CrystalSystem,
    stresses::AbstractVector{<:TensorStress},
    strains::AbstractVector{<:TensorStrain},
)
    sᵢⱼ = solve_elastic_constants(
        system,
        EngineeringStrain.(strains),
        EngineeringStress.(stresses),
    )
    return ComplianceTensor(sᵢⱼ)
end
solve_elastic_constants(strains_or_stresses, stresses_or_strains) =
    solve_elastic_constants(Triclinic(), strains_or_stresses, stresses_or_strains)

solve_stiffnesses(
    system::CrystalSystem,
    strains::AbstractVector{<:EngineeringStrain},
    stresses::AbstractVector{<:EngineeringStress},
) = solve_elastic_constants(system, strains, stresses)
solve_stiffnesses(
    system::CrystalSystem,
    strains::AbstractVector{<:TensorStrain},
    stresses::AbstractVector{<:TensorStress},
) = solve_elastic_constants(system, strains, stresses)
solve_stiffnesses(strains, stresses) = solve_stiffnesses(Triclinic(), strains, stresses)

solve_compliances(
    system::CrystalSystem,
    stresses::AbstractVector{<:EngineeringStress},
    strains::AbstractVector{<:EngineeringStrain},
) = solve_elastic_constants(system, stresses, strains)
solve_compliances(
    system::CrystalSystem,
    stresses::AbstractVector{<:TensorStress},
    strains::AbstractVector{<:TensorStrain},
) = solve_elastic_constants(system, stresses, strains)
solve_compliances(stresses, strains) = solve_compliances(Triclinic(), stresses, strains)

minimal_npairs(::Cubic) = 1
minimal_npairs(::Hexagonal) = 2
minimal_npairs(::Trigonal) = 2
minimal_npairs(::Tetragonal) = 2
minimal_npairs(::Orthorhombic) = 3
minimal_npairs(::Monoclinic) = 5
minimal_npairs(::Triclinic) = 6


text:
import Base: getindex

include("point.jl")

abstract type GeoObject end

struct Construction{T <: Number}
    objects::Vector{GeoObject}
end

Construction{T}() where T = Construction{T}(Vector{GeoObject}())

struct GeoPointer{T}
    index::Int
end

function getindex(construction::Construction, pointer::GeoPointer{T}) where T
    return state(construction.objects[pointer.index])::T
end

struct FreePoint{T} <: GeoObject
    position::Point{T}
end

state(free::FreePoint) = free.position

function new_free_point!(construction::Construction{T}, position::Point) where T
    free = FreePoint(convert(Point{T}, position))
    push!(construction.objects, free)

    return GeoPointer{Point{T}}(length(construction.objects))
end

update(::Construction, p::FreePoint) = p

update_param(::Construction, ::FreePoint, p::Point) = FreePoint(p)

function interpolate(::Construction, t::Real, point::FreePoint, target::Point)
    t2 = t * t
    dt = 0.5 / (1 + t2)
    tc = complex((2 * t) * dt + 0.5, (1 - t2) * dt)

    source = point.position

    return source + tc * (target - source)
end

struct Distance{T,S <: Real} <: GeoObject
    a::GeoPointer{Point{T}}
    b::GeoPointer{Point{T}}

    d::S
end

state(dist::Distance) = dist.d

function new_distance!(construction::Construction{T}, a::GeoPointer{Point{T}}, b::GeoPointer{Point{T}}) where T
    dist = Distance(a, b, distance(construction[a], construction[b]))

    push!(construction.objects, dist)

    return GeoPointer{real(T)}(length(construction.objects))
end

update(construction::Construction, dist::Distance) = Distance(dist.a, dist.b, distance(construction[dist.a], construction[dist.b]))

struct Join{T} <: GeoObject
    a::GeoPointer{Point{T}}
    b::GeoPointer{Point{T}}

    line::Point{T}
end

state(j::Join) = j.line

function new_join!(construction::Construction{T}, a::GeoPointer{Point{T}}, b::GeoPointer{Point{T}}) where T
    join = Join(a, b, cross(construction[a], construction[b]))
    push!(construction.objects, join)

    return GeoPointer{Point{T}}(length(construction.objects))
end

update(construction::Construction, j::Join) = Join(j.a, j.b, cross(construction[j.a], construction[j.b]))

struct IntersectionLineLine{T} <: GeoObject
    l1::GeoPointer{Point{T}}
    l2::GeoPointer{Point{T}}

    point::Point{T}
end

state(ll::IntersectionLineLine) = ll.point

function new_intersect_line_line!(construction::Construction{T}, l1::GeoPointer{Point{T}}, l2::GeoPointer{Point{T}}) where T
    ll = IntersectionLineLine(l1, l2, cross(construction[l1], construction[l2]))

    push!(construction.objects, ll)

    return GeoPointer{Point{T}}(length(construction.objects))
end

function update(construction::Construction, ll::IntersectionLineLine)
    return IntersectionLineLine(ll.l1, ll.l2, cross(construction[ll.l1], construction[ll.l2]))
end

struct IntersectionLineCircle{T} <: GeoObject
    line::GeoPointer{Point{T}}
    circle::GeoPointer{CircleMatrix{T}}

    p1::Point{T}
    p2::Point{T}
end

state(lc::IntersectionLineCircle) = (lc.p1, lc.p2)

function new_intersect_line_circle!(construction::Construction{T}, line::GeoPointer{Point{T}}, circle::GeoPointer{CircleMatrix{T}}) where T
    lc = IntersectionLineCircle(line, circle, intersect_line_circle(construction[line], construction[circle])...)

    push!(construction.objects, lc)

    return GeoPointer{Tuple{Point{T},Point{T}}}(length(construction.objects))
end

function update(construction::Construction, lc::IntersectionLineCircle)
    result = tracing2((lc.p1, lc.p2), intersect_line_circle(construction[lc.line], construction[lc.circle]))

    if result === nothing
        return nothing
    end

    return IntersectionLineCircle(lc.line, lc.circle, result...)
end

struct PointOnLine{T} <: GeoObject
    line::GeoPointer{Point{T}}

    point::Point{T}
end

state(pol::PointOnLine) = pol.point

function new_point_on_line!(construction::Construction{T}, line::GeoPointer{Point{T}}, nearest::Point) where T
    nearest = convert(Point{T}, nearest)

    pol = PointOnLine(line, project_point_to_line(nearest, construction[line]))
    push!(construction.objects, pol)

    return GeoPointer{Point{T}}(length(construction.objects))
end

function update(construction::Construction, pol::PointOnLine, orig_construction::Construction, orig_pol::PointOnLine)
    new_line = construction[pol.line]
    old_point = pol.point

    real_line = orig_construction[orig_pol.line]
    real_point = orig_pol.point

    center = cross(real_line, new_line)

    # TODO: handle line stay almost the same
    @assert abs(center) > 3eps(Float64)

    circle = circle_from_mp(center, real_point)

    candidates = intersect_line_circle(new_line, circle)

    old_antipode = point_reflection(center, old_point)

    result = tracing2((old_point, old_antipode), candidates)

    if result === nothing
        return nothing
    end

    return PointOnLine(pol.line, result[1])
end

function update_param(construction::Construction, pol::PointOnLine, param::Point)
    return PointOnLine(pol.line, param)
end

function interpolate(construction::Construction, t::Real, pol::PointOnLine, target::Point)
    t2 = t * t
    dt = 0.5 / (1 + t2)
    tc = complex((2 * t) * dt + 0.5, (1 - t2) * dt)

    target = project_point_to_line(target, construction[pol.line])
    source = pol.point

    return source + tc * (target - source)
end

struct CircleMP{T} <: GeoObject
    m::GeoPointer{Point{T}}
    p::GeoPointer{Point{T}}

    matrix::CircleMatrix{T}
end

state(circle::CircleMP) = circle.matrix

function new_circle_mp!(construction::Construction{T}, m::GeoPointer{Point{T}}, p::GeoPointer{Point{T}}) where T
    circle = CircleMP(m, p, circle_from_mp(construction[m], construction[p]))
    push!(construction.objects, circle)

    return GeoPointer{CircleMatrix{T}}(length(construction.objects))
end

update(construction::Construction, c::CircleMP) = CircleMP(c.m, c.p, circle_from_mp(construction[c.m], construction[c.p]))

struct CircleMr{T,S <: Real} <: GeoObject
    m::GeoPointer{Point{T}}
    r::GeoPointer{S}

    matrix::CircleMatrix{T}
end

state(circle::CircleMr) = circle.matrix

function new_circle_mr!(construction::Construction{T}, m::GeoPointer{Point{T}}, r::GeoPointer{S}) where {T,S <: Real}
    @assert S == real(T)

    circle = CircleMr(m, r, circle_from_mr(construction[m], construction[r]))
    push!(construction.objects, circle)

    return GeoPointer{CircleMatrix{T}}(length(construction.objects))
end

update(construction::Construction, c::CircleMr) = CircleMr(c.m, c.r, circle_from_mr(construction[c.m], construction[c.r]))

struct IntersectionCircleCircle{T} <: GeoObject
    c1::GeoPointer{CircleMatrix{T}}
    c2::GeoPointer{CircleMatrix{T}}

    p1::Point{T}
    p2::Point{T}
end

state(cc::IntersectionCircleCircle) = (cc.p1, cc.p2)

function new_intersect_circle_circle!(construction::Construction{T}, c1::GeoPointer{CircleMatrix{T}}, c2::GeoPointer{CircleMatrix{T}}) where T
    cc = IntersectionCircleCircle(c1, c2, intersect_circle_circle(construction[c1], construction[c2])...)
    push!(construction.objects, cc)

    return GeoPointer{Tuple{Point{T},Point{T}}}(length(construction.objects))
end

function update(construction::Construction, cc::IntersectionCircleCircle)
    result = tracing2((cc.p1, cc.p2), intersect_circle_circle(construction[cc.c1], construction[cc.c2]))

    if result === nothing
        return nothing
    end

    return IntersectionCircleCircle(cc.c1, cc.c2, result...)
end

struct SelectPoint{T} <: GeoObject
    collection::GeoPointer{Tuple{Point{T},Point{T}}}

    index::Int
    p::Point{T}
end

state(select::SelectPoint) = select.p

function new_select_point!(construction::Construction{T}, collection::GeoPointer{Tuple{Point{T},Point{T}}}, index::Int) where T
    select = SelectPoint(collection, index, construction[collection][index])
    push!(construction.objects, select)

    return GeoPointer{Point{T}}(length(construction.objects))
end

function update(construction::Construction, select::SelectPoint)
    return SelectPoint(select.collection, select.index, construction[select.collection][select.index])
end

struct PointOnCircle{T} <: GeoObject
    circle::GeoPointer{CircleMatrix{T}}

    param::Point{T} # far point polar to the diameter through the point
    point::Point{T}
    antipodal::Point{T}
end

state(poc::PointOnCircle) = poc.point

function new_point_on_circle!(construction::Construction{T}, circle::GeoPointer{CircleMatrix{T}}, nearest::Point) where T
    nearest = convert(Point{T}, nearest)

    circle_mat = construction[circle]

    pos = normalize_z(nearest)
    mid = normalize_z(center_of_circle(circle_mat))

    dir = pos - mid

    param = Point(dir.y, -dir.x, 0)

    diameter = cross(pos, mid)

    candidates = intersect_line_circle(diameter, circle_mat)

    d1 = projective_dist_min_scal(pos, candidates[1])
    d2 = projective_dist_min_scal(pos, candidates[2])

    if d2 < d1
        point = candidates[2]
        antipodal = candidates[1]
    else
        point = candidates[1]
        antipodal = candidates[2]
    end

    poc = PointOnCircle(circle, param, point, antipodal)
    push!(construction.objects, poc)

    return GeoPointer{Point{T}}(length(construction.objects))
end

function update_param(construction::Construction, poc::PointOnCircle, param::Point)
    circle_mat = construction[poc.circle]

    diameter = circle_mat * param

    candidates = intersect_line_circle(diameter, circle_mat)

    result = tracing2((poc.point, poc.antipodal), candidates)

    if result === nothing
        return nothing
    end

    return PointOnCircle(poc.circle, param, result...)
end

update(construction::Construction, poc::PointOnCircle) = update_param(construction, poc, poc.param)

function poc_param_from_nearest(poc::PointOnCircle, circle::CircleMatrix, nearest::Point)
    mid = normalize_z(center_of_circle(circle))
    dir = nearest - mid

    oldparam = poc.param
    oldpos = normalize_z(poc.point)
    olddir = oldpos - mid

    oldsign = oldparam.x * olddir.y - oldparam.y * olddir.x

    if real(oldsign) < 0
        dir = -dir
    end

    return Point(dir.y, -dir.x, 0)
end

function interpolate(construction::Construction, tr::Real, poc::PointOnCircle, target::Point)
    src = poc.param
    dst = poc_param_from_nearest(poc, construction[poc.circle], target)

    sp = dot(src, dst)

    if real(sp) >= 0
        t2 = tr * tr
        dt = 0.5 / (1 + t2)
        tc = complex((2 * tr) * dt + 0.5, (1 - t2) * dt)

        return src + tc * (dst - src)
    end

    mid = Point(src.y - dst.y, dst.x - src.x, 0)

    sp = dot(src, mid)

    if real(sp) < 0
        mid = -mid
    end

    if tr < 0
        tr = 2 * tr + 1
        t2 = tr * tr
        dt = 0.25 / (1 + t2)
        tc = complex((2 * tr) * dt + 0.25, (1 - t2) * dt)
    else
        tr = 2 * tr - 1
        t2 = tr * tr
        dt = 0.25 / (1 + t2)
        tc = complex((2 * tr) * dt + 0.75, (1 - t2) * dt)
    end

    uc = 1 - tc

    tc2 = tc * tc

    uc2 = uc * uc

    tuc = tc * uc

    return uc2 * src + tuc * mid + tc2 * dst
end

function move_to!(construction::Construction{T}, pointer::GeoPointer{Point{T}}, target::Point) where T
    target = convert(Point{T}, target)

    # we are trying to reach a new good state
    # thus we start by working on a copy of our Construction
    # and apply the changes only after we have successfully reached a new good state
    sandbox = Construction{T}(copy(construction.objects))

    # also remember last good state
    lastgood = Construction{T}(copy(construction.objects))

    last = -1.0
    step = 0.9

    t = last + step

    mover_index = pointer.index
    mover = construction.objects[mover_index]

    while last != t
        refine = false

        result = update_param(sandbox, sandbox.objects[mover_index], interpolate(construction, t, mover, target))

        if result !== nothing
            sandbox.objects[mover_index] = result

            # update all dependent objects
            for i = mover_index+1:length(sandbox.objects)
                if sandbox.objects[i] isa PointOnLine # special snowflake...
                    result = update(sandbox, sandbox.objects[i], construction, construction.objects[i])
                else
                    result = update(sandbox, sandbox.objects[i])
                end

                if result === nothing
                    refine = true
                    break
                end

                sandbox.objects[i] = result
            end
        else
            refine = true
        end

        if refine
            # this is bad. we could not update position for this object
            # throw current sandbox away, create new sandbox from lastgood and refine
            sandbox = Construction{T}(copy(lastgood.objects))

            step *= 0.5 # reduce step size
            t = last + step

            continue
        end

        last = t # successfully traced up to t
        step *= 1.25
        t += step

        if t >= 1
            t = 1.0
        end

        lastgood = sandbox
        sandbox = Construction{T}(copy(lastgood.objects))
    end

    if t != 1
        error("could not trace to end")
    end

    copyto!(construction.objects, lastgood.objects)
    return nothing
end