text:
#Reference:
#须寅, 龙驭球. 采用广义协调条件构造具有旋转自由度的四边形膜元[J]. 工程力学, 1993, 10(3): 27―36.
function K_GQ12(elm::Quad)::Matrix{Float64}
    E₀,ν₀=elm.material.E,elm.material.ν
    center=elm.center
    t=elm.t
    T=elm.T[1:3,1:3]
    x₁,y₁,z₁=T*(elm.node1.loc.-center)
    x₂,y₂,z₂=T*(elm.node2.loc.-center)
    x₃,y₃,z₃=T*(elm.node3.loc.-center)
    x₄,y₄,z₄=T*(elm.node4.loc.-center)
    K=Matrix{Float64}(undef,12,12)
    J=Matrix{Float64}(undef,2,2)
    function BtDB(x)
            ξ,η=x[1],x[2]
            K[1,1]=0.0625*(1 - η)^2*E₀/(1 - ν₀^2) + 0.03125*(1 - ξ)^2*E₀*(1 - ν₀)/(1 - ν₀^2)
            K[1,2]=K[2,1]=0.0625*(1 - ξ)*(1 - η)*E₀*ν₀/(1 - ν₀^2) + 0.03125*(1 - ξ)*(1 - η)*E₀*(1 - ν₀)/(1 - ν₀^2)
            K[1,3]=K[3,1]=-0.03125*(1 - η)*E₀*((1 - η^2)*(0.25*(-y₁ - y₂ + y₃ + y₄) - 0.25*(y₁ - y₂ + y₃ - y₄)) + 2*ξ*(1 - η)*(0.25*(-y₁ + y₂ + y₃ - y₄) - 0.25*(y₁ - y₂ + y₃ - y₄)))/(1 - ν₀^2) + 0.03125*(1 - η)*E₀*ν₀*((1 - ξ^2)*(-0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ + x₂ + x₃ - x₄)) + 2*η*(1 - ξ)*(-0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ - x₂ + x₃ + x₄)))/(1 - ν₀^2) - 0.125*(1 - ξ)*E₀*(1 - ν₀)*(-0.125*((1 - η^2)*(-0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ - x₂ + x₃ + x₄)) + 2*ξ*(1 - η)*(-0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ + x₂ + x₃ - x₄))) + 0.125*((1 - ξ^2)*(0.25*(-y₁ + y₂ + y₃ - y₄) - 0.25*(y₁ - y₂ + y₃ - y₄)) + 2*η*(1 - ξ)*(0.25*(-y₁ - y₂ + y₃ + y₄) - 0.25*(y₁ - y₂ + y₃ - y₄))))/(1 - ν₀^2)
            K[1,4]=K[4,1]=-0.0625*(1 - η)^2*E₀/(1 - ν₀^2) + 0.03125*(1 - ξ)*(1 + ξ)*E₀*(1 - ν₀)/(1 - ν₀^2)
            K[1,5]=K[5,1]=-0.03125*(1 - ξ)*(1 - η)*E₀*(1 - ν₀)/(1 - ν₀^2) + 0.0625*(1 + ξ)*(1 - η)*E₀*ν₀/(1 - ν₀^2)
            K[1,6]=K[6,1]=-0.03125*(1 - η)*E₀*(-(0.25*(-y₁ - y₂ + y₃ + y₄) + 0.25*(y₁ - y₂ + y₃ - y₄))*(1 - η^2) - 2*ξ*(1 - η)*(0.25*(-y₁ + y₂ + y₃ - y₄) - 0.25*(y₁ - y₂ + y₃ - y₄)))/(1 - ν₀^2) + 0.03125*(1 - η)*E₀*ν₀*(-(1 - ξ^2)*(-0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ + x₂ + x₃ - x₄)) + 2*(0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ - x₂ + x₃ + x₄))*η*(1 + ξ))/(1 - ν₀^2) - 0.125*(1 - ξ)*E₀*(1 - ν₀)*(0.125*(-(1 - ξ^2)*(0.25*(-y₁ + y₂ + y₃ - y₄) - 0.25*(y₁ - y₂ + y₃ - y₄)) + 2*(0.25*(-y₁ - y₂ + y₃ + y₄) + 0.25*(y₁ - y₂ + y₃ - y₄))*η*(1 + ξ)) - 0.125*(-(0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ - x₂ + x₃ + x₄))*(1 - η^2) - 2*ξ*(1 - η)*(-0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ + x₂ + x₃ - x₄))))/(1 - ν₀^2)
            K[1,7]=K[7,1]=-0.0625*(1 - η)*(1 + η)*E₀/(1 - ν₀^2) - 0.03125*(1 - ξ)*(1 + ξ)*E₀*(1 - ν₀)/(1 - ν₀^2)
            K[1,8]=K[8,1]=-0.03125*(1 - ξ)*(1 + η)*E₀*(1 - ν₀)/(1 - ν₀^2) - 0.0625*(1 + ξ)*(1 - η)*E₀*ν₀/(1 - ν₀^2)
            K[1,9]=K[9,1]=-0.03125*(1 - η)*E₀*((0.25*(-y₁ - y₂ + y₃ + y₄) + 0.25*(y₁ - y₂ + y₃ - y₄))*(1 - η^2) - 2*(0.25*(-y₁ + y₂ + y₃ - y₄) + 0.25*(y₁ - y₂ + y₃ - y₄))*ξ*(1 + η))/(1 - ν₀^2) + 0.03125*(1 - η)*E₀*ν₀*((0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ + x₂ + x₃ - x₄))*(1 - ξ^2) - 2*(0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ - x₂ + x₃ + x₄))*η*(1 + ξ))/(1 - ν₀^2) - 0.125*(1 - ξ)*E₀*(1 - ν₀)*(-0.125*((0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ - x₂ + x₃ + x₄))*(1 - η^2) - 2*(0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ + x₂ + x₃ - x₄))*ξ*(1 + η)) + 0.125*((0.25*(-y₁ + y₂ + y₃ - y₄) + 0.25*(y₁ - y₂ + y₃ - y₄))*(1 - ξ^2) - 2*(0.25*(-y₁ - y₂ + y₃ + y₄) + 0.25*(y₁ - y₂ + y₃ - y₄))*η*(1 + ξ)))/(1 - ν₀^2)
            K[1,10]=K[10,1]=0.0625*(1 - η)*(1 + η)*E₀/(1 - ν₀^2) - 0.03125*(1 - ξ)^2*E₀*(1 - ν₀)/(1 - ν₀^2)
            K[1,11]=K[11,1]=-0.0625*(1 - ξ)*(1 - η)*E₀*ν₀/(1 - ν₀^2) + 0.03125*(1 - ξ)*(1 + η)*E₀*(1 - ν₀)/(1 - ν₀^2)
            K[1,12]=K[12,1]=-0.03125*(1 - η)*E₀*(-(1 - η^2)*(0.25*(-y₁ - y₂ + y₃ + y₄) - 0.25*(y₁ - y₂ + y₃ - y₄)) + 2*(0.25*(-y₁ + y₂ + y₃ - y₄) + 0.25*(y₁ - y₂ + y₃ - y₄))*ξ*(1 + η))/(1 - ν₀^2) + 0.03125*(1 - η)*E₀*ν₀*(-(0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ + x₂ + x₃ - x₄))*(1 - ξ^2) - 2*η*(1 - ξ)*(-0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ - x₂ + x₃ + x₄)))/(1 - ν₀^2) - 0.125*(1 - ξ)*E₀*(1 - ν₀)*(-0.125*(-(1 - η^2)*(-0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ - x₂ + x₃ + x₄)) + 2*(0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ + x₂ + x₃ - x₄))*ξ*(1 + η)) + 0.125*(-(0.25*(-y₁ + y₂ + y₃ - y₄) + 0.25*(y₁ - y₂ + y₃ - y₄))*(1 - ξ^2) - 2*η*(1 - ξ)*(0.25*(-y₁ - y₂ + y₃ + y₄) - 0.25*(y₁ - y₂ + y₃ - y₄))))/(1 - ν₀^2)
            K[2,2]=0.0625*(1 - ξ)^2*E₀/(1 - ν₀^2) + 0.03125*(1 - η)^2*E₀*(1 - ν₀)/(1 - ν₀^2)
            K[2,3]=K[3,2]=0.03125*(1 - ξ)*E₀*((1 - ξ^2)*(-0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ + x₂ + x₃ - x₄)) + 2*η*(1 - ξ)*(-0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ - x₂ + x₃ + x₄)))/(1 - ν₀^2) - 0.125*(1 - η)*E₀*(1 - ν₀)*(-0.125*((1 - η^2)*(-0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ - x₂ + x₃ + x₄)) + 2*ξ*(1 - η)*(-0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ + x₂ + x₃ - x₄))) + 0.125*((1 - ξ^2)*(0.25*(-y₁ + y₂ + y₃ - y₄) - 0.25*(y₁ - y₂ + y₃ - y₄)) + 2*η*(1 - ξ)*(0.25*(-y₁ - y₂ + y₃ + y₄) - 0.25*(y₁ - y₂ + y₃ - y₄))))/(1 - ν₀^2) - 0.03125*(1 - ξ)*E₀*ν₀*((1 - η^2)*(0.25*(-y₁ - y₂ + y₃ + y₄) - 0.25*(y₁ - y₂ + y₃ - y₄)) + 2*ξ*(1 - η)*(0.25*(-y₁ + y₂ + y₃ - y₄) - 0.25*(y₁ - y₂ + y₃ - y₄)))/(1 - ν₀^2)
            K[2,4]=K[4,2]=-0.0625*(1 - ξ)*(1 - η)*E₀*ν₀/(1 - ν₀^2) + 0.03125*(1 + ξ)*(1 - η)*E₀*(1 - ν₀)/(1 - ν₀^2)
            K[2,5]=K[5,2]=-0.03125*(1 - η)^2*E₀*(1 - ν₀)/(1 - ν₀^2) + 0.0625*(1 - ξ)*(1 + ξ)*E₀/(1 - ν₀^2)
            K[2,6]=K[6,2]=0.03125*(1 - ξ)*E₀*(-(1 - ξ^2)*(-0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ + x₂ + x₃ - x₄)) + 2*(0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ - x₂ + x₃ + x₄))*η*(1 + ξ))/(1 - ν₀^2) - 0.125*(1 - η)*E₀*(1 - ν₀)*(0.125*(-(1 - ξ^2)*(0.25*(-y₁ + y₂ + y₃ - y₄) - 0.25*(y₁ - y₂ + y₃ - y₄)) + 2*(0.25*(-y₁ - y₂ + y₃ + y₄) + 0.25*(y₁ - y₂ + y₃ - y₄))*η*(1 + ξ)) - 0.125*(-(0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ - x₂ + x₃ + x₄))*(1 - η^2) - 2*ξ*(1 - η)*(-0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ + x₂ + x₃ - x₄))))/(1 - ν₀^2) - 0.03125*(1 - ξ)*E₀*ν₀*(-(0.25*(-y₁ - y₂ + y₃ + y₄) + 0.25*(y₁ - y₂ + y₃ - y₄))*(1 - η^2) - 2*ξ*(1 - η)*(0.25*(-y₁ + y₂ + y₃ - y₄) - 0.25*(y₁ - y₂ + y₃ - y₄)))/(1 - ν₀^2)
            K[2,7]=K[7,2]=-0.0625*(1 - ξ)*(1 + η)*E₀*ν₀/(1 - ν₀^2) - 0.03125*(1 + ξ)*(1 - η)*E₀*(1 - ν₀)/(1 - ν₀^2)
            K[2,8]=K[8,2]=-0.0625*(1 - ξ)*(1 + ξ)*E₀/(1 - ν₀^2) - 0.03125*(1 - η)*(1 + η)*E₀*(1 - ν₀)/(1 - ν₀^2)
            K[2,9]=K[9,2]=0.03125*(1 - ξ)*E₀*((0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ + x₂ + x₃ - x₄))*(1 - ξ^2) - 2*(0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ - x₂ + x₃ + x₄))*η*(1 + ξ))/(1 - ν₀^2) - 0.125*(1 - η)*E₀*(1 - ν₀)*(-0.125*((0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ - x₂ + x₃ + x₄))*(1 - η^2) - 2*(0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ + x₂ + x₃ - x₄))*ξ*(1 + η)) + 0.125*((0.25*(-y₁ + y₂ + y₃ - y₄) + 0.25*(y₁ - y₂ + y₃ - y₄))*(1 - ξ^2) - 2*(0.25*(-y₁ - y₂ + y₃ + y₄) + 0.25*(y₁ - y₂ + y₃ - y₄))*η*(1 + ξ)))/(1 - ν₀^2) - 0.03125*(1 - ξ)*E₀*ν₀*((0.25*(-y₁ - y₂ + y₃ + y₄) + 0.25*(y₁ - y₂ + y₃ - y₄))*(1 - η^2) - 2*(0.25*(-y₁ + y₂ + y₃ - y₄) + 0.25*(y₁ - y₂ + y₃ - y₄))*ξ*(1 + η))/(1 - ν₀^2)
            K[2,10]=K[10,2]=-0.03125*(1 - ξ)*(1 - η)*E₀*(1 - ν₀)/(1 - ν₀^2) + 0.0625*(1 - ξ)*(1 + η)*E₀*ν₀/(1 - ν₀^2)
            K[2,11]=K[11,2]=-0.0625*(1 - ξ)^2*E₀/(1 - ν₀^2) + 0.03125*(1 - η)*(1 + η)*E₀*(1 - ν₀)/(1 - ν₀^2)
            K[2,12]=K[12,2]=0.03125*(1 - ξ)*E₀*(-(0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ + x₂ + x₃ - x₄))*(1 - ξ^2) - 2*η*(1 - ξ)*(-0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ - x₂ + x₃ + x₄)))/(1 - ν₀^2) - 0.125*(1 - η)*E₀*(1 - ν₀)*(-0.125*(-(1 - η^2)*(-0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ - x₂ + x₃ + x₄)) + 2*(0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ + x₂ + x₃ - x₄))*ξ*(1 + η)) + 0.125*(-(0.25*(-y₁ + y₂ + y₃ - y₄) + 0.25*(y₁ - y₂ + y₃ - y₄))*(1 - ξ^2) - 2*η*(1 - ξ)*(0.25*(-y₁ - y₂ + y₃ + y₄) - 0.25*(y₁ - y₂ + y₃ - y₄))))/(1 - ν₀^2) - 0.03125*(1 - ξ)*E₀*ν₀*(-(1 - η^2)*(0.25*(-y₁ - y₂ + y₃ + y₄) - 0.25*(y₁ - y₂ + y₃ - y₄)) + 2*(0.25*(-y₁ + y₂ + y₃ - y₄) + 0.25*(y₁ - y₂ + y₃ - y₄))*ξ*(1 + η))/(1 - ν₀^2)
            K[3,3]=0.125*((1 - η^2)*(0.25*(-y₁ - y₂ + y₃ + y₄) - 0.25*(y₁ - y₂ + y₃ - y₄)) + 2*ξ*(1 - η)*(0.25*(-y₁ + y₂ + y₃ - y₄) - 0.25*(y₁ - y₂ + y₃ - y₄)))*(0.125*E₀*((1 - η^2)*(0.25*(-y₁ - y₂ + y₃ + y₄) - 0.25*(y₁ - y₂ + y₃ - y₄)) + 2*ξ*(1 - η)*(0.25*(-y₁ + y₂ + y₃ - y₄) - 0.25*(y₁ - y₂ + y₃ - y₄)))/(1 - ν₀^2) - 0.125*E₀*ν₀*((1 - ξ^2)*(-0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ + x₂ + x₃ - x₄)) + 2*η*(1 - ξ)*(-0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ - x₂ + x₃ + x₄)))/(1 - ν₀^2)) - 0.125*((1 - ξ^2)*(-0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ + x₂ + x₃ - x₄)) + 2*η*(1 - ξ)*(-0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ - x₂ + x₃ + x₄)))*(-0.125*E₀*((1 - ξ^2)*(-0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ + x₂ + x₃ - x₄)) + 2*η*(1 - ξ)*(-0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ - x₂ + x₃ + x₄)))/(1 - ν₀^2) + 0.125*E₀*ν₀*((1 - η^2)*(0.25*(-y₁ - y₂ + y₃ + y₄) - 0.25*(y₁ - y₂ + y₃ - y₄)) + 2*ξ*(1 - η)*(0.25*(-y₁ + y₂ + y₃ - y₄) - 0.25*(y₁ - y₂ + y₃ - y₄)))/(1 - ν₀^2)) + (1/2)*E₀*(1 - ν₀)*(-0.125*((1 - η^2)*(-0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ - x₂ + x₃ + x₄)) + 2*ξ*(1 - η)*(-0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ + x₂ + x₃ - x₄))) + 0.125*((1 - ξ^2)*(0.25*(-y₁ + y₂ + y₃ - y₄) - 0.25*(y₁ - y₂ + y₃ - y₄)) + 2*η*(1 - ξ)*(0.25*(-y₁ - y₂ + y₃ + y₄) - 0.25*(y₁ - y₂ + y₃ - y₄))))^2/(1 - ν₀^2)
            K[3,4]=K[4,3]=0.25*(1 - η)*(0.125*E₀*((1 - η^2)*(0.25*(-y₁ - y₂ + y₃ + y₄) - 0.25*(y₁ - y₂ + y₃ - y₄)) + 2*ξ*(1 - η)*(0.25*(-y₁ + y₂ + y₃ - y₄) - 0.25*(y₁ - y₂ + y₃ - y₄)))/(1 - ν₀^2) - 0.125*E₀*ν₀*((1 - ξ^2)*(-0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ + x₂ + x₃ - x₄)) + 2*η*(1 - ξ)*(-0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ - x₂ + x₃ + x₄)))/(1 - ν₀^2)) - 0.125*(1 + ξ)*E₀*(1 - ν₀)*(-0.125*((1 - η^2)*(-0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ - x₂ + x₃ + x₄)) + 2*ξ*(1 - η)*(-0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ + x₂ + x₃ - x₄))) + 0.125*((1 - ξ^2)*(0.25*(-y₁ + y₂ + y₃ - y₄) - 0.25*(y₁ - y₂ + y₃ - y₄)) + 2*η*(1 - ξ)*(0.25*(-y₁ - y₂ + y₃ + y₄) - 0.25*(y₁ - y₂ + y₃ - y₄))))/(1 - ν₀^2)
            K[3,5]=K[5,3]=-0.25*(1 + ξ)*(-0.125*E₀*((1 - ξ^2)*(-0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ + x₂ + x₃ - x₄)) + 2*η*(1 - ξ)*(-0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ - x₂ + x₃ + x₄)))/(1 - ν₀^2) + 0.125*E₀*ν₀*((1 - η^2)*(0.25*(-y₁ - y₂ + y₃ + y₄) - 0.25*(y₁ - y₂ + y₃ - y₄)) + 2*ξ*(1 - η)*(0.25*(-y₁ + y₂ + y₃ - y₄) - 0.25*(y₁ - y₂ + y₃ - y₄)))/(1 - ν₀^2)) + 0.125*(1 - η)*E₀*(1 - ν₀)*(-0.125*((1 - η^2)*(-0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ - x₂ + x₃ + x₄)) + 2*ξ*(1 - η)*(-0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ + x₂ + x₃ - x₄))) + 0.125*((1 - ξ^2)*(0.25*(-y₁ + y₂ + y₃ - y₄) - 0.25*(y₁ - y₂ + y₃ - y₄)) + 2*η*(1 - ξ)*(0.25*(-y₁ - y₂ + y₃ + y₄) - 0.25*(y₁ - y₂ + y₃ - y₄))))/(1 - ν₀^2)
            K[3,6]=K[6,3]=-0.125*(-(1 - ξ^2)*(-0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ + x₂ + x₃ - x₄)) + 2*(0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ - x₂ + x₃ + x₄))*η*(1 + ξ))*(-0.125*E₀*((1 - ξ^2)*(-0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ + x₂ + x₃ - x₄)) + 2*η*(1 - ξ)*(-0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ - x₂ + x₃ + x₄)))/(1 - ν₀^2) + 0.125*E₀*ν₀*((1 - η^2)*(0.25*(-y₁ - y₂ + y₃ + y₄) - 0.25*(y₁ - y₂ + y₃ - y₄)) + 2*ξ*(1 - η)*(0.25*(-y₁ + y₂ + y₃ - y₄) - 0.25*(y₁ - y₂ + y₃ - y₄)))/(1 - ν₀^2)) + 0.125*(-(0.25*(-y₁ - y₂ + y₃ + y₄) + 0.25*(y₁ - y₂ + y₃ - y₄))*(1 - η^2) - 2*ξ*(1 - η)*(0.25*(-y₁ + y₂ + y₃ - y₄) - 0.25*(y₁ - y₂ + y₃ - y₄)))*(0.125*E₀*((1 - η^2)*(0.25*(-y₁ - y₂ + y₃ + y₄) - 0.25*(y₁ - y₂ + y₃ - y₄)) + 2*ξ*(1 - η)*(0.25*(-y₁ + y₂ + y₃ - y₄) - 0.25*(y₁ - y₂ + y₃ - y₄)))/(1 - ν₀^2) - 0.125*E₀*ν₀*((1 - ξ^2)*(-0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ + x₂ + x₃ - x₄)) + 2*η*(1 - ξ)*(-0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ - x₂ + x₃ + x₄)))/(1 - ν₀^2)) + (1/2)*E₀*(1 - ν₀)*(-0.125*((1 - η^2)*(-0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ - x₂ + x₃ + x₄)) + 2*ξ*(1 - η)*(-0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ + x₂ + x₃ - x₄))) + 0.125*((1 - ξ^2)*(0.25*(-y₁ + y₂ + y₃ - y₄) - 0.25*(y₁ - y₂ + y₃ - y₄)) + 2*η*(1 - ξ)*(0.25*(-y₁ - y₂ + y₃ + y₄) - 0.25*(y₁ - y₂ + y₃ - y₄))))*(0.125*(-(1 - ξ^2)*(0.25*(-y₁ + y₂ + y₃ - y₄) - 0.25*(y₁ - y₂ + y₃ - y₄)) + 2*(0.25*(-y₁ - y₂ + y₃ + y₄) + 0.25*(y₁ - y₂ + y₃ - y₄))*η*(1 + ξ)) - 0.125*(-(0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ - x₂ + x₃ + x₄))*(1 - η^2) - 2*ξ*(1 - η)*(-0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ + x₂ + x₃ - x₄))))/(1 - ν₀^2)
            K[3,7]=K[7,3]=0.25*(1 + η)*(0.125*E₀*((1 - η^2)*(0.25*(-y₁ - y₂ + y₃ + y₄) - 0.25*(y₁ - y₂ + y₃ - y₄)) + 2*ξ*(1 - η)*(0.25*(-y₁ + y₂ + y₃ - y₄) - 0.25*(y₁ - y₂ + y₃ - y₄)))/(1 - ν₀^2) - 0.125*E₀*ν₀*((1 - ξ^2)*(-0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ + x₂ + x₃ - x₄)) + 2*η*(1 - ξ)*(-0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ - x₂ + x₃ + x₄)))/(1 - ν₀^2)) + 0.125*(1 + ξ)*E₀*(1 - ν₀)*(-0.125*((1 - η^2)*(-0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ - x₂ + x₃ + x₄)) + 2*ξ*(1 - η)*(-0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ + x₂ + x₃ - x₄))) + 0.125*((1 - ξ^2)*(0.25*(-y₁ + y₂ + y₃ - y₄) - 0.25*(y₁ - y₂ + y₃ - y₄)) + 2*η*(1 - ξ)*(0.25*(-y₁ - y₂ + y₃ + y₄) - 0.25*(y₁ - y₂ + y₃ - y₄))))/(1 - ν₀^2)
            K[3,8]=K[8,3]=0.25*(1 + ξ)*(-0.125*E₀*((1 - ξ^2)*(-0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ + x₂ + x₃ - x₄)) + 2*η*(1 - ξ)*(-0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ - x₂ + x₃ + x₄)))/(1 - ν₀^2) + 0.125*E₀*ν₀*((1 - η^2)*(0.25*(-y₁ - y₂ + y₃ + y₄) - 0.25*(y₁ - y₂ + y₃ - y₄)) + 2*ξ*(1 - η)*(0.25*(-y₁ + y₂ + y₃ - y₄) - 0.25*(y₁ - y₂ + y₃ - y₄)))/(1 - ν₀^2)) + 0.125*(1 + η)*E₀*(1 - ν₀)*(-0.125*((1 - η^2)*(-0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ - x₂ + x₃ + x₄)) + 2*ξ*(1 - η)*(-0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ + x₂ + x₃ - x₄))) + 0.125*((1 - ξ^2)*(0.25*(-y₁ + y₂ + y₃ - y₄) - 0.25*(y₁ - y₂ + y₃ - y₄)) + 2*η*(1 - ξ)*(0.25*(-y₁ - y₂ + y₃ + y₄) - 0.25*(y₁ - y₂ + y₃ - y₄))))/(1 - ν₀^2)
            K[3,9]=K[9,3]=0.125*((0.25*(-y₁ - y₂ + y₃ + y₄) + 0.25*(y₁ - y₂ + y₃ - y₄))*(1 - η^2) - 2*(0.25*(-y₁ + y₂ + y₃ - y₄) + 0.25*(y₁ - y₂ + y₃ - y₄))*ξ*(1 + η))*(0.125*E₀*((1 - η^2)*(0.25*(-y₁ - y₂ + y₃ + y₄) - 0.25*(y₁ - y₂ + y₃ - y₄)) + 2*ξ*(1 - η)*(0.25*(-y₁ + y₂ + y₃ - y₄) - 0.25*(y₁ - y₂ + y₃ - y₄)))/(1 - ν₀^2) - 0.125*E₀*ν₀*((1 - ξ^2)*(-0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ + x₂ + x₃ - x₄)) + 2*η*(1 - ξ)*(-0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ - x₂ + x₃ + x₄)))/(1 - ν₀^2)) - 0.125*(-0.125*E₀*((1 - ξ^2)*(-0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ + x₂ + x₃ - x₄)) + 2*η*(1 - ξ)*(-0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ - x₂ + x₃ + x₄)))/(1 - ν₀^2) + 0.125*E₀*ν₀*((1 - η^2)*(0.25*(-y₁ - y₂ + y₃ + y₄) - 0.25*(y₁ - y₂ + y₃ - y₄)) + 2*ξ*(1 - η)*(0.25*(-y₁ + y₂ + y₃ - y₄) - 0.25*(y₁ - y₂ + y₃ - y₄)))/(1 - ν₀^2))*((0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ + x₂ + x₃ - x₄))*(1 - ξ^2) - 2*(0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ - x₂ + x₃ + x₄))*η*(1 + ξ)) + (1/2)*E₀*(1 - ν₀)*(-0.125*((1 - η^2)*(-0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ - x₂ + x₃ + x₄)) + 2*ξ*(1 - η)*(-0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ + x₂ + x₃ - x₄))) + 0.125*((1 - ξ^2)*(0.25*(-y₁ + y₂ + y₃ - y₄) - 0.25*(y₁ - y₂ + y₃ - y₄)) + 2*η*(1 - ξ)*(0.25*(-y₁ - y₂ + y₃ + y₄) - 0.25*(y₁ - y₂ + y₃ - y₄))))*(-0.125*((0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ - x₂ + x₃ + x₄))*(1 - η^2) - 2*(0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ + x₂ + x₃ - x₄))*ξ*(1 + η)) + 0.125*((0.25*(-y₁ + y₂ + y₃ - y₄) + 0.25*(y₁ - y₂ + y₃ - y₄))*(1 - ξ^2) - 2*(0.25*(-y₁ - y₂ + y₃ + y₄) + 0.25*(y₁ - y₂ + y₃ - y₄))*η*(1 + ξ)))/(1 - ν₀^2)
            K[3,10]=K[10,3]=-0.25*(1 + η)*(0.125*E₀*((1 - η^2)*(0.25*(-y₁ - y₂ + y₃ + y₄) - 0.25*(y₁ - y₂ + y₃ - y₄)) + 2*ξ*(1 - η)*(0.25*(-y₁ + y₂ + y₃ - y₄) - 0.25*(y₁ - y₂ + y₃ - y₄)))/(1 - ν₀^2) - 0.125*E₀*ν₀*((1 - ξ^2)*(-0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ + x₂ + x₃ - x₄)) + 2*η*(1 - ξ)*(-0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ - x₂ + x₃ + x₄)))/(1 - ν₀^2)) + 0.125*(1 - ξ)*E₀*(1 - ν₀)*(-0.125*((1 - η^2)*(-0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ - x₂ + x₃ + x₄)) + 2*ξ*(1 - η)*(-0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ + x₂ + x₃ - x₄))) + 0.125*((1 - ξ^2)*(0.25*(-y₁ + y₂ + y₃ - y₄) - 0.25*(y₁ - y₂ + y₃ - y₄)) + 2*η*(1 - ξ)*(0.25*(-y₁ - y₂ + y₃ + y₄) - 0.25*(y₁ - y₂ + y₃ - y₄))))/(1 - ν₀^2)
            K[3,11]=K[11,3]=0.25*(1 - ξ)*(-0.125*E₀*((1 - ξ^2)*(-0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ + x₂ + x₃ - x₄)) + 2*η*(1 - ξ)*(-0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ - x₂ + x₃ + x₄)))/(1 - ν₀^2) + 0.125*E₀*ν₀*((1 - η^2)*(0.25*(-y₁ - y₂ + y₃ + y₄) - 0.25*(y₁ - y₂ + y₃ - y₄)) + 2*ξ*(1 - η)*(0.25*(-y₁ + y₂ + y₃ - y₄) - 0.25*(y₁ - y₂ + y₃ - y₄)))/(1 - ν₀^2)) - 0.125*(1 + η)*E₀*(1 - ν₀)*(-0.125*((1 - η^2)*(-0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ - x₂ + x₃ + x₄)) + 2*ξ*(1 - η)*(-0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ + x₂ + x₃ - x₄))) + 0.125*((1 - ξ^2)*(0.25*(-y₁ + y₂ + y₃ - y₄) - 0.25*(y₁ - y₂ + y₃ - y₄)) + 2*η*(1 - ξ)*(0.25*(-y₁ - y₂ + y₃ + y₄) - 0.25*(y₁ - y₂ + y₃ - y₄))))/(1 - ν₀^2)
            K[3,12]=K[12,3]=0.125*(-(1 - η^2)*(0.25*(-y₁ - y₂ + y₃ + y₄) - 0.25*(y₁ - y₂ + y₃ - y₄)) + 2*(0.25*(-y₁ + y₂ + y₃ - y₄) + 0.25*(y₁ - y₂ + y₃ - y₄))*ξ*(1 + η))*(0.125*E₀*((1 - η^2)*(0.25*(-y₁ - y₂ + y₃ + y₄) - 0.25*(y₁ - y₂ + y₃ - y₄)) + 2*ξ*(1 - η)*(0.25*(-y₁ + y₂ + y₃ - y₄) - 0.25*(y₁ - y₂ + y₃ - y₄)))/(1 - ν₀^2) - 0.125*E₀*ν₀*((1 - ξ^2)*(-0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ + x₂ + x₃ - x₄)) + 2*η*(1 - ξ)*(-0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ - x₂ + x₃ + x₄)))/(1 - ν₀^2)) - 0.125*(-(0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ + x₂ + x₃ - x₄))*(1 - ξ^2) - 2*η*(1 - ξ)*(-0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ - x₂ + x₃ + x₄)))*(-0.125*E₀*((1 - ξ^2)*(-0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ + x₂ + x₃ - x₄)) + 2*η*(1 - ξ)*(-0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ - x₂ + x₃ + x₄)))/(1 - ν₀^2) + 0.125*E₀*ν₀*((1 - η^2)*(0.25*(-y₁ - y₂ + y₃ + y₄) - 0.25*(y₁ - y₂ + y₃ - y₄)) + 2*ξ*(1 - η)*(0.25*(-y₁ + y₂ + y₃ - y₄) - 0.25*(y₁ - y₂ + y₃ - y₄)))/(1 - ν₀^2)) + (1/2)*E₀*(1 - ν₀)*(-0.125*((1 - η^2)*(-0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ - x₂ + x₃ + x₄)) + 2*ξ*(1 - η)*(-0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ + x₂ + x₃ - x₄))) + 0.125*((1 - ξ^2)*(0.25*(-y₁ + y₂ + y₃ - y₄) - 0.25*(y₁ - y₂ + y₃ - y₄)) + 2*η*(1 - ξ)*(0.25*(-y₁ - y₂ + y₃ + y₄) - 0.25*(y₁ - y₂ + y₃ - y₄))))*(-0.125*(-(1 - η^2)*(-0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ - x₂ + x₃ + x₄)) + 2*(0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ + x₂ + x₃ - x₄))*ξ*(1 + η)) + 0.125*(-(0.25*(-y₁ + y₂ + y₃ - y₄) + 0.25*(y₁ - y₂ + y₃ - y₄))*(1 - ξ^2) - 2*η*(1 - ξ)*(0.25*(-y₁ - y₂ + y₃ + y₄) - 0.25*(y₁ - y₂ + y₃ - y₄))))/(1 - ν₀^2)
            K[4,4]=0.0625*(1 - η)^2*E₀/(1 - ν₀^2) + 0.03125*(1 + ξ)^2*E₀*(1 - ν₀)/(1 - ν₀^2)
            K[4,5]=K[5,4]=-0.0625*(1 + ξ)*(1 - η)*E₀*ν₀/(1 - ν₀^2) - 0.03125*(1 + ξ)*(1 - η)*E₀*(1 - ν₀)/(1 - ν₀^2)
            K[4,6]=K[6,4]=0.03125*(1 - η)*E₀*(-(0.25*(-y₁ - y₂ + y₃ + y₄) + 0.25*(y₁ - y₂ + y₃ - y₄))*(1 - η^2) - 2*ξ*(1 - η)*(0.25*(-y₁ + y₂ + y₃ - y₄) - 0.25*(y₁ - y₂ + y₃ - y₄)))/(1 - ν₀^2) - 0.03125*(1 - η)*E₀*ν₀*(-(1 - ξ^2)*(-0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ + x₂ + x₃ - x₄)) + 2*(0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ - x₂ + x₃ + x₄))*η*(1 + ξ))/(1 - ν₀^2) - 0.125*(1 + ξ)*E₀*(1 - ν₀)*(0.125*(-(1 - ξ^2)*(0.25*(-y₁ + y₂ + y₃ - y₄) - 0.25*(y₁ - y₂ + y₃ - y₄)) + 2*(0.25*(-y₁ - y₂ + y₃ + y₄) + 0.25*(y₁ - y₂ + y₃ - y₄))*η*(1 + ξ)) - 0.125*(-(0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ - x₂ + x₃ + x₄))*(1 - η^2) - 2*ξ*(1 - η)*(-0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ + x₂ + x₃ - x₄))))/(1 - ν₀^2)
            K[4,7]=K[7,4]=0.0625*(1 - η)*(1 + η)*E₀/(1 - ν₀^2) - 0.03125*(1 + ξ)^2*E₀*(1 - ν₀)/(1 - ν₀^2)
            K[4,8]=K[8,4]=0.0625*(1 + ξ)*(1 - η)*E₀*ν₀/(1 - ν₀^2) - 0.03125*(1 + ξ)*(1 + η)*E₀*(1 - ν₀)/(1 - ν₀^2)
            K[4,9]=K[9,4]=0.03125*(1 - η)*E₀*((0.25*(-y₁ - y₂ + y₃ + y₄) + 0.25*(y₁ - y₂ + y₃ - y₄))*(1 - η^2) - 2*(0.25*(-y₁ + y₂ + y₃ - y₄) + 0.25*(y₁ - y₂ + y₃ - y₄))*ξ*(1 + η))/(1 - ν₀^2) - 0.03125*(1 - η)*E₀*ν₀*((0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ + x₂ + x₃ - x₄))*(1 - ξ^2) - 2*(0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ - x₂ + x₃ + x₄))*η*(1 + ξ))/(1 - ν₀^2) - 0.125*(1 + ξ)*E₀*(1 - ν₀)*(-0.125*((0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ - x₂ + x₃ + x₄))*(1 - η^2) - 2*(0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ + x₂ + x₃ - x₄))*ξ*(1 + η)) + 0.125*((0.25*(-y₁ + y₂ + y₃ - y₄) + 0.25*(y₁ - y₂ + y₃ - y₄))*(1 - ξ^2) - 2*(0.25*(-y₁ - y₂ + y₃ + y₄) + 0.25*(y₁ - y₂ + y₃ - y₄))*η*(1 + ξ)))/(1 - ν₀^2)
            K[4,10]=K[10,4]=-0.0625*(1 - η)*(1 + η)*E₀/(1 - ν₀^2) - 0.03125*(1 - ξ)*(1 + ξ)*E₀*(1 - ν₀)/(1 - ν₀^2)
            K[4,11]=K[11,4]=0.0625*(1 - ξ)*(1 - η)*E₀*ν₀/(1 - ν₀^2) + 0.03125*(1 + ξ)*(1 + η)*E₀*(1 - ν₀)/(1 - ν₀^2)
            K[4,12]=K[12,4]=0.03125*(1 - η)*E₀*(-(1 - η^2)*(0.25*(-y₁ - y₂ + y₃ + y₄) - 0.25*(y₁ - y₂ + y₃ - y₄)) + 2*(0.25*(-y₁ + y₂ + y₃ - y₄) + 0.25*(y₁ - y₂ + y₃ - y₄))*ξ*(1 + η))/(1 - ν₀^2) - 0.03125*(1 - η)*E₀*ν₀*(-(0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ + x₂ + x₃ - x₄))*(1 - ξ^2) - 2*η*(1 - ξ)*(-0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ - x₂ + x₃ + x₄)))/(1 - ν₀^2) - 0.125*(1 + ξ)*E₀*(1 - ν₀)*(-0.125*(-(1 - η^2)*(-0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ - x₂ + x₃ + x₄)) + 2*(0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ + x₂ + x₃ - x₄))*ξ*(1 + η)) + 0.125*(-(0.25*(-y₁ + y₂ + y₃ - y₄) + 0.25*(y₁ - y₂ + y₃ - y₄))*(1 - ξ^2) - 2*η*(1 - ξ)*(0.25*(-y₁ - y₂ + y₃ + y₄) - 0.25*(y₁ - y₂ + y₃ - y₄))))/(1 - ν₀^2)
            K[5,5]=0.0625*(1 + ξ)^2*E₀/(1 - ν₀^2) + 0.03125*(1 - η)^2*E₀*(1 - ν₀)/(1 - ν₀^2)
            K[5,6]=K[6,5]=0.03125*(1 + ξ)*E₀*(-(1 - ξ^2)*(-0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ + x₂ + x₃ - x₄)) + 2*(0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ - x₂ + x₃ + x₄))*η*(1 + ξ))/(1 - ν₀^2) + 0.125*(1 - η)*E₀*(1 - ν₀)*(0.125*(-(1 - ξ^2)*(0.25*(-y₁ + y₂ + y₃ - y₄) - 0.25*(y₁ - y₂ + y₃ - y₄)) + 2*(0.25*(-y₁ - y₂ + y₃ + y₄) + 0.25*(y₁ - y₂ + y₃ - y₄))*η*(1 + ξ)) - 0.125*(-(0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ - x₂ + x₃ + x₄))*(1 - η^2) - 2*ξ*(1 - η)*(-0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ + x₂ + x₃ - x₄))))/(1 - ν₀^2) - 0.03125*(1 + ξ)*E₀*ν₀*(-(0.25*(-y₁ - y₂ + y₃ + y₄) + 0.25*(y₁ - y₂ + y₃ - y₄))*(1 - η^2) - 2*ξ*(1 - η)*(0.25*(-y₁ + y₂ + y₃ - y₄) - 0.25*(y₁ - y₂ + y₃ - y₄)))/(1 - ν₀^2)
            K[5,7]=K[7,5]=0.03125*(1 + ξ)*(1 - η)*E₀*(1 - ν₀)/(1 - ν₀^2) - 0.0625*(1 + ξ)*(1 + η)*E₀*ν₀/(1 - ν₀^2)
            K[5,8]=K[8,5]=-0.0625*(1 + ξ)^2*E₀/(1 - ν₀^2) + 0.03125*(1 - η)*(1 + η)*E₀*(1 - ν₀)/(1 - ν₀^2)
            K[5,9]=K[9,5]=0.03125*(1 + ξ)*E₀*((0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ + x₂ + x₃ - x₄))*(1 - ξ^2) - 2*(0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ - x₂ + x₃ + x₄))*η*(1 + ξ))/(1 - ν₀^2) + 0.125*(1 - η)*E₀*(1 - ν₀)*(-0.125*((0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ - x₂ + x₃ + x₄))*(1 - η^2) - 2*(0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ + x₂ + x₃ - x₄))*ξ*(1 + η)) + 0.125*((0.25*(-y₁ + y₂ + y₃ - y₄) + 0.25*(y₁ - y₂ + y₃ - y₄))*(1 - ξ^2) - 2*(0.25*(-y₁ - y₂ + y₃ + y₄) + 0.25*(y₁ - y₂ + y₃ - y₄))*η*(1 + ξ)))/(1 - ν₀^2) - 0.03125*(1 + ξ)*E₀*ν₀*((0.25*(-y₁ - y₂ + y₃ + y₄) + 0.25*(y₁ - y₂ + y₃ - y₄))*(1 - η^2) - 2*(0.25*(-y₁ + y₂ + y₃ - y₄) + 0.25*(y₁ - y₂ + y₃ - y₄))*ξ*(1 + η))/(1 - ν₀^2)
            K[5,10]=K[10,5]=0.03125*(1 - ξ)*(1 - η)*E₀*(1 - ν₀)/(1 - ν₀^2) + 0.0625*(1 + ξ)*(1 + η)*E₀*ν₀/(1 - ν₀^2)
            K[5,11]=K[11,5]=-0.0625*(1 - ξ)*(1 + ξ)*E₀/(1 - ν₀^2) - 0.03125*(1 - η)*(1 + η)*E₀*(1 - ν₀)/(1 - ν₀^2)
            K[5,12]=K[12,5]=0.03125*(1 + ξ)*E₀*(-(0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ + x₂ + x₃ - x₄))*(1 - ξ^2) - 2*η*(1 - ξ)*(-0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ - x₂ + x₃ + x₄)))/(1 - ν₀^2) + 0.125*(1 - η)*E₀*(1 - ν₀)*(-0.125*(-(1 - η^2)*(-0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ - x₂ + x₃ + x₄)) + 2*(0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ + x₂ + x₃ - x₄))*ξ*(1 + η)) + 0.125*(-(0.25*(-y₁ + y₂ + y₃ - y₄) + 0.25*(y₁ - y₂ + y₃ - y₄))*(1 - ξ^2) - 2*η*(1 - ξ)*(0.25*(-y₁ - y₂ + y₃ + y₄) - 0.25*(y₁ - y₂ + y₃ - y₄))))/(1 - ν₀^2) - 0.03125*(1 + ξ)*E₀*ν₀*(-(1 - η^2)*(0.25*(-y₁ - y₂ + y₃ + y₄) - 0.25*(y₁ - y₂ + y₃ - y₄)) + 2*(0.25*(-y₁ + y₂ + y₃ - y₄) + 0.25*(y₁ - y₂ + y₃ - y₄))*ξ*(1 + η))/(1 - ν₀^2)
            K[6,6]=0.125*(-(0.25*(-y₁ - y₂ + y₃ + y₄) + 0.25*(y₁ - y₂ + y₃ - y₄))*(1 - η^2) - 2*ξ*(1 - η)*(0.25*(-y₁ + y₂ + y₃ - y₄) - 0.25*(y₁ - y₂ + y₃ - y₄)))*(0.125*E₀*(-(0.25*(-y₁ - y₂ + y₃ + y₄) + 0.25*(y₁ - y₂ + y₃ - y₄))*(1 - η^2) - 2*ξ*(1 - η)*(0.25*(-y₁ + y₂ + y₃ - y₄) - 0.25*(y₁ - y₂ + y₃ - y₄)))/(1 - ν₀^2) - 0.125*E₀*ν₀*(-(1 - ξ^2)*(-0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ + x₂ + x₃ - x₄)) + 2*(0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ - x₂ + x₃ + x₄))*η*(1 + ξ))/(1 - ν₀^2)) - 0.125*(-0.125*E₀*(-(1 - ξ^2)*(-0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ + x₂ + x₃ - x₄)) + 2*(0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ - x₂ + x₃ + x₄))*η*(1 + ξ))/(1 - ν₀^2) + 0.125*E₀*ν₀*(-(0.25*(-y₁ - y₂ + y₃ + y₄) + 0.25*(y₁ - y₂ + y₃ - y₄))*(1 - η^2) - 2*ξ*(1 - η)*(0.25*(-y₁ + y₂ + y₃ - y₄) - 0.25*(y₁ - y₂ + y₃ - y₄)))/(1 - ν₀^2))*(-(1 - ξ^2)*(-0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ + x₂ + x₃ - x₄)) + 2*(0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ - x₂ + x₃ + x₄))*η*(1 + ξ)) + (1/2)*E₀*(1 - ν₀)*(0.125*(-(1 - ξ^2)*(0.25*(-y₁ + y₂ + y₃ - y₄) - 0.25*(y₁ - y₂ + y₃ - y₄)) + 2*(0.25*(-y₁ - y₂ + y₃ + y₄) + 0.25*(y₁ - y₂ + y₃ - y₄))*η*(1 + ξ)) - 0.125*(-(0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ - x₂ + x₃ + x₄))*(1 - η^2) - 2*ξ*(1 - η)*(-0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ + x₂ + x₃ - x₄))))^2/(1 - ν₀^2)
            K[6,7]=K[7,6]=0.25*(1 + η)*(0.125*E₀*(-(0.25*(-y₁ - y₂ + y₃ + y₄) + 0.25*(y₁ - y₂ + y₃ - y₄))*(1 - η^2) - 2*ξ*(1 - η)*(0.25*(-y₁ + y₂ + y₃ - y₄) - 0.25*(y₁ - y₂ + y₃ - y₄)))/(1 - ν₀^2) - 0.125*E₀*ν₀*(-(1 - ξ^2)*(-0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ + x₂ + x₃ - x₄)) + 2*(0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ - x₂ + x₃ + x₄))*η*(1 + ξ))/(1 - ν₀^2)) + 0.125*(1 + ξ)*E₀*(1 - ν₀)*(0.125*(-(1 - ξ^2)*(0.25*(-y₁ + y₂ + y₃ - y₄) - 0.25*(y₁ - y₂ + y₃ - y₄)) + 2*(0.25*(-y₁ - y₂ + y₃ + y₄) + 0.25*(y₁ - y₂ + y₃ - y₄))*η*(1 + ξ)) - 0.125*(-(0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ - x₂ + x₃ + x₄))*(1 - η^2) - 2*ξ*(1 - η)*(-0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ + x₂ + x₃ - x₄))))/(1 - ν₀^2)
            K[6,8]=K[8,6]=0.25*(1 + ξ)*(-0.125*E₀*(-(1 - ξ^2)*(-0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ + x₂ + x₃ - x₄)) + 2*(0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ - x₂ + x₃ + x₄))*η*(1 + ξ))/(1 - ν₀^2) + 0.125*E₀*ν₀*(-(0.25*(-y₁ - y₂ + y₃ + y₄) + 0.25*(y₁ - y₂ + y₃ - y₄))*(1 - η^2) - 2*ξ*(1 - η)*(0.25*(-y₁ + y₂ + y₃ - y₄) - 0.25*(y₁ - y₂ + y₃ - y₄)))/(1 - ν₀^2)) + 0.125*(1 + η)*E₀*(1 - ν₀)*(0.125*(-(1 - ξ^2)*(0.25*(-y₁ + y₂ + y₃ - y₄) - 0.25*(y₁ - y₂ + y₃ - y₄)) + 2*(0.25*(-y₁ - y₂ + y₃ + y₄) + 0.25*(y₁ - y₂ + y₃ - y₄))*η*(1 + ξ)) - 0.125*(-(0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ - x₂ + x₃ + x₄))*(1 - η^2) - 2*ξ*(1 - η)*(-0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ + x₂ + x₃ - x₄))))/(1 - ν₀^2)
            K[6,9]=K[9,6]=0.125*(0.125*E₀*(-(0.25*(-y₁ - y₂ + y₃ + y₄) + 0.25*(y₁ - y₂ + y₃ - y₄))*(1 - η^2) - 2*ξ*(1 - η)*(0.25*(-y₁ + y₂ + y₃ - y₄) - 0.25*(y₁ - y₂ + y₃ - y₄)))/(1 - ν₀^2) - 0.125*E₀*ν₀*(-(1 - ξ^2)*(-0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ + x₂ + x₃ - x₄)) + 2*(0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ - x₂ + x₃ + x₄))*η*(1 + ξ))/(1 - ν₀^2))*((0.25*(-y₁ - y₂ + y₃ + y₄) + 0.25*(y₁ - y₂ + y₃ - y₄))*(1 - η^2) - 2*(0.25*(-y₁ + y₂ + y₃ - y₄) + 0.25*(y₁ - y₂ + y₃ - y₄))*ξ*(1 + η)) - 0.125*(-0.125*E₀*(-(1 - ξ^2)*(-0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ + x₂ + x₃ - x₄)) + 2*(0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ - x₂ + x₃ + x₄))*η*(1 + ξ))/(1 - ν₀^2) + 0.125*E₀*ν₀*(-(0.25*(-y₁ - y₂ + y₃ + y₄) + 0.25*(y₁ - y₂ + y₃ - y₄))*(1 - η^2) - 2*ξ*(1 - η)*(0.25*(-y₁ + y₂ + y₃ - y₄) - 0.25*(y₁ - y₂ + y₃ - y₄)))/(1 - ν₀^2))*((0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ + x₂ + x₃ - x₄))*(1 - ξ^2) - 2*(0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ - x₂ + x₃ + x₄))*η*(1 + ξ)) + (1/2)*E₀*(1 - ν₀)*(-0.125*((0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ - x₂ + x₃ + x₄))*(1 - η^2) - 2*(0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ + x₂ + x₃ - x₄))*ξ*(1 + η)) + 0.125*((0.25*(-y₁ + y₂ + y₃ - y₄) + 0.25*(y₁ - y₂ + y₃ - y₄))*(1 - ξ^2) - 2*(0.25*(-y₁ - y₂ + y₃ + y₄) + 0.25*(y₁ - y₂ + y₃ - y₄))*η*(1 + ξ)))*(0.125*(-(1 - ξ^2)*(0.25*(-y₁ + y₂ + y₃ - y₄) - 0.25*(y₁ - y₂ + y₃ - y₄)) + 2*(0.25*(-y₁ - y₂ + y₃ + y₄) + 0.25*(y₁ - y₂ + y₃ - y₄))*η*(1 + ξ)) - 0.125*(-(0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ - x₂ + x₃ + x₄))*(1 - η^2) - 2*ξ*(1 - η)*(-0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ + x₂ + x₃ - x₄))))/(1 - ν₀^2)
            K[6,10]=K[10,6]=-0.25*(1 + η)*(0.125*E₀*(-(0.25*(-y₁ - y₂ + y₃ + y₄) + 0.25*(y₁ - y₂ + y₃ - y₄))*(1 - η^2) - 2*ξ*(1 - η)*(0.25*(-y₁ + y₂ + y₃ - y₄) - 0.25*(y₁ - y₂ + y₃ - y₄)))/(1 - ν₀^2) - 0.125*E₀*ν₀*(-(1 - ξ^2)*(-0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ + x₂ + x₃ - x₄)) + 2*(0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ - x₂ + x₃ + x₄))*η*(1 + ξ))/(1 - ν₀^2)) + 0.125*(1 - ξ)*E₀*(1 - ν₀)*(0.125*(-(1 - ξ^2)*(0.25*(-y₁ + y₂ + y₃ - y₄) - 0.25*(y₁ - y₂ + y₃ - y₄)) + 2*(0.25*(-y₁ - y₂ + y₃ + y₄) + 0.25*(y₁ - y₂ + y₃ - y₄))*η*(1 + ξ)) - 0.125*(-(0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ - x₂ + x₃ + x₄))*(1 - η^2) - 2*ξ*(1 - η)*(-0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ + x₂ + x₃ - x₄))))/(1 - ν₀^2)
            K[6,11]=K[11,6]=0.25*(1 - ξ)*(-0.125*E₀*(-(1 - ξ^2)*(-0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ + x₂ + x₃ - x₄)) + 2*(0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ - x₂ + x₃ + x₄))*η*(1 + ξ))/(1 - ν₀^2) + 0.125*E₀*ν₀*(-(0.25*(-y₁ - y₂ + y₃ + y₄) + 0.25*(y₁ - y₂ + y₃ - y₄))*(1 - η^2) - 2*ξ*(1 - η)*(0.25*(-y₁ + y₂ + y₃ - y₄) - 0.25*(y₁ - y₂ + y₃ - y₄)))/(1 - ν₀^2)) - 0.125*(1 + η)*E₀*(1 - ν₀)*(0.125*(-(1 - ξ^2)*(0.25*(-y₁ + y₂ + y₃ - y₄) - 0.25*(y₁ - y₂ + y₃ - y₄)) + 2*(0.25*(-y₁ - y₂ + y₃ + y₄) + 0.25*(y₁ - y₂ + y₃ - y₄))*η*(1 + ξ)) - 0.125*(-(0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ - x₂ + x₃ + x₄))*(1 - η^2) - 2*ξ*(1 - η)*(-0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ + x₂ + x₃ - x₄))))/(1 - ν₀^2)
            K[6,12]=K[12,6]=0.125*(0.125*E₀*(-(0.25*(-y₁ - y₂ + y₃ + y₄) + 0.25*(y₁ - y₂ + y₃ - y₄))*(1 - η^2) - 2*ξ*(1 - η)*(0.25*(-y₁ + y₂ + y₃ - y₄) - 0.25*(y₁ - y₂ + y₃ - y₄)))/(1 - ν₀^2) - 0.125*E₀*ν₀*(-(1 - ξ^2)*(-0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ + x₂ + x₃ - x₄)) + 2*(0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ - x₂ + x₃ + x₄))*η*(1 + ξ))/(1 - ν₀^2))*(-(1 - η^2)*(0.25*(-y₁ - y₂ + y₃ + y₄) - 0.25*(y₁ - y₂ + y₃ - y₄)) + 2*(0.25*(-y₁ + y₂ + y₃ - y₄) + 0.25*(y₁ - y₂ + y₃ - y₄))*ξ*(1 + η)) - 0.125*(-0.125*E₀*(-(1 - ξ^2)*(-0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ + x₂ + x₃ - x₄)) + 2*(0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ - x₂ + x₃ + x₄))*η*(1 + ξ))/(1 - ν₀^2) + 0.125*E₀*ν₀*(-(0.25*(-y₁ - y₂ + y₃ + y₄) + 0.25*(y₁ - y₂ + y₃ - y₄))*(1 - η^2) - 2*ξ*(1 - η)*(0.25*(-y₁ + y₂ + y₃ - y₄) - 0.25*(y₁ - y₂ + y₃ - y₄)))/(1 - ν₀^2))*(-(0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ + x₂ + x₃ - x₄))*(1 - ξ^2) - 2*η*(1 - ξ)*(-0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ - x₂ + x₃ + x₄))) + (1/2)*E₀*(1 - ν₀)*(-0.125*(-(1 - η^2)*(-0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ - x₂ + x₃ + x₄)) + 2*(0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ + x₂ + x₃ - x₄))*ξ*(1 + η)) + 0.125*(-(0.25*(-y₁ + y₂ + y₃ - y₄) + 0.25*(y₁ - y₂ + y₃ - y₄))*(1 - ξ^2) - 2*η*(1 - ξ)*(0.25*(-y₁ - y₂ + y₃ + y₄) - 0.25*(y₁ - y₂ + y₃ - y₄))))*(0.125*(-(1 - ξ^2)*(0.25*(-y₁ + y₂ + y₃ - y₄) - 0.25*(y₁ - y₂ + y₃ - y₄)) + 2*(0.25*(-y₁ - y₂ + y₃ + y₄) + 0.25*(y₁ - y₂ + y₃ - y₄))*η*(1 + ξ)) - 0.125*(-(0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ - x₂ + x₃ + x₄))*(1 - η^2) - 2*ξ*(1 - η)*(-0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ + x₂ + x₃ - x₄))))/(1 - ν₀^2)
            K[7,7]=0.0625*(1 + η)^2*E₀/(1 - ν₀^2) + 0.03125*(1 + ξ)^2*E₀*(1 - ν₀)/(1 - ν₀^2)
            K[7,8]=K[8,7]=0.0625*(1 + ξ)*(1 + η)*E₀*ν₀/(1 - ν₀^2) + 0.03125*(1 + ξ)*(1 + η)*E₀*(1 - ν₀)/(1 - ν₀^2)
            K[7,9]=K[9,7]=0.03125*(1 + η)*E₀*((0.25*(-y₁ - y₂ + y₃ + y₄) + 0.25*(y₁ - y₂ + y₃ - y₄))*(1 - η^2) - 2*(0.25*(-y₁ + y₂ + y₃ - y₄) + 0.25*(y₁ - y₂ + y₃ - y₄))*ξ*(1 + η))/(1 - ν₀^2) - 0.03125*(1 + η)*E₀*ν₀*((0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ + x₂ + x₃ - x₄))*(1 - ξ^2) - 2*(0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ - x₂ + x₃ + x₄))*η*(1 + ξ))/(1 - ν₀^2) + 0.125*(1 + ξ)*E₀*(1 - ν₀)*(-0.125*((0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ - x₂ + x₃ + x₄))*(1 - η^2) - 2*(0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ + x₂ + x₃ - x₄))*ξ*(1 + η)) + 0.125*((0.25*(-y₁ + y₂ + y₃ - y₄) + 0.25*(y₁ - y₂ + y₃ - y₄))*(1 - ξ^2) - 2*(0.25*(-y₁ - y₂ + y₃ + y₄) + 0.25*(y₁ - y₂ + y₃ - y₄))*η*(1 + ξ)))/(1 - ν₀^2)
            K[7,10]=K[10,7]=-0.0625*(1 + η)^2*E₀/(1 - ν₀^2) + 0.03125*(1 - ξ)*(1 + ξ)*E₀*(1 - ν₀)/(1 - ν₀^2)
            K[7,11]=K[11,7]=0.0625*(1 - ξ)*(1 + η)*E₀*ν₀/(1 - ν₀^2) - 0.03125*(1 + ξ)*(1 + η)*E₀*(1 - ν₀)/(1 - ν₀^2)
            K[7,12]=K[12,7]=0.03125*(1 + η)*E₀*(-(1 - η^2)*(0.25*(-y₁ - y₂ + y₃ + y₄) - 0.25*(y₁ - y₂ + y₃ - y₄)) + 2*(0.25*(-y₁ + y₂ + y₃ - y₄) + 0.25*(y₁ - y₂ + y₃ - y₄))*ξ*(1 + η))/(1 - ν₀^2) - 0.03125*(1 + η)*E₀*ν₀*(-(0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ + x₂ + x₃ - x₄))*(1 - ξ^2) - 2*η*(1 - ξ)*(-0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ - x₂ + x₃ + x₄)))/(1 - ν₀^2) + 0.125*(1 + ξ)*E₀*(1 - ν₀)*(-0.125*(-(1 - η^2)*(-0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ - x₂ + x₃ + x₄)) + 2*(0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ + x₂ + x₃ - x₄))*ξ*(1 + η)) + 0.125*(-(0.25*(-y₁ + y₂ + y₃ - y₄) + 0.25*(y₁ - y₂ + y₃ - y₄))*(1 - ξ^2) - 2*η*(1 - ξ)*(0.25*(-y₁ - y₂ + y₃ + y₄) - 0.25*(y₁ - y₂ + y₃ - y₄))))/(1 - ν₀^2)
            K[8,8]=0.0625*(1 + ξ)^2*E₀/(1 - ν₀^2) + 0.03125*(1 + η)^2*E₀*(1 - ν₀)/(1 - ν₀^2)
            K[8,9]=K[9,8]=-0.03125*(1 + ξ)*E₀*((0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ + x₂ + x₃ - x₄))*(1 - ξ^2) - 2*(0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ - x₂ + x₃ + x₄))*η*(1 + ξ))/(1 - ν₀^2) + 0.125*(1 + η)*E₀*(1 - ν₀)*(-0.125*((0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ - x₂ + x₃ + x₄))*(1 - η^2) - 2*(0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ + x₂ + x₃ - x₄))*ξ*(1 + η)) + 0.125*((0.25*(-y₁ + y₂ + y₃ - y₄) + 0.25*(y₁ - y₂ + y₃ - y₄))*(1 - ξ^2) - 2*(0.25*(-y₁ - y₂ + y₃ + y₄) + 0.25*(y₁ - y₂ + y₃ - y₄))*η*(1 + ξ)))/(1 - ν₀^2) + 0.03125*(1 + ξ)*E₀*ν₀*((0.25*(-y₁ - y₂ + y₃ + y₄) + 0.25*(y₁ - y₂ + y₃ - y₄))*(1 - η^2) - 2*(0.25*(-y₁ + y₂ + y₃ - y₄) + 0.25*(y₁ - y₂ + y₃ - y₄))*ξ*(1 + η))/(1 - ν₀^2)
            K[8,10]=K[10,8]=0.03125*(1 - ξ)*(1 + η)*E₀*(1 - ν₀)/(1 - ν₀^2) - 0.0625*(1 + ξ)*(1 + η)*E₀*ν₀/(1 - ν₀^2)
            K[8,11]=K[11,8]=-0.03125*(1 + η)^2*E₀*(1 - ν₀)/(1 - ν₀^2) + 0.0625*(1 - ξ)*(1 + ξ)*E₀/(1 - ν₀^2)
            K[8,12]=K[12,8]=-0.03125*(1 + ξ)*E₀*(-(0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ + x₂ + x₃ - x₄))*(1 - ξ^2) - 2*η*(1 - ξ)*(-0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ - x₂ + x₃ + x₄)))/(1 - ν₀^2) + 0.125*(1 + η)*E₀*(1 - ν₀)*(-0.125*(-(1 - η^2)*(-0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ - x₂ + x₃ + x₄)) + 2*(0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ + x₂ + x₃ - x₄))*ξ*(1 + η)) + 0.125*(-(0.25*(-y₁ + y₂ + y₃ - y₄) + 0.25*(y₁ - y₂ + y₃ - y₄))*(1 - ξ^2) - 2*η*(1 - ξ)*(0.25*(-y₁ - y₂ + y₃ + y₄) - 0.25*(y₁ - y₂ + y₃ - y₄))))/(1 - ν₀^2) + 0.03125*(1 + ξ)*E₀*ν₀*(-(1 - η^2)*(0.25*(-y₁ - y₂ + y₃ + y₄) - 0.25*(y₁ - y₂ + y₃ - y₄)) + 2*(0.25*(-y₁ + y₂ + y₃ - y₄) + 0.25*(y₁ - y₂ + y₃ - y₄))*ξ*(1 + η))/(1 - ν₀^2)
            K[9,9]=-0.125*((0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ + x₂ + x₃ - x₄))*(1 - ξ^2) - 2*(0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ - x₂ + x₃ + x₄))*η*(1 + ξ))*(-0.125*E₀*((0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ + x₂ + x₃ - x₄))*(1 - ξ^2) - 2*(0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ - x₂ + x₃ + x₄))*η*(1 + ξ))/(1 - ν₀^2) + 0.125*E₀*ν₀*((0.25*(-y₁ - y₂ + y₃ + y₄) + 0.25*(y₁ - y₂ + y₃ - y₄))*(1 - η^2) - 2*(0.25*(-y₁ + y₂ + y₃ - y₄) + 0.25*(y₁ - y₂ + y₃ - y₄))*ξ*(1 + η))/(1 - ν₀^2)) + 0.125*(0.125*E₀*((0.25*(-y₁ - y₂ + y₃ + y₄) + 0.25*(y₁ - y₂ + y₃ - y₄))*(1 - η^2) - 2*(0.25*(-y₁ + y₂ + y₃ - y₄) + 0.25*(y₁ - y₂ + y₃ - y₄))*ξ*(1 + η))/(1 - ν₀^2) - 0.125*E₀*ν₀*((0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ + x₂ + x₃ - x₄))*(1 - ξ^2) - 2*(0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ - x₂ + x₃ + x₄))*η*(1 + ξ))/(1 - ν₀^2))*((0.25*(-y₁ - y₂ + y₃ + y₄) + 0.25*(y₁ - y₂ + y₃ - y₄))*(1 - η^2) - 2*(0.25*(-y₁ + y₂ + y₃ - y₄) + 0.25*(y₁ - y₂ + y₃ - y₄))*ξ*(1 + η)) + (1/2)*E₀*(1 - ν₀)*(-0.125*((0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ - x₂ + x₃ + x₄))*(1 - η^2) - 2*(0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ + x₂ + x₃ - x₄))*ξ*(1 + η)) + 0.125*((0.25*(-y₁ + y₂ + y₃ - y₄) + 0.25*(y₁ - y₂ + y₃ - y₄))*(1 - ξ^2) - 2*(0.25*(-y₁ - y₂ + y₃ + y₄) + 0.25*(y₁ - y₂ + y₃ - y₄))*η*(1 + ξ)))^2/(1 - ν₀^2)
            K[9,10]=K[10,9]=-0.25*(1 + η)*(0.125*E₀*((0.25*(-y₁ - y₂ + y₃ + y₄) + 0.25*(y₁ - y₂ + y₃ - y₄))*(1 - η^2) - 2*(0.25*(-y₁ + y₂ + y₃ - y₄) + 0.25*(y₁ - y₂ + y₃ - y₄))*ξ*(1 + η))/(1 - ν₀^2) - 0.125*E₀*ν₀*((0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ + x₂ + x₃ - x₄))*(1 - ξ^2) - 2*(0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ - x₂ + x₃ + x₄))*η*(1 + ξ))/(1 - ν₀^2)) + 0.125*(1 - ξ)*E₀*(1 - ν₀)*(-0.125*((0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ - x₂ + x₃ + x₄))*(1 - η^2) - 2*(0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ + x₂ + x₃ - x₄))*ξ*(1 + η)) + 0.125*((0.25*(-y₁ + y₂ + y₃ - y₄) + 0.25*(y₁ - y₂ + y₃ - y₄))*(1 - ξ^2) - 2*(0.25*(-y₁ - y₂ + y₃ + y₄) + 0.25*(y₁ - y₂ + y₃ - y₄))*η*(1 + ξ)))/(1 - ν₀^2)
            K[9,11]=K[11,9]=0.25*(1 - ξ)*(-0.125*E₀*((0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ + x₂ + x₃ - x₄))*(1 - ξ^2) - 2*(0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ - x₂ + x₃ + x₄))*η*(1 + ξ))/(1 - ν₀^2) + 0.125*E₀*ν₀*((0.25*(-y₁ - y₂ + y₃ + y₄) + 0.25*(y₁ - y₂ + y₃ - y₄))*(1 - η^2) - 2*(0.25*(-y₁ + y₂ + y₃ - y₄) + 0.25*(y₁ - y₂ + y₃ - y₄))*ξ*(1 + η))/(1 - ν₀^2)) - 0.125*(1 + η)*E₀*(1 - ν₀)*(-0.125*((0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ - x₂ + x₃ + x₄))*(1 - η^2) - 2*(0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ + x₂ + x₃ - x₄))*ξ*(1 + η)) + 0.125*((0.25*(-y₁ + y₂ + y₃ - y₄) + 0.25*(y₁ - y₂ + y₃ - y₄))*(1 - ξ^2) - 2*(0.25*(-y₁ - y₂ + y₃ + y₄) + 0.25*(y₁ - y₂ + y₃ - y₄))*η*(1 + ξ)))/(1 - ν₀^2)
            K[9,12]=K[12,9]=-0.125*(-(0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ + x₂ + x₃ - x₄))*(1 - ξ^2) - 2*η*(1 - ξ)*(-0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ - x₂ + x₃ + x₄)))*(-0.125*E₀*((0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ + x₂ + x₃ - x₄))*(1 - ξ^2) - 2*(0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ - x₂ + x₃ + x₄))*η*(1 + ξ))/(1 - ν₀^2) + 0.125*E₀*ν₀*((0.25*(-y₁ - y₂ + y₃ + y₄) + 0.25*(y₁ - y₂ + y₃ - y₄))*(1 - η^2) - 2*(0.25*(-y₁ + y₂ + y₃ - y₄) + 0.25*(y₁ - y₂ + y₃ - y₄))*ξ*(1 + η))/(1 - ν₀^2)) + 0.125*(0.125*E₀*((0.25*(-y₁ - y₂ + y₃ + y₄) + 0.25*(y₁ - y₂ + y₃ - y₄))*(1 - η^2) - 2*(0.25*(-y₁ + y₂ + y₃ - y₄) + 0.25*(y₁ - y₂ + y₃ - y₄))*ξ*(1 + η))/(1 - ν₀^2) - 0.125*E₀*ν₀*((0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ + x₂ + x₃ - x₄))*(1 - ξ^2) - 2*(0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ - x₂ + x₃ + x₄))*η*(1 + ξ))/(1 - ν₀^2))*(-(1 - η^2)*(0.25*(-y₁ - y₂ + y₃ + y₄) - 0.25*(y₁ - y₂ + y₃ - y₄)) + 2*(0.25*(-y₁ + y₂ + y₃ - y₄) + 0.25*(y₁ - y₂ + y₃ - y₄))*ξ*(1 + η)) + (1/2)*E₀*(1 - ν₀)*(-0.125*((0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ - x₂ + x₃ + x₄))*(1 - η^2) - 2*(0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ + x₂ + x₃ - x₄))*ξ*(1 + η)) + 0.125*((0.25*(-y₁ + y₂ + y₃ - y₄) + 0.25*(y₁ - y₂ + y₃ - y₄))*(1 - ξ^2) - 2*(0.25*(-y₁ - y₂ + y₃ + y₄) + 0.25*(y₁ - y₂ + y₃ - y₄))*η*(1 + ξ)))*(-0.125*(-(1 - η^2)*(-0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ - x₂ + x₃ + x₄)) + 2*(0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ + x₂ + x₃ - x₄))*ξ*(1 + η)) + 0.125*(-(0.25*(-y₁ + y₂ + y₃ - y₄) + 0.25*(y₁ - y₂ + y₃ - y₄))*(1 - ξ^2) - 2*η*(1 - ξ)*(0.25*(-y₁ - y₂ + y₃ + y₄) - 0.25*(y₁ - y₂ + y₃ - y₄))))/(1 - ν₀^2)
            K[10,10]=0.0625*(1 + η)^2*E₀/(1 - ν₀^2) + 0.03125*(1 - ξ)^2*E₀*(1 - ν₀)/(1 - ν₀^2)
            K[10,11]=K[11,10]=-0.0625*(1 - ξ)*(1 + η)*E₀*ν₀/(1 - ν₀^2) - 0.03125*(1 - ξ)*(1 + η)*E₀*(1 - ν₀)/(1 - ν₀^2)
            K[10,12]=K[12,10]=-0.03125*(1 + η)*E₀*(-(1 - η^2)*(0.25*(-y₁ - y₂ + y₃ + y₄) - 0.25*(y₁ - y₂ + y₃ - y₄)) + 2*(0.25*(-y₁ + y₂ + y₃ - y₄) + 0.25*(y₁ - y₂ + y₃ - y₄))*ξ*(1 + η))/(1 - ν₀^2) + 0.03125*(1 + η)*E₀*ν₀*(-(0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ + x₂ + x₃ - x₄))*(1 - ξ^2) - 2*η*(1 - ξ)*(-0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ - x₂ + x₃ + x₄)))/(1 - ν₀^2) + 0.125*(1 - ξ)*E₀*(1 - ν₀)*(-0.125*(-(1 - η^2)*(-0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ - x₂ + x₃ + x₄)) + 2*(0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ + x₂ + x₃ - x₄))*ξ*(1 + η)) + 0.125*(-(0.25*(-y₁ + y₂ + y₃ - y₄) + 0.25*(y₁ - y₂ + y₃ - y₄))*(1 - ξ^2) - 2*η*(1 - ξ)*(0.25*(-y₁ - y₂ + y₃ + y₄) - 0.25*(y₁ - y₂ + y₃ - y₄))))/(1 - ν₀^2)
            K[11,11]=0.0625*(1 - ξ)^2*E₀/(1 - ν₀^2) + 0.03125*(1 + η)^2*E₀*(1 - ν₀)/(1 - ν₀^2)
            K[11,12]=K[12,11]=-0.03125*(1 - ξ)*E₀*(-(0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ + x₂ + x₃ - x₄))*(1 - ξ^2) - 2*η*(1 - ξ)*(-0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ - x₂ + x₃ + x₄)))/(1 - ν₀^2) - 0.125*(1 + η)*E₀*(1 - ν₀)*(-0.125*(-(1 - η^2)*(-0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ - x₂ + x₃ + x₄)) + 2*(0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ + x₂ + x₃ - x₄))*ξ*(1 + η)) + 0.125*(-(0.25*(-y₁ + y₂ + y₃ - y₄) + 0.25*(y₁ - y₂ + y₃ - y₄))*(1 - ξ^2) - 2*η*(1 - ξ)*(0.25*(-y₁ - y₂ + y₃ + y₄) - 0.25*(y₁ - y₂ + y₃ - y₄))))/(1 - ν₀^2) + 0.03125*(1 - ξ)*E₀*ν₀*(-(1 - η^2)*(0.25*(-y₁ - y₂ + y₃ + y₄) - 0.25*(y₁ - y₂ + y₃ - y₄)) + 2*(0.25*(-y₁ + y₂ + y₃ - y₄) + 0.25*(y₁ - y₂ + y₃ - y₄))*ξ*(1 + η))/(1 - ν₀^2)
            K[12,12]=-0.125*(-(0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ + x₂ + x₃ - x₄))*(1 - ξ^2) - 2*η*(1 - ξ)*(-0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ - x₂ + x₃ + x₄)))*(-0.125*E₀*(-(0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ + x₂ + x₃ - x₄))*(1 - ξ^2) - 2*η*(1 - ξ)*(-0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ - x₂ + x₃ + x₄)))/(1 - ν₀^2) + 0.125*E₀*ν₀*(-(1 - η^2)*(0.25*(-y₁ - y₂ + y₃ + y₄) - 0.25*(y₁ - y₂ + y₃ - y₄)) + 2*(0.25*(-y₁ + y₂ + y₃ - y₄) + 0.25*(y₁ - y₂ + y₃ - y₄))*ξ*(1 + η))/(1 - ν₀^2)) + 0.125*(0.125*E₀*(-(1 - η^2)*(0.25*(-y₁ - y₂ + y₃ + y₄) - 0.25*(y₁ - y₂ + y₃ - y₄)) + 2*(0.25*(-y₁ + y₂ + y₃ - y₄) + 0.25*(y₁ - y₂ + y₃ - y₄))*ξ*(1 + η))/(1 - ν₀^2) - 0.125*E₀*ν₀*(-(0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ + x₂ + x₃ - x₄))*(1 - ξ^2) - 2*η*(1 - ξ)*(-0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ - x₂ + x₃ + x₄)))/(1 - ν₀^2))*(-(1 - η^2)*(0.25*(-y₁ - y₂ + y₃ + y₄) - 0.25*(y₁ - y₂ + y₃ - y₄)) + 2*(0.25*(-y₁ + y₂ + y₃ - y₄) + 0.25*(y₁ - y₂ + y₃ - y₄))*ξ*(1 + η)) + (1/2)*E₀*(1 - ν₀)*(-0.125*(-(1 - η^2)*(-0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ - x₂ + x₃ + x₄)) + 2*(0.25*(x₁ - x₂ + x₃ - x₄) + 0.25*(-x₁ + x₂ + x₃ - x₄))*ξ*(1 + η)) + 0.125*(-(0.25*(-y₁ + y₂ + y₃ - y₄) + 0.25*(y₁ - y₂ + y₃ - y₄))*(1 - ξ^2) - 2*η*(1 - ξ)*(0.25*(-y₁ - y₂ + y₃ + y₄) - 0.25*(y₁ - y₂ + y₃ - y₄))))^2/(1 - ν₀^2)

            J[1,1]=-0.25*(1 - η)*x₁ + 0.25*(1 - η)*x₂ + 0.25*(1 + η)*x₃ - 0.25*(1 + η)*x₄
            J[1,2]=-0.25*(1 - η)*y₁ + 0.25*(1 - η)*y₂ + 0.25*(1 + η)*y₃ - 0.25*(1 + η)*y₄
            J[2,1]=-0.25*(1 - ξ)*x₁ + 0.25*(1 - ξ)*x₄ - 0.25*(1 + ξ)*x₂ + 0.25*(1 + ξ)*x₃
            J[2,2]=-0.25*(1 - ξ)*y₁ + 0.25*(1 - ξ)*y₄ - 0.25*(1 + ξ)*y₂ + 0.25*(1 + ξ)*y₃

        return K*det(J)*t
    end
    K=hcubature(BtDB,[-1,-1],[1,1])[1]
    #12x12 to 24x24
    I=1:12
    J=[1,2,6,7,8,12,13,14,18,19,20,24]
    L=sparse(I,J,1.,12,24)
    return L'*K*L
