Voigt

import LinAlg.norm

abstract AbstractVoigt{T}

immutable VoigtStress <: AbstractVoigt
  xx::Float64
  yy::Float64
  zz::Float64
  yz::Float64
  xz::Float64
  xy::Float64
end

immutable VoigtStrain <: AbstractVoigt
  xx::Float64
  yy::Float64
  zz::Float64
  yz::Float64
  xz::Float64
  xy::Float64
end

# Double contraction
function (*){T <: AbstractVoigt, P <: AbstractVoigt}(s::T, e::P)
    return sqrt(s.xx*e.xx + s.yy*e.yy + s.zz*e.zz +
                s.yz*e.yz + s.xz*e.xz + s.xy*e.xy)
end

function (*){T <: AbstractVoigt}(voigt::T, vec::AbstractArray)
    tmp = typevec!(voigt) * vec
    T(tmp...)
end

function (*){T <: AbstractVoigt}(vec::AbstractArray, voigt::T)
    tmp = vec * typevec!(voigt)
    T(tmp...)
end


function dev{T <: AbstractVoigt}(s::T)
  xx = 1/3 * (2s.xx - s.yy - s.zz)
  yy = 1/3 * (2s.yy - s.xx - s.zz)
  zz = 1/3 * (2s.zz - s.xx - s.yy)
  return T(xx, yy, zz, s.yz, s.xz, s.xy)
end

function norm(s::VoigtStress)
return sqrt(s.xx*s.xx + s.yy*s.yy + s.zz*s.zz +
            s.yz*s.yz + s.xz*s.xz + s.xy*s.xy)
end

function norm(e::VoigtStrain)
return sqrt(e.xx*e.xx + e.yy*e.yy + e.zz*e.zz +
            0.25*(e.yz*e.yz + e.xz*e.xz + e.xy*e.xy))
end

function vonmises(s::VoigtStress)
  s_dev = dev(s)
  return sqrt(3/2) * (s_dev * s_dev)
end


typevec!(v::AbstractVoigt) = pointer_to_array(
    convert(Ptr{typeof(v.xx)}, pointer_from_objref(v))+sizeof(typeof(v.xx)), 6)

s = VoigtStress(1.0,2.0, 3.0,1.0,2.0, 3.0)
e = VoigtStrain(1.0,2.0, 3.0,1.0,2.0, 3.0)