end

function P_GQ12(elm::Quad)::Vector{Float64}
    return zeros(12)
end


text:
"""
    SnpArrays.filter(s)

Filter a SnpArray according to genotyping success rate, minor allele frequencies,
and/or Hardy-Weinberg test.

# Input
- `s`: a SnpArray or Plink file name without the bim, fam, bed suffix.

# Keyword argument
- `min_success_rate_per_row`: Threshold for SNP genotyping success rate. Default 0.98.
- `min_success_rate_per_col`: Threshold for person genotyping success rate. Default 0.98.
- `min_maf`: Minimum minor allele frequency. Default 0.01.
- `min_hwe_pval`: Minimum p-value for Hardy-Weinberg test. Default 0 (not filter HWE).
- `maxiters`: Maximum number of filtering iterations. Default is 5.
- `verbose`: Show progress.

# Output
- `rmask`: BitVector indicating rows after filtering.
- `cmask`: BitVector indicating columns after filtering.
"""
function filter(
    s::SnpArray;
    min_success_rate_per_row::Real = 0.98,
    min_success_rate_per_col::Real = 0.98,
    min_maf::Real = 0.01,
    min_hwe_pval::Real = 0,
    maxiters::Integer = 5,
    verbose::Bool=false)
    m, n = size(s)
    # row-wise counts of missing genotypes
    rmissings = zeros(Int, m)
    # create temporary missings to add to rmissings if chromosome is counted (not filtered)
    temprmissings = zeros(Int, m)
    # columnwise counts: n00, missing, n01, n11
    cc = zeros(Int, 4)
    rmask, cmask = trues(m), trues(n)
    rmissrate, cmissrate = 1 - min_success_rate_per_row, 1 - min_success_rate_per_col
    for iter in 1:maxiters
        verbose && @info ("iteration $iter")
        # number of remaining rows
        rows = count(rmask)
        # maximum allowed missing genotypes each row and col
        cmisses = cmissrate * rows
        fill!(rmissings, 0)
        @inbounds for j in 1:n
            verbose && j % 1000 == 0 && @info "filtering: $j-th of $n SNPs"
            cmask[j] || continue
            # accumulate row and column counts
            fill!(cc, 0)
            fill!(temprmissings, 0)
            for i in 1:m
                rmask[i] || continue
                sij = s[i, j]
                cc[sij + 1] += 1
                sij == 0x01 && (temprmissings[i] += 1)
            end
            # if too many missing genotypes, filter out
            if cc[2] > cmisses
                cmask[j] = false; continue
            end
            # if maf too low, filter out
            maf = (cc[3] + 2cc[4]) / 2(cc[1] + cc[3] + cc[4])
            maf = maf ≤ 0.5 ? maf : 1 - maf
            if maf < min_maf
                cmask[j] = false; continue
            end
            # if HWE p-value too low, filter out
            if min_hwe_pval > 0 && hwe(cc[1], cc[3], cc[4]) < min_hwe_pval
                cmask[j] = false; continue
            end
            # if snp passes all filters, then count the row missings
            rmissings += temprmissings
        end
        # number of remaining cols
        cols = count(cmask)
        # maximum allowed missing genotypes each row
        rmisses = rmissrate * cols
        # filter rows/samples
        @inbounds for i in 1:m
            rmask[i] = rmask[i] && rmissings[i] < rmisses
        end
        # if no change in filter results, done
        count(cmask) == cols && count(rmask) == rows && break
        iter == maxiters && @warn("success rate not satisfied; consider increase maxiters")
    end
    rmask, cmask
end

# function filter(src::AbstractString; des::AbstractString = src * ".filtered", kwargs...)
#     dirname = splitdir(src)[1]
#     srcbedfile = readdir(glob"src.bed", dirname)[1]
#     srcfamfile = readdir(glob"src.fam", dirname)[1]
#     srcm = makestream(srcfamfile) do stream
#         countlines(stream)
#     end
#     s = SnpArray(srcbedfile, srcm)
#     rowmask, colmask = SnpArrays.filter(s; kwargs...)
#     SnpArrays.filter(src, rowmask, colmask; des=des)
# end

"""
    SnpArrays.filter(src, rowinds, colinds; des = src * ".filtered")

Filter `src` Plink files according to row indices `rowinds` and column indices
`colinds` and write to a new set of Plink files `des`.

# Input
- `src`: source Plink file name without suffix ".bed", ".fam" or ".bim".
- `rowinds`: row indices.
- `colinds`: column indices.

# Keyword arguments
- `des`: output Plink file name; default is `src * ".filtered"`.
"""
function filter(
    src::AbstractString,
    rowinds::AbstractVector{<:Integer},
    colinds::AbstractVector{<:Integer};
    des::AbstractString = src * ".filtered")
    # source bed, fam, bim file names
    dirname, filename = splitdir(src)
    srcbedfile = glob(filename * ".bed", dirname)[1]
    srcbimfile = glob(filename * ".bim", dirname)[1]
    srcfamfile = glob(filename * ".fam", dirname)[1]
    # check source plink files
    isfile(srcbedfile) || throw(ArgumentError("$src.bed file not found"))
    isfile(srcbimfile) || throw(ArgumentError("$src.bim file not found"))
    isfile(srcfamfile) || throw(ArgumentError("$src.fam file not found"))
    # destination bed, fam, bim file names
    desbedfile = replace(srcbedfile, src => des)
    desbimfile = replace(srcbimfile, src => des)
    desfamfile = replace(srcfamfile, src => des)
    # numbers of samples and SNPs in src
    srcm = makestream(srcfamfile) do stream
        countlines(stream)
    end
    srcn = makestream(srcbimfile) do stream
        countlines(stream)
    end
    # create row and column masks
    if eltype(rowinds) == Bool
        rmask = rowinds
    else
        rmask = falses(countlines(src * ".fam"))
        rmask[rowinds] .= true
    end
    if eltype(colinds) == Bool
        cmask = colinds
    else
        cmask = falses(countlines(src * ".bim"))
        cmask[colinds] .= true
    end
    desm, desn = count(rmask), count(cmask)
    # write filtered bed file
    bfsrc = SnpArray(srcbedfile)
    makestream(desbedfile, "w+") do io
        write(io, 0x1b6c)
        write(io, 0x01)
        write(io, Matrix{UInt8}(undef, (desm + 3) >> 2, desn))
    end
    bfdes  = SnpArray(desbedfile, desm, "r+")
    bfdes .= @view bfsrc[rmask, cmask]
    # write filtered fam file
    makestream(desfamfile, "w") do io
        for (i, line) in enumerate(eachline(srcfamfile))
            rmask[i] && println(io, line)
        end
    end
    # write filtered bim file
    makestream(desbimfile, "w") do io
        for (j, line) in enumerate(eachline(srcbimfile))
            cmask[j] && println(io, line)
        end
    end
    # output SnpArray
    bfdes
end

@inline ccdf_chisq_1(x) = gamma_inc(1/2, x/2, 0)[2]
"""
    hwe(n00, n01, n11)

Hardy-Weinberg equilibrium test. `n00`, `n01`, `n11` are counts of homozygotes
and heterozygoes respectively. Output is the p-value of type Float64.
"""

function hwe(n00::Integer, n01::Integer, n11::Integer)
    n = n00 + n01 + n11
    n == 0 && return 1.0
    p0 = (n01 + 2n00) / 2n
    (p0 ≤ 0.0 || p0 ≥ 1.0) && return 1.0
    p1 = 1 - p0
    # Pearson's Chi-squared test
    e00 = n * p0 * p0
    e01 = 2n * p0 * p1
    e11 = n * p1 * p1
    ts = (n00 - e00)^2 / e00 + (n01 - e01)^2 / e01 + (n11 - e11)^2 / e11
    #pval = ccdf(Chisq(1), ts)
    pval = ccdf_chisq_1(ts)

    # TODO Fisher exact test
    return pval
end

"""
    indexin_sorted(v, w)

Same as `indexin_general(v, w)` but assumes `v` and `w` are sorted.
"""
function indexin_sorted(v::Union{AbstractArray, Tuple}, w)
    viter, witer = keys(v), eachindex(w)
    vind, wmask = eltype(viter)[], falses(length(w))
    vy, wy = iterate(viter), iterate(witer)
    if vy === nothing || wy === nothing
        return vind
    end
    viteri, i = vy
    witerj, j = wy
    @inbounds begin
        vi, wj = v[viteri], w[witerj]
        while true
            if isless(vi, wj)
                # advance vi
                vy = iterate(viter, i)
                if vy === nothing
                    break
                end
                viteri, i = vy
                vi        = v[viteri]
            elseif isless(wj, vi)
                # advance wj
                wy = iterate(witer, j)
                if wy === nothing
                    break
                end
                witerj, j = wy
                wj        = w[witerj]
            else
                push!(vind, viteri)
                wmask[witerj] = true
                # advance vi
                vy = iterate(viter, i)
                if vy === nothing
                    break
                end
                viteri, i = vy
                vi        = v[viteri]
                # advance wj
                wy = iterate(witer, j)
                if wy === nothing
                    break
                end
                witerj, j = wy
                wj        = w[witerj]
            end
        end
    end
    return vind, wmask
end

"""
    indexin_general(v, w)

Returns an index vector `vind` and `wmask` such that `v[vind]` is the subset
of `v` that appear in `w` and `wmask` is a bitvector such that `v[vind] .== w[wmask]`.
Repeated matches in `v` will be ignored. Only the first match is kept.

# Examples
```jldoctest
julia> a = ['b', 'c', 'd'];

julia> b = ['a', 'b', 'c'];

julia> SnpArrays.indexin_general(a, b)
([1, 2], Bool[false, true, true])
```
"""
function indexin_general(v::Union{AbstractArray, Tuple}, w)
    if issorted(v) && issorted(w)
        return indexin_sorted(v, w)
    else
        Iv, Iw = sortperm(v), sortperm(w)
        vind, wmask = indexin_sorted(v[Iv], w[Iw])
        invpermute!(wmask, Iw)
        if v[Iv[vind]] == w[wmask]
            return Iv[vind], wmask
        else
            return sort!(Iv[vind]), wmask
        end
    end
end

"""
    SnpArrays.filter(srcbedfile, srcbimfile, srcfamfile, rowinds, colinds; des = src * ".filtered")

Filter Plink files  with .gz format or differently named bim and bed files according to row indices
`rowinds` and column indices `colinds` and write to a new set of Plink files `des`.

# Input
- `srcbedfile`: bed file name with suffix such as .bed or .bed.gz.
- `srcbimfile`: bed file name with suffix such as .bim or .bim.gz.
- `srcfamfile`: bed file name with suffix such as .fam or .fam.gz.
- `rowinds`: row indices.
- `colinds`: column indices.

# Keyword arguments
- `des`: output Plink file name; default is `src * ".filtered"`.
"""
function filter(
    srcbedfile::AbstractString,
    srcbimfile::AbstractString,
    srcfamfile::AbstractString,
    rowinds::AbstractVector{<:Integer},
    colinds::AbstractVector{<:Integer};
    des::Union{Nothing, AbstractString} = nothing)
    srcbed = split(srcbedfile, ".bed")[1]
    srcbim = split(srcbimfile, ".bim")[1]
    srcfam = split(srcfamfile, ".fam")[1]
    # check source plink files
    isfile(srcbedfile) || throw(ArgumentError("$srcbedfile file not found"))
    isfile(srcbimfile) || throw(ArgumentError("$srcbimfile file not found"))
    isfile(srcfamfile) || throw(ArgumentError("$srcfamfile file not found"))
    # destination bed, fam, bim file names
    if des == nothing
        desbedfile = srcbed * ".filtered.bed"
        desbimfile = srcbim * ".filtered.bim"
        desfamfile = srcfam * ".filtered.fam"
    else
        desbedfile = des * ".bed"
        desbimfile = des * ".bim"
        desfamfile = des * ".fam"
    end
    # numbers of samples and SNPs in src
    srcm = makestream(srcfamfile) do stream
        countlines(stream)
    end
    srcn = makestream(srcbimfile) do stream
        countlines(stream)
    end
    # create row and column masks
    if eltype(rowinds) == Bool
        rmask = rowinds
    else
        rmask = falses(srcm)
        rmask[rowinds] .= true
    end
    if eltype(colinds) == Bool
        cmask = colinds
    else
        cmask = falses(srcn)
        cmask[colinds] .= true
    end
    desm, desn = count(rmask), count(cmask)
    # write filtered bed file
    bfsrc = SnpArray(srcbedfile, srcm)
    makestream(desbedfile, "w+") do io
        write(io, 0x1b6c)
        write(io, 0x01)
        write(io, Matrix{UInt8}(undef, (desm + 3) >> 2, desn))
    end
    bfdes  = SnpArray(desbedfile, desm, "r+")
    bfdes .= @view bfsrc[rmask, cmask]
    # write filtered fam file
    makestream(desfamfile, "w") do io
        if endswith(srcfamfile, ".gz")
            gzio = GzipDecompressorStream(open(srcfamfile, "r"))
            for (i, line) in enumerate(eachline(gzio))
                rmask[i] && println(io, line)
            end
        else
            for (i, line) in enumerate(eachline(srcfamfile))
                rmask[i] && println(io, line)
            end
        end
    end
    # write filtered bim file
    makestream(desbimfile, "w") do io
        if endswith(srcbimfile, ".gz")
            gzio = GzipDecompressorStream(open(srcbimfile, "r"))
            for (j, line) in enumerate(eachline(gzio))
                cmask[j] && println(io, line)
            end
        else
            for (j, line) in enumerate(eachline(srcbimfile))
                cmask[j] && println(io, line)
            end
        end
    end
    # output SnpArray
    bfdes
end

ext: jl


text:
function fibsq32()
#
#  test polynomial suggested by Goedecker
#  square of Fibocacci polynomial
#
    n = 32;
    p = fib(n)
    p = conv(p,p)
    z = [-.9665892387165332;
         -.9485092065867169-.1867222182129376*im;
         -.9485092065867169+.1867222182129376*im;
         -.8949223987510440-.3665688458356083*im;
         -.8949223987510440+.3665688458356083*im;
         -.8077651569138829-.5329094943253877*im;
         -.8077651569138829+.5329094943253877*im;
         -.6901871862102008-.6795950061744879*im;
         -.6901871862102008+.6795950061744879*im;
         -.5464383280299697-.8011754158225671*im;
         -.5464383280299697+.8011754158225671*im;
         -.3817161237254617-.8930914814237271*im;
         -.3817161237254617+.8930914814237271*im;
         -.2019801534661303-.9518321713734013*im;
         -.2019801534661303+.9518321713734013*im;
         -.1374081948052525e-1-.9750513807247308*im;
         -.1374081948052525e-1+.9750513807247308*im;
         .1761678077924798-.9616381291177512*im;
         .1761678077924798+.9616381291177512*im;
         .3608186423377573-.9117356125589683*im;
         .3608186423377573+.9117356125589683*im;
         .5334081767292508-.8267064135873833*im;
         .5334081767292508+.8267064135873833*im;
         .6874310990889808-.7090471025782057*im;
         .6874310990889808+.7090471025782057*im;
         .8167660514846921-.5622778345619681*im;
         .8167660514846921+.5622778345619681*im;
         .9156129648696007-.3909102172601782*im;
         .9156129648696007+.3909102172601782*im;
         .9783492503358520-.2007808902007197*im;
         .9783492503358520+.2007808902007197*im;
         1.999999999767169];
    z = [z 2*ones(n)]
    p, PolyZeros(z)
end


text:
module TwoDotsModels
    import TwoDots: game_over!, is_adjecent, Dot, Game
    export TwoDotsModel, run!, game_over!, decide!, search_all, search
    export is_in_bounds, connected_dots, can_connect, color_matches

    include("structs.jl")
    include("functions.jl")
end


text:
mutable struct Data{T<:Integer}
    x::Vector{T}
    n::Vector{T}
    function Data(x::Vector{T}, n::Vector{T}) where T<:Integer
        @assert length(x) == length(n) "Inputs x and n must have the same length."
        @assert all(x .<= n) "xᵢ <= nᵢ for all elements i = 1,...,N"
        @assert all(x .>= 0) "xᵢ >= 0 for all elements i = 1,...,N"
        new{T}(x, n)
    end
end

mutable struct Prior
    μ::Beta
    σ::Beta
end

mutable struct Candidate
    μ::Normal
    σ::Normal
end

struct Chain{T<:AbstractFloat}
    μ::Vector{T}
    σ::Vector{T}
    θ::Matrix{T}
    acceptance::Matrix{<:Integer}
end


text:
using LinearAlgebra

function householder(a,b)
    n = size(a,1);

    w = a-b;
    Id = Matrix{Float64}(I,n,n);
    Q =  Id - 2(w*w')/dot(w,w);

    return Q;
end


text:
function intice()
    # Nothing because weebs are already attracted to underage anime girls
end

function callthefbi(age)
    if age < 15
        run(`call_911`)
    else
        println("Probably just something weird happened")
    end
end

println("What\'s the girl\'s age?")
age = parse(Int32, readline())
println(age)
callthefbi(age)

max_stars_repo_name: REWERK/lolicon

max_issues_repo_name: seaguli/lollipop

max_forks_repo_name: seaguli/lollipop


text:
#
# This file is part of the Guards.jl Julia package,
# MIT license, part of https://github.com/JuliaActors
#

using Guards, Test
import Actors: newLink

me = newLink()
gd = guard([1,2,3])
@test gd isa Guards.Guard{Array{Int64,1}}
@test call(gd) == [1,2,3]
call(gd, me)
@test receive(me).y == [1,2,3]
@test push!(call(gd), 4) == [1,2,3,4]
@test call(gd) == [1,2,3]
@test call(gd, push!, 4) == [1,2,3,4]
@test call(gd) == [1,2,3,4]
call(gd, me, push!, 5)
@test receive(me).y == [1,2,3,4,5]
cast(gd, push!, 6)
@test (@grd gd) == [1,2,3,4,5,6]
@test (@grd pop!(gd)) == 6
@test pop!(update!(gd, [6,7,8])) == 8
@test (@grd gd) == [6,7,8]

@test (@grd 1) == "@grd: not a symbol or a call!"
@test_throws AssertionError update!(gd, [10.0, 11.0])


text:
using Test
using NQCDynamics
using NQCDynamics.InitialConditions
using Random
using Distributions
using HDF5
using ComponentArrays

a = [rand(Float64, 3, 2) for i=1:100]
d = DynamicalDistribution(a, a, (3,2))
@test (3,2) == size(d)
@test rand(d) isa ComponentVector
@test lastindex(d) == 100

d = DynamicalDistribution(a, Normal(), (3,2))
@test rand(d) isa ComponentVector
@test lastindex(d) == 100

d = DynamicalDistribution(1, Normal(), (3,2))
@test rand(d) isa ComponentVector
@test lastindex(d) == 1

@testset "BoltzmannVelocityDistribution" begin
    boltz = NonadiabaticDistributions.BoltzmannVelocityDistribution(1, [1, 2, 3], (2,3))
    @test size(rand(boltz)) === boltz.size
    boltz = NonadiabaticDistributions.BoltzmannVelocityDistribution(1, [1, 2, 3], (2,3,3))
    @test size(rand(boltz)) === boltz.size
end

@testset "HDF5 distribution" begin
    filename = "dat.h5"
    NonadiabaticDistributions.write_hdf5(filename, a)
    @test NonadiabaticDistributions.select_item(filename, 1, (3,2)) ∈ a

    d = DynamicalDistribution(filename, 1, (3,2))
    u = rand(d)
    @test all(u.r .== 1)
    @test u.v ∈ a
    @test lastindex(d) == 100
    rm(filename)
end


text:
"""
    deleteRows(x, rowss)

Delete rows from a 2x2 matrix. Use Float or Int array for x.
"""
function deleteRows(xIn::Array{Float64, 2}, rows::Array{Int64, 1})

    # get dimensions of input
    (numRows, numCols) = size(xIn)

    # for calculating new dimensions
    numRowDelete = length(rows)

    # uninitialized array with new dimensions
    xOut = Array{Float64}(undef, (numRows - numRowDelete), numCols)

    # dynamically build new matrix without columns to delete
    xIndexIn  = 1
    xIndexOut = 1
    for ii in rows
        xRowsIn =  xIndexIn : (ii - 1)
        xRowsOut =  xIndexOut : xIndexOut + length(xRowsIn) - 1

        # setting xOut cols to xIn, xColsOut always <= xColsIn
        xOut[xRowsOut, :] = xIn[xRowsIn, :]

        xIndexIn = ii + 1
        xIndexOut = xIndexOut + length(xRowsIn)
    end

    # need one more iteration to get last section of cols
    if xIndexIn <= numRows
        xRowsIn = xIndexIn : numRows
        xRowsOut = xIndexOut : xIndexOut + length(xRowsIn) - 1

        xOut[xRowsOut, :] = xIn[xRowsIn, :]
    end

    return xOut
end


text:
# Test the gaussian pre-image code
# Danny Perry (dperry@cs.utah.edu)
# May 2015

using Debug
using PyPlot
using Kernel

@debug function main()

	n = 100
	rad = 10
	thetas = 2*pi*rand(n)
	X = rad*[cos(thetas) sin(thetas)]
	Xtest = mean(X,1)

	figure()
	plot(X[:,1],X[:,2],"ro")
	hold(:on)
	plot(Xtest[1],Xtest[2],"bo")

	iters = 1
	neibsize = 2
	sigma = 10.0
	#result,neibs = GaussianKernelDenoise(X, Xtest, iters, neibsize, sigma)
	K = GaussianKernel(X,sigma)
	Ktest = GaussianKernel(X,Xtest,sigma)
	result,neibs = GaussianKernelPreImage(X, K, Ktest, neibsize, sigma)
	println("==================")
	println("neibsize: ", neibsize)
	println("input: ", Xtest)
	println("result: ", result)
	println("diff: ", norm(Xtest-result))
	println("==================")

	figure()
	plot(X[:,1],X[:,2],"ro")
	hold(:on)
	plot([Xtest[1],result[1]],[Xtest[2],result[2]],"bx--")
	plot(result[1],result[2],"gs")
	plot(X[vec(neibs[1,:]),1],X[vec(neibs[1,:]),2],"y^")
	title(@sprintf("%d",neibsize))

	neibsize = 4
	result,neibs = GaussianKernelDenoise(X, Xtest, iters, neibsize, sigma)
	println("==================")
	println("neibsize: ", neibsize)
	println("input: ", Xtest)
	println("result: ", result)
	println("diff: ", norm(Xtest-result))
	println("==================")

	figure()
	plot(X[:,1],X[:,2],"ro")
	hold(:on)
	plot([Xtest[1],result[1]],[Xtest[2],result[2]],"bx--")
	plot(result[1],result[2],"gs")
	plot(X[vec(neibs[1,:]),1],X[vec(neibs[1,:]),2],"y^")
	title(@sprintf("%d",neibsize))

	neibsize = 10
	result,neibs = GaussianKernelDenoise(X, Xtest, iters, neibsize, sigma)
	println("==================")
	println("neibsize: ", neibsize)
	println("input: ", Xtest)
	println("result: ", result)
	println("diff: ", norm(Xtest-result))
	println("==================")

	figure()
	plot(X[:,1],X[:,2],"ro")
	hold(:on)
	plot([Xtest[1],result[1]],[Xtest[2],result[2]],"bx--")
	plot(result[1],result[2],"gs")
	plot(X[vec(neibs[1,:]),1],X[vec(neibs[1,:]),2],"y^")
	title(@sprintf("%d",neibsize))

	theta = rand()*2*pi
	Xtest = rad*[cos(theta) sin(theta)]
	neibsize = 2
	result,neibs = GaussianKernelDenoise(X, Xtest, iters, neibsize, sigma)
	println("neibsize: ", neibsize)
	println("input: ", Xtest)
	println("result: ", result)
	println("diff: ", norm(Xtest-result))

	figure()
	plot(X[:,1],X[:,2],"ro")
	hold(:on)
	plot([Xtest[1],result[1]],[Xtest[2],result[2]],"bx--")
	plot(result[1],result[2],"gs")
	plot(X[vec(neibs[1,:]),1],X[vec(neibs[1,:]),2],"y^")
	title(@sprintf("%d",neibsize))

	neibsize = 4
	result,neibs = GaussianKernelDenoise(X, Xtest, iters, neibsize, sigma)
	println("neibsize: ", neibsize)
	println("input: ", Xtest)
	println("result: ", result)
	println("diff: ", norm(Xtest-result))

	figure()
	plot(X[:,1],X[:,2],"ro")
	hold(:on)
	plot([Xtest[1],result[1]],[Xtest[2],result[2]],"bx--")
	plot(result[1],result[2],"gs")
	plot(X[vec(neibs[1,:]),1],X[vec(neibs[1,:]),2],"y^")
	title(@sprintf("%d",neibsize))

	neibsize = 10
	result,neibs = GaussianKernelDenoise(X, Xtest, iters, neibsize, sigma)
	println("neibsize: ", neibsize)
	println("exact: ", Xtest)
	println("result: ", result)
	println("diff: ", norm(Xtest-result))

	figure()
	plot(X[:,1],X[:,2],"ro")
	hold(:on)
	plot([Xtest[1],result[1]],[Xtest[2],result[2]],"bx--")
	plot(result[1],result[2],"gs")
	plot(X[vec(neibs[1,:]),1],X[vec(neibs[1,:]),2],"y^")
	title(@sprintf("%d",neibsize))

	PyPlot.plt.show()
	@bp

end

function testdata()
	X = reshape([1:33],11,3)
	Xtest = X[5,:]
	X = X[[1:4,6:end],:]

	println("Xtest: ", Xtest)
	println("X:\n", X)

	return X,Xtest
end

@debug function simple(neibsize = 2)
	println("~~~~~~~~~~~~~~~~~~~")
	println(" neighbors approach")
	println("~~~~~~~~~~~~~~~~~~~")

	X,Xtest = testdata()

	sigma = 1.0
	K = GaussianKernel(X,sigma)
	Ktest = GaussianKernel(X,Xtest,sigma)
	result,neibs = GaussianKernelPreImage(X, K, Ktest, neibsize, sigma)

	println("neibs: ", neibs)

	println("===================")
	println("neibsize: ", neibsize)
	println("exact: ", Xtest)
	println("result: ", result)
	println("diff: ", norm(Xtest-result))
	println("===================")

	return result
end

@debug function simple2()
	println("~~~~~~~~~~~~~~~~~~~")
	println(" fixed point approach")
	println("~~~~~~~~~~~~~~~~~~~")

	X,Xtest = testdata()

	sigma = 1.0
	K = GaussianKernel(X,sigma)
	U,S,V = svd(K)
	Ktest = GaussianKernel(X,Xtest,sigma)
	y = vec(Ktest' * V)
	z_init = mean(X,1)
	result,iters = GaussianKernelPreImage(y,V, X,z_init, sigma, 1e-3,100)

	println("===================")
	println("exact: ", Xtest)
	println("iters: ", iters)
	println("result: ", result)
	println("diff: ", norm(Xtest-result))
	println("===================")
end

@debug function simple3()
	println("~~~~~~~~~~~~~~~~~~~")
	println(" both approach")
	println("~~~~~~~~~~~~~~~~~~~")

	X,Xtest = testdata()

	neibsize = 3
	sigma = 1.0
	K = GaussianKernel(X,sigma)
	U,S,V = svd(K)
	Ktest = GaussianKernel(X,Xtest,sigma)
	y = vec(Ktest' * V)
	z_init = simple(neibsize)
	println("zinit: ", z_init)
	println("diff: ", norm(Xtest-z_init))
	result,iters = GaussianKernelPreImage(y,V, X,z_init, sigma, 1e-20,10000)

	println("===================")
	println("neibsize: ", neibsize)
	println("exact: ", Xtest)
	println("iters: ", iters)
	println("result: ", result)
	println("diff: ", norm(Xtest-result))
	println("===================")
end

simple() # method 1
simple2() # method 2
simple3() # method 1 and 2
#main()


text:
module ParameterisedModule
using MLStyle

export @sig, @structure, @open, mk_signature, mk_structure, open_module, module_get, module_types, module_values

abstract type MLModule end

extract_tvar(var :: Union{Symbol, Expr})::Symbol =
    @match var begin
        :($a <: $_) => a
        :($a >: $_) => a
        :($_ >: $a >: $_) => a
        :($_ <: $a <: $_) => a
        a::Symbol         => a
    end

function module_get end
function module_values end
function module_types end

function mk_signature(__source__, __module__, ex)
    @when :(struct $modtype; $(decls...) end) = ex begin
        typedecls = Symbol[]
        valsigs = Tuple{Symbol, Any}[]
        line = __source__
        (modn, modps) = @match modtype begin
            :($modn{$(modps...)}) => (modn, modps)
            modn::Symbol => (modn, [])
            a => error("invalid signature name $a")
        end
        foreach(decls) do decl
            @match decl begin
                :(struct $(n::Symbol); $(_...) end) => push!(typedecls,  n)
                :($x :: $t) =>
                    if x isa Symbol
                        push!(valsigs, (x, t))
                    else
                        error("invalid value declaration $x, expect a Symbol.")
                    end
                l::LineNumberNode => (line = l)
                a => error("invalid declaration $a")
            end
        end

        valns   = map(x -> x[1], valsigs)
        valts   = map(x->x[2], valsigs)
        ex_valts = Expr(:curly, :Tuple, valts...)
        use_modps = map(extract_tvar, modps)

        val_type_conflicts = [n in typedecls for n in valns]

        any(val_type_conflicts) &&
            error("Name conflicts: names $val_type_conflicts in both types and values.")

        @gensym Values
        @gensym values

        fieldgetter = Expr[]

        append!(fieldgetter, [
            :(
                $ParameterisedModule.module_get($values :: $modn, ::Val{$(QuoteNode(n))}) =
                    getfield($values, $(QuoteNode(values)))[$i]
            )
            for (i, n) in enumerate(valns)
        ])

        append!(fieldgetter, [
            :(
                $ParameterisedModule.module_get(
                    $values :: $modn{$(use_modps...),
                                     $(typedecls...)},
                    $Values::Val{$(QuoteNode(n))}) where {$(use_modps...), $(typedecls...)} = $n
            )
            for n in typedecls
        ])
        quote
            $__source__

            struct $modn{
                        $(modps...),
                        $(typedecls...),
                        $Values <: $ex_valts
                    } <: $MLModule
                $values :: $Values
            end

            $ParameterisedModule.module_types(::Type{$modn}) = $(Tuple(typedecls))
            $ParameterisedModule.module_values(::Type{$modn}) = $(Tuple(valns))

            Base.getproperty($values :: $modn, $Values::Symbol) =
                $ParameterisedModule.module_get($values, Val{$Values}())

            $(fieldgetter...)
        end
    @otherwise
        error("invalid signature")
    end
end

function mk_structure(__source__, __module__, ex)
    @when :(struct $modtype; $(decls...) end) = ex begin
        line = __source__
        (modn, modps) = @match modtype begin
            :($modn{$(modps...)}) => (modn, modps)
            modn::Symbol => (modn, [])
            a => error("invalid signature name $a")
        end
        mod = __module__.eval(modn)
        type_ns = module_types(mod)
        value_ns = module_values(mod)
        body = Any[]
        for decl in decls
            @match decl begin
                Expr(:abstract, _...) || Expr(:struct, _...) => __module__.eval(decl)
                _ => push!(body, decl)
            end
        end

        vals = Expr(:tuple, value_ns...)
        @gensym val
        push!(body,
            :(
                let $val = $vals
                    $mod{$(modps...), $(type_ns...), typeof($val)}($val)
                end
            ))
        Expr(:let, Expr(:block), Expr(:block, body...))
    @otherwise
        error("invalid structure")
    end
end

macro sig(ex)
    __module__.eval(mk_signature(__source__, __module__, ex))
end

macro structure(ex)
    esc(mk_structure(__source__, __module__, ex))
end

function open_module(::Type{M}, m) where M <: MLModule
    value_ns = module_values(M)
    type_ns = module_types(M)
    ret = Expr[]
    for each in value_ns
        push!(ret, :($each = $m.$each))
    end
    for each in type_ns
        push!(ret, :($each = $m.$each))
    end
    ret
end

macro open(mtype, m)
    esc(Expr(:block, open_module(__module__.eval(mtype), m)...))
end

macro open(mtype, m, inner)
    esc(
        Expr(:let,
            Expr(:block, open_module(__module__.eval(mtype), m)...),
            inner
        ))
end

end # module


text:
# Copyright 2021 Google LLC
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
#
# Executes the identifiability experiement with some instances of the netlib
# database.
# To run this script execute
# $ julia --project run_netlib_identifiability_experiment.jl <path_where_to_store_results> <path_to_directory_containing_mps_files>
#
# By default the script runs a specific set of twelve instances. The script also
# optionally accepts a specific instance to run and whether to generate the
# iteration plot based on assuming differences are better than average and
# current iterates or not, for example,
# $ julia --project run_netlib_identifiability_experiment.jl <path_where_to_store_results> <path_to_directory_containing_mps_files> box1 true
# Note that you must specify either both the instance name and the boolean, or
# neither of them.

import FirstOrderLp

include("primal_dual_hybrid_gradient.jl")
include("utils.jl")

function run_netlib_instance(instance_name::String, difference_is_better::Bool)
  params = PdhgParameters(
    nothing, # step_size (forces the solver to use a provably correct step size)
    1.0, # primal_weight
    1000, # printing frequency
    false, # verbose
    1000000, # iteration limit
    nothing, # two_dimensional_subspace
    0.0, # eps_primal_infeasible
    0.0, # eps_dual_infeasible
    false, # store_iterates
    nothing, # reference_point
  )
  print("About to start solving " * instance_name * "\n")
  instance_path = joinpath(ARGS[2], (instance_name * ".mps"))
  lp = FirstOrderLp.qps_reader_to_standard_form(instance_path)
  output = optimize(params, lp)
  factor = difference_is_better ? 10 : 2
  params = PdhgParameters(
    nothing, # step_size (forces the solver to use a provably correct step size)
    1.0, # primal_weight
    1000, # printing frequency
    false, # verbose
    min(1000000, max(2000, factor * output.last_active_set_change)), # iteration limit
    nothing, # two_dimensional_subspace
    0.0, # eps_primal_infeasible
    0.0, # eps_dual_infeasible
    false, # store_iterates
    [output.primal_delta; output.dual_delta], # reference_point
  )
  output_with_distances = optimize(params, lp)
  generate_infeasibility_convergence_plot(
    output_with_distances,
    joinpath(ARGS[1], instance_name * ".pdf"),
  )
end

function main()
  if length(ARGS) == 4
    run_netlib_instance(ARGS[3], parse(Bool, ARGS[4]))
  else
    for difference_is_better in [false, true]
      if difference_is_better
        instance_names =
          ["box1", "galenet", "itest2", "itest6", "woodinfe", "ex72a", "ex73a"]
      else
        instance_names = ["bgdbg1", "bgindy", "cplex1", "chemcom", "mondou2"]
      end

      for instance_name in instance_names
        run_netlib_instance(instance_name, difference_is_better)
      end
    end
  end
end

main()


text:
using HorizonSideRobots
include("../basics.jl")

md"""
            На прямоугольном поле произвольных размеров расставить маркеры в виде "шахматных" клеток, начиная с юго-западного угла поля, когда каждая
            отдельная "шахматная" клетка имеет размер n x n клеток поля (n - это параметр функции). Начальное положение Робота - произвольное,
            конечное - совпадает с начальным. Клетки на севере и востоке могут получаться "обрезанными" - зависит от соотношения размеров поля
            и "шахматных" клеток. (Подсказка: здесь могут быть полезными две глобальных переменных, в которых будут содержаться текущие декартовы
            координаты Робота относительно начала координат в левом нижнем углу поля, например)
"""

 """Решение задачи 12"""
 function solve_12!(r, n)
    (M, N), (I, J) = get_field_size_and_position!(r)
    cells = Set()
    visited = Set()
    count = 0

    # Отбирает клетки поля подходящие для раскраски и пушит в {cells}
    for t in 1:M
        i = 1 + M - t
        if (M-i) % (2*n) < n
            count = 0
        else
            count = n
        end
        for j in 1:N
            if count < n
                push!(cells, (i, j))
            elseif count == 2*n
                push!(cells, (i, j))
                count = 0
            end
            count += 1
        end
    end

    mark_some_cells_with_dfs!(r, visited, (I, J), false, cells)

    show(r)
end


text:
"""
Fast (non-allocating) version of norm(x-y,2)^2
"""
function normdiff2(x,y)
    s = 0.0
    for i in eachindex(x)
        s += abs2(x[i]-y[i])
    end
    return s
end

"""
Fast (non-allocating) version of norm(x-y,2)
"""
normdiff(x,y) = sqrt(normdiff2(x,y))


text:
using Documenter, Arpack

makedocs(
    format = Documenter.HTML(
        canonical = "https://julialinearalgebra.github.io/Arpack.jl/stable/",
    ),
    sitename = "Arpack.jl",
    modules = [Arpack],
    pages = [
        "Home" => "index.md",
        "Standard Eigen Decomposition" => "eigs.md",
        "Generalized Eigen Decomposition" => "eigs_gen.md",
        "Singular Value Decomposition" => "svds.md",
        "API Reference" => "api.md",
    ]
)

deploydocs(
    repo = "github.com/JuliaLinearAlgebra/Arpack.jl.git",
    target = "build",
    deps = nothing,
    make = nothing
)


text:
import Base:   ==
import ACE: read_dict, write_dict,
       transform, transform_d, transform_dd, inv_transform

export polytransform, morsetransform, agnesitransform

abstract type DistanceTransform end

# ----- new transforms implementation
import ACE: λ

@deprecate PolyTransform(p, r0) polytransform(p, r0)
polytransform(p, r0) = λ("r -> ((1+$r0)/(1+r))^$p")

@deprecate IdTransform() idtransform()
idtransform() = λ("r -> r")

@deprecate MorseTransform(lambda, r0) morsetransform(lambda, r0)
morsetransform(lambda, r0) = λ("r -> exp(- $lambda * (r / $r0 - 1))")

@deprecate AgnesiTransform(args...) agnesitransform(args...)
agnesitransform(r0, p=2, a=(p-1)/(p+1)) = λ("r -> 1 / (1 + $a * (r / $r0)^$p)")

# ------------------------------------------------------
#  implementation of inverse transform

import Roots: find_zero
function inv_transform(trans, x)
   # solve trans(r) = x
   r = find_zero(r -> trans(r) - x, 1.0)
   if !(trans(r) ≈ x)
      @warn("inv_transform via find_zero didn't find a good solution")
   end
   return r
end

# ------------------------------------------------------
# generic ad codes for distance transforms

import ACE: rrule_evaluate, frule_evaluate, evaluate

evaluate(trans::DistanceTransform, r::Number) = transform(trans, r)

function frule_evaluate(trans::DistanceTransform, r::Number, ∇r)
   # ∇r will be an abstract array or an abstract array encapsulated in a state
   dt = transform_d(trans, r)
   return dt * ∇r
end

function rrule_evaluate(trans::DistanceTransform, r::Number, w::AbstractVector)
   # here, the typical scenario is that evaluate -> B(T(r)), so the rrule
   # should give ∑_n u_n d/dr B_n(T(r)) = ∑_n u_n B_n'(T(r)) * T'(r)
   # i.e. w = ∑_n u_n B_n'(T(r)) which is scalar-like.
   return transform_d(trans, r) * w
end


text:
using jInv.Mesh
using jInv.Utils
using Test

# setup 3D mesh
nc = [5, 7, 8]
x0 = rand(3)
domain = [x0[1], 4, x0[2], 2, x0[3], 6]
h   = (domain[2:2:end]-domain[1:2:end])./nc
h1  = h[1]*ones(nc[1])
h2  = h[2]*ones(nc[2])
h3  = h[3]*ones(nc[3])

Mt = getTensorMesh3D(h1,h2,h3,x0)
Mt2 = getTensorMesh3D(h1,h2,h3,x0)
Mr = getRegularMesh(domain,nc)
Mr2 = getRegularMesh(domain,nc)

# test if Mt and Mt2 are equal
@test Mt == Mt2
@test Mr == Mr2

Meshes = (Mt,Mt2,Mr,Mr2)

function buildOps!(M)
	D = getDivergenceMatrix(M)
	G = getNodalGradientMatrix(M)
	C = getCurlMatrix(M)
	Af = getFaceAverageMatrix(M)
	Ae = getEdgeAverageMatrix(M)
	An = getNodalAverageMatrix(M)
	V  = getVolume(M)
	F  = getFaceArea(M)
	L  = getLength(M)
	Vi = getVolumeInv(M)
	Fi = getFaceAreaInv(M)
	Li = getLengthInv(M)
	Lap = getNodalLaplacianMatrix(M)
end

# first all fields should be empty
for k=length(Meshes)
	@test isempty(Meshes[k].Div )
	@test isempty(Meshes[k].Grad)
	@test isempty(Meshes[k].Curl)
	@test isempty(Meshes[k].Af  )
	@test isempty(Meshes[k].Ae  )
	@test isempty(Meshes[k].An  )
	@test isempty(Meshes[k].V   )
	@test isempty(Meshes[k].F   )
	@test isempty(Meshes[k].L   )
	@test isempty(Meshes[k].Vi  )
	@test isempty(Meshes[k].Fi  )
	@test isempty(Meshes[k].Li  )
	@test isempty(Meshes[k].nLap)
end

# build all matrices
buildOps!(Mt)
buildOps!(Mt2)
buildOps!(Mr)
buildOps!(Mr2)
@test Mt == Mt2
@test Mr == Mr2

# now all fields should be non-empty
for k=length(Meshes)
	@test !isempty(Meshes[k].Div )
	@test !isempty(Meshes[k].Grad)
	@test !isempty(Meshes[k].Curl)
	@test !isempty(Meshes[k].Af  )
	@test !isempty(Meshes[k].Ae  )
	@test !isempty(Meshes[k].An  )
	@test !isempty(Meshes[k].V   )
	@test !isempty(Meshes[k].F   )
	@test !isempty(Meshes[k].L   )
	@test !isempty(Meshes[k].Vi  )
	@test !isempty(Meshes[k].Fi  )
	@test !isempty(Meshes[k].Li  )
	@test !isempty(Meshes[k].nLap)
end

clear!(Mt)
clear!(Mt2)
clear!(Mr)
clear!(Mr2)

@test Mt == Mt2
@test Mr == Mr2

# now all matrices should be gone
for k=length(Meshes)
	@test isempty(Meshes[k].Div )
	@test isempty(Meshes[k].Grad)
	@test isempty(Meshes[k].Curl)
	@test isempty(Meshes[k].Af  )
	@test isempty(Meshes[k].Ae  )
	@test isempty(Meshes[k].An  )
	@test isempty(Meshes[k].V   )
	@test isempty(Meshes[k].F   )
	@test isempty(Meshes[k].L   )
	@test isempty(Meshes[k].Vi  )
	@test isempty(Meshes[k].Fi  )
	@test isempty(Meshes[k].Li  )
	@test isempty(Meshes[k].nLap)
end


text:
mutable struct Line <: Shape
    position::Tuple{Luxor.Point,Luxor.Point}
    color::Union{Blending,Coloring}
    style::Union{Symbol, Vector{Float64}}
end

function Line( ; position, color, style=:stroke, kwargs...)
    return Line(position, color, style)
end


text:
using AnyMOD, Gurobi

# create model
anyM = anyModel("testModel1","results")
createOptModel!(anyM)
setObjective!(:costs,anyM)

# analyze created model
plotTree(:region,anyM)
plotEnergyFlow(:graph,anyM)

# solve model
set_optimizer(anyM.optModel,Gurobi.Optimizer) # select a solver
optimize!(anyM.optModel) # solve the model

# analyze results
reportResults(:summary,anyM)
reportTimeSeries(:electricity, anyM)

# create sankey plots
plotEnergyFlow(:sankey,anyM)


text:
"""
    AbstractPopulation

Abstract base type for population.

Keeps track of which generation it is and allows for iterating over its members.
"""
abstract type AbstractPopulation end

generation(p::AbstractPopulation) = generation(p.base)

Base.iterate(p::AbstractPopulation, s...) = Base.iterate(wrappedpop(p), s...)

Base.length(p::AbstractPopulation) = length(wrappedpop(p))
Base.eltype(p::AbstractPopulation) = eltype(wrappedpop(p))
Base.size(p::AbstractPopulation) = size(wrappedpop(p))
Base.IteratorSize(p::AbstractPopulation) = Base.IteratorSize(wrappedpop(p))
Base.IteratorEltype(p::AbstractPopulation) = Base.IteratorEltype(wrappedpop(p))
Base.keys(p::AbstractPopulation) = keys(wrappedpop(p))
Base.getindex(p::AbstractPopulation, args...) = getindex(wrappedpop(p), args...)

generation_filename(dir) = joinpath(dir, "generation.txt")

"""
    struct Population{N, P}
    Population(gen, members)
    Population(members)

Basic population type which just adds generation counting to its members.

Evolving the population returns a new `Population` with the generation counter incremented.
"""
struct Population{N<:Integer, P} <: AbstractPopulation
    gen::N
    members::P
end
Population(members) = Population(1, members)
Population(members::PersistentArray) = Population(members.savedir, members)
function Population(savedir::AbstractString, members)
    gfile = generation_filename(savedir)
    gen = !isfile(gfile) ? 1 : parse(Int, readline(gfile))
    return Population(gen, members)
end

wrappedpop(p::Population) = p.members

generation(p::Population) = p.gen

fitness(p::Population, f::AbstractFitness) = Population(p.gen, map(wrappedpop(p)) do cand
    FittedCandidate(cand, f, p.gen)
end)
evolve(e::AbstractEvolution, p::Population) = Population(p.gen + 1, evolve(e, wrappedpop(p)))

"""
    evolve(e::AbstractEvolution, f::AbstractFitness, p::Population)

Return a new population fitted by `f` and evolved by `e`.

This is done by first replacing each member of `p` with a [`FittedCandidate`](@ref) with fitness computed by `f`.

Then evolve `population` into a new population using `e`. New population may or may not contain same individuals as before.
"""
evolve(e::AbstractEvolution, f::AbstractFitness, p::Population) = evolve(e, fitness(p, f))

function persist(p::Population{N, <:PersistentArray}) where N
    persist(wrappedpop(p))
    open(io -> write(io, string(p.gen)), generation_filename(wrappedpop(p).savedir); write=true)
end


text:
using JuLIP.Testing, ASE, PyCall

@info("These tests to compare JuLIP vs ASE implementations of some potentials")

h3("Compare JuLIP vs ASE: EAM")
# JuLIP's EMT implementation
at = set_pbc!( bulk(:Cu, cubic=true) * (2,2,2), (true,false,false) )
rattle!(at, 0.02)
emt = EMT(at)

@info("Test JuLIP vs ASE EMT implementation")
pyemt = ASE.Models.EMTCalculator()
print("   energy: ")
println(@test abs(energy(emt, at) - energy(pyemt, at)) < 1e-10)
print("   forces: ")
println(@test norm(forces(pyemt, at) - forces(emt, at), Inf) < 1e-10)

# ------------------------------------------------------------------------

h3("Compare JuLIP EAM Implementation against ASE EAM Implementation")

@pyimport ase.calculators.eam as eam

pot_file = joinpath(dirname(pathof(JuLIP)), "..", "data", "w_eam4.fs")

@info("Generate the ASE potential")
eam4_ase = eam.EAM(potential=pot_file) |> ASECalculator

@info("Generate low-, med-, high-accuracy JuLIP potential")
eam4_jl1 = EAM(pot_file)
eam4_jl2 = EAM(pot_file; s = 1e-4)
eam4_jl3 = EAM(pot_file; s = 1e-6)

at1 = rattle!(bulk(:W, cubic=true) * 3, 0.1)
at2 = deleteat!(bulk(:W, cubic=true) * 3, 1)
at1_ase = ASEAtoms(at1)
at2_ase = ASEAtoms(at2)

for (i, (at, at_ase)) in enumerate(zip([at1, at2], [at1_ase, at2_ase]))
   @info("Test $i")
   err_low = (energy(eam4_ase, at_ase) - energy(eam4_jl1, at)) / length(at)
   err_med = (energy(eam4_ase, at_ase) - energy(eam4_jl2, at)) / length(at)
   err_hi = (energy(eam4_ase, at_ase) - energy(eam4_jl3, at)) / length(at)
   print("   Low Accuracy energy error:", err_low, "; ",
         (@test abs(err_low) < 0.02))
   print("   Low Accuracy energy error:", err_med, "; ",
         (@test abs(err_med) < 0.006))
   print("   Low Accuracy energy error:", err_hi, "; ",
         (@test abs(err_hi) < 0.0004))
end

# ------------------------------------------------------------------------