module Ela

    importall Mat, Geo, Sec, Ele, Bc

    export Domain, currentDomain,
            addNode, addRefpoint, addMat, addSec, 
            addEle, addFix, addCons, addLoad, 
            setTimeSeries,
            buildModel, saveModel,
            assembleSystemMatrix,
            # assembleStiffMatrix, 
            # assembleMassMatrix, 
            # assembleDampMatrix,
            applyCons, applyLoad,
            getSolution,
            getReactionAndStiff, getReaction, 
            updateHistoryPara,updateStiffHistPara

    typealias R Float64
    global currentDomain

    type LinearAlgebraSystem{S<:Integer}

        nDof::S
        K::SparseMatrixCSC{R}  # 刚度矩阵
        M::SparseMatrixCSC{R}  # 质量矩阵
        C::SparseMatrixCSC{R}  # 阻尼矩阵
        dP::Vector{R}  # 增量荷载向量
        dU::Vector{R}  # 不平衡力向量
        UF::Vector{R}  # 不平衡力向量

        function LinearAlgebraSystem(nDof::S)

            K  = spzeros(nDof,nDof)
            M  = spzeros(nDof,nDof)
            C  = spzeros(nDof,nDof)
            dP = zeros(nDof)
            dU = zeros(nDof)
            UF = zeros(nDof)

            return new(nDof,K,M,C,dP,dU,UF)
        end
    end

    LinearAlgebraSystem{S<:Signed}(nDof::S) = LinearAlgebraSystem{S}(nDof)

    function linSolve(linalgsys::LinearAlgebraSystem)
        linalgsys.dU = linalgsys.dF\linalgsys.K
        return nothing
    end

    type Domain{S<:Integer}

        dim::S #// 维度
        dof::S #// 节点自由度

        name::String

        dofIndex::Dict{Vector{S},S} #// 节点自由度索引([节点标签, 节点局部自由度标签]=>全局自由度标签)
        IdofIndex::Dict{S,Vector{S}} #// 节点自由度索引(全局自由度标签=>[节点标签, 节点局部自由度标签])
        Nodes::Dict{S,Node}
        Refps::Dict{S,RefPoint}
        Mats::Dict{S,Material}
        Secs::Dict{S,Section}
        Eles::Dict{S,Element}
        Cons::Dict{Vector{S},Boundary}
        Loads::Dict{Vector{S},Boundary}

        nNode::S
        nEle::S
        nDof::S
        nActiveDof::S

        alength::R #// 单元平均长度
        time::Vector{R} #// 伪时间序列
        cstep::S #// 当前荷载步

        P::Vector{R}   # 节点荷载向量
        dP::Vector{R}  # 节点增量荷载向量
        U::Vector{R}   # 节点位移向量
        dU::Vector{R}  # 节点增量位移向量
        UF::Vector{R}  # 节点不平衡内力向量

        K::Matrix{R}  # 刚度矩阵
        M::Matrix{R}  # 质量矩阵
        C::Matrix{R}  # 阻尼矩阵

        linalgsys::LinearAlgebraSystem

        alpha::R
        beta::R

        function Domain(dim::S,dof::S,name::String)

            dofIndex = Dict{Vector{S},S}()
            IdofIndex = Dict{S,Vector{S}}()

            Nodes = Dict{S,Node}()
            Refps = Dict{S,RefPoint}()
            Mats = Dict{S,Material}()
            Secs = Dict{S,Section}()
            Eles = Dict{S,Element}()
            Cons = Dict{Vector{S},Boundary}()
            Loads = Dict{Vector{S},Boundary}()

            nNode = 0
            nEle  = 0
            nDof  = 0
            nActiveDof  = 0
            alength = 0.0 # 单元平均长度

            time = R[]
            cstep = 1

            P  = R[]   # 节点荷载向量
            dP = R[]  # 节点增量荷载向量
            U  = R[]   # 节点位移向量
            dU = R[]  # 节点增量位移向量
            UF = R[]  # 节点不平衡内力向量

            K = zeros(nDof,nDof)  # 刚度矩阵
            M = zeros(nDof,nDof)  # 质量矩阵
            C = zeros(nDof,nDof)  # 阻尼矩阵

            linalgsys = LinearAlgebraSystem(nDof)

            alpha = 0.0
            beta = 0.0

            domain = new( dim, dof, name, 
                        dofIndex, IdofIndex,
                        Nodes, Refps, Mats, Secs, Eles, 
                        Cons, Loads, 
                        nNode, nEle, nDof, nActiveDof, 
                        alength,
                        time, cstep, 
                        P , dP, U , dU, UF,
                        K, M, C,
                        linalgsys,
                        alpha, beta
                    )

            setCurrentDomain(domain)
            return domain
        end
    end

    Domain{S<:Integer}(dim::S,dof::S,name::String="model")=Domain{S}(dim,dof,name)

    function setCurrentDomain(domain::Domain)
        global currentDomain = domain
    end

    function addRefpoint(domain::Domain,tag::Int64,coord::Vector{Float64})

        rp = RefPoint(tag,coord)
        domain.Refps[rp.tag] = rp
        nothing

    end

    function addNode(domain::Domain,tag::Int64,coord::Vector{Float64},mass::Vector{Float64})
        # -- 添加节点 -- #
        node = Node(tag,coord,mass,domain.dim,domain.dof)
        domain.Nodes[node.tag] = node
        # domain.nDof += node.dof
        nothing
    end

    function addNode(domain::Domain,tag::Int64,coord::Vector{Float64})
        # -- 添加节点 -- #
        mass = zeros(domain.dof)
        node = Node(tag,coord,mass,domain.dim,domain.dof)
        domain.Nodes[node.tag] = node
        # domain.nDof += node.dof
        nothing
    end

    function addNode(tag::Int64,coord::Vector{Float64},mass::Vector{Float64})
        # -- 添加节点 -- #
        addNode(currentDomain,tag,coord,mass)
    end

    function addNode(tag::Int64,coord::Vector{Float64})
        # -- 添加节点 -- #
        addNode(currentDomain,tag,coord)
    end

    function setNodeMass(domain::Domain,nodeTag::Int64,mass::Vector{Float64})
        domain.Nodes[nodeTag].mass = mass
    end

    function setNodeMass(domain::Domain,nodeTag::Int64,mass::Float64)
        node = domain.Nodes[nodeTag]
        node.mass = mass.*ones(node.dof)
        if node.dof == 3
            node.mass[3] = eps(1.0)*mass
        elseif node.dof == 6
            node.mass[4:6] = eps(1.0)*mass
        end
    end

    function addMat(domain::Domain,matType::DataType,tag::Int64,arg...)
        mat = matType(tag,arg...)
        domain.Mats[mat.tag] = mat
    end

    function addMat(domain::Domain,mat::Material)
        domain.Mats[mat.tag] = mat
    end

    function addMat(matType::DataType,tag::Int64,arg...)
        addMat(currentDomain,matType,tag,arg...)
    end

    function addMat(mat::Material)
        addMat(currentDomain,mat)
    end

    function addSec(domain::Domain,secType::DataType,tag::Int64,arg...)
        sec = secType(tag,arg...)
        domain.Secs[sec.tag] = sec
    end

    function addSec(domain::Domain,sec::Section)
        domain.Secs[sec.tag] = sec
    end

    function addSec(secType::DataType,tag::Int64,arg...)
        addSec(currentDomain,secType,tag,arg...)
    end

    function addSec(sec::Section)
        addSec(currentDomain,sec)
    end

    function addEle(domain::Domain,eleType::DataType,tag::Int64,nodeTags::Vector{Int64},arg...)
        node = Node[]
        for nodeTag in nodeTags
            push!(node,domain.Nodes[nodeTag])
        end
        ele = eleType(tag,node,arg...)
        domain.alength += ele.len

        domain.Eles[ele.tag] = ele
    end

    function addEle(domain::Domain,ele::Element)
        domain.alength += ele.len
        domain.Eles[ele.tag] = ele
    end

    function addEle(eleType::DataType,tag::Int64,nodeTags::Vector{Int64},arg...)
        addEle(currentDomain,eleType,tag,nodeTags,arg...)
    end

    function addEle(ele::Element)
        addEle(currentDomain,ele)
    end

    function setTimeSeries(domain::Domain,time::Vector{Float64})
        domain.time = time
    end

    function setTimeSeries(domain::Domain,nstage::Int64,dts::Vector{Float64},nsteps::Vector{Int64})

        time = Float64[]
        for i in 1:nstage
            dt = dts[i]
            nstep = nsteps[i]
            t = linspace(dt,dt*nstep,nstep)
            append!(time,t)
        end
        domain.time = time
    end

    function setTimeSeries(time::Vector{Float64})
        setTimeSeries(currentDomain,time)
    end

    function setTimeSeries(domain::Domain,nstage::Int64,dts::Vector{Float64},nsteps::Vector{Int64})
        setTimeSeries(currentDomain,nstage,dts,nsteps)
    end

    function addFix(domain::Domain,node::Node,dofTag::Int64)
        # 添加节点位移约束
        cons = CNodeDisp(node,dofTag,0.0)
        domain.Cons[[node.tag,dofTag]] = cons
    end

    function addFix(domain::Domain,nodeTag::Int64,dofTag::Int64)
        # 添加节点位移约束
        node = domain.Nodes[nodeTag]
        cons = CNodeDisp(node,dofTag,0.0)
        domain.Cons[[node.tag,dofTag]] = cons
    end

    function addFix(node::Node,dofTag::Int64)
        # 添加节点位移约束
        addFix(currentDomain,node,dofTag)
    end

    function addFix(nodeTag::Int64,dofTag::Int64)
        # 添加节点位移约束
        addFix(currentDomain,nodeTag,dofTag)
    end

    function addCons(domain::Domain,node::Node,dofTag::Int64,dofValue::Float64)
        # 添加节点位移约束
        cons = CNodeDisp(node,dofTag,dofValue)
        domain.Cons[[node.tag,dofTag]] = cons
    end

    function addCons(domain::Domain,nodeTag::Int64,dofTag::Int64,dofValue::Float64)
        # 添加节点位移约束
        node = domain.Nodes[nodeTag]
        cons = CNodeDisp(node,dofTag,dofValue)
        domain.Cons[[node.tag,dofTag]] = cons
    end

    function addCons(node::Node,dofTag::Int64,dofValue::Float64)
        # 添加节点位移约束
        addCons(currentDomain,node,dofTag,dofValue)
    end

    function addCons(nodeTag::Int64,dofTag::Int64,dofValue::Float64)
        # 添加节点位移约束
        addCons(currentDomain,nodeTag,dofTag,dofValue)
    end

    function addLoad(domain::Domain,nodeTag::Int64,dofTag::Int64,dofValue::Float64)
        # 添加节点集中荷载
        load = CNodeForce(node,dofTag,dofValue)
        domain.Loads[[node.tag,dofTag]] = load
    end

    function addLoad(domain::Domain,node::Node,dofTag::Int64,dofValue::Float64)
        # 添加节点集中荷载
        node = domain.Nodes[nodeTag]
        load = CNodeForce(node,dofTag,dofValue)
        domain.Loads[[node.tag,dofTag]] = load
    end

    function addLoad(node::Node,dofTag::Int64,dofValue::Float64)
        # 添加节点集中荷载
        addLoad(currentDomain,node,dofTag,dofValue)
    end

    function addLoad(nodeTag::Int64,dofTag::Int64,dofValue::Float64)
        # 添加节点集中荷载
        addLoad(currentDomain,nodeTag,dofTag,dofValue)
    end

    function buildModel(domain::Domain)

        domain.alength /= domain.nEle
        domain.nNode = length(domain.Nodes)
        domain.nEle = length(domain.Eles)

        nodeDofNumbering(domain::Domain)

        # 初始化荷载向量 P, 节点变形向量 U, 节点不平衡力向量 UF
        domain.P = zeros(domain.nDof)
        domain.dP = zeros(domain.nDof)
        domain.U = zeros(domain.nDof)
        domain.dU = zeros(domain.nDof)
        domain.UF = zeros(domain.nDof)

        linalgsys = LinearAlgebraSystem(domain.nActiveDof)

        assembleSystemMatrix(domain)
        applyLoad(domain)
        applyCons(domain)

    end

    function nodeDofNumbering(domain::Domain)
        # # 节点自由度编号
        # > 非约束(Active)自由度与约束(Constraint)自由度编号分开
        # >> Active Dof range     --- 1:domain.nActiveDof
        # >> Constraint Dof range --- domain.nActiveDof+1:domain.nDof
        
        nodeTags = collect(keys(domain.Nodes)) #// 节点标签序列
        sort!(nodeTags)
        # domain.nDof = 0
        # for nodeTag in nodeTags
        #     node = domain.Nodes[nodeTag]
        #     domain.nDof += node.dof
        # end
        domain.nDof = domain.nNode*domain.dof

        nConsDof = length(domain.Cons) #// 受约束的自由度个数
        domain.nActiveDof = domain.nDof - nConsDof

        aOrder = 1 #// 未约束自由度编号计数
        cOrder = domain.nActiveDof + 1 #// 约束自由度编号计数

        for nodeTag in nodeTags
            node = domain.Nodes[nodeTag]
            for nodeDof in 1:node.dof
                if haskey(domain.Cons,[nodeTag,nodeDof])
                    domain.dofIndex[[nodeTag,nodeDof]] = cOrder
                    domain.IdofIndex[cOrder] = [nodeTag,nodeDof]
                    setDofOrder(node,nodeDof,cOrder)
                    cOrder += 1
                else
                    domain.dofIndex[[nodeTag,nodeDof]] = aOrder
                    domain.IdofIndex[aOrder] = [nodeTag,nodeDof]
                    setDofOrder(node,nodeDof,aOrder)
                    aOrder += 1
                end
            end
        end
    end

    function assembleSystemMatrix(domain::Domain)

        domain.K = zeros(domain.nDof,domain.nDof)
        domain.M = zeros(domain.nDof,domain.nDof)
        domain.C = zeros(domain.nDof,domain.nDof)

        for ele in values(domain.Eles)
            getStiff(ele)
            getStiffMatrix(ele)
            getMassMatrix(ele)
            getDofMatrix(ele)

            Rows = ele.dofMatrix
            Cols = (ele.dofMatrix).'
            for i in 1:length(ele.Kg)
                I,J = Rows[i],Cols[i]
                domain.K[I,J] += ele.Kg[i]
                domain.M[I,J] += ele.Mg[i]
                domain.C[I,J] += (domain.alpha*ele.Kg[i]+domain.beta*ele.Mg[i])
            end
        end
    end

    function assembleStiffMatrix(domain::Domain)

        domain.K = zeros(domain.nDof,domain.nDof)
        for ele in values(domain.Eles)
            getStiff(ele)
            getStiffMatrix(ele)
            # getDofMatrix(ele)
            Rows = ele.dofMatrix
            Cols = (ele.dofMatrix).'
            for i in 1:length(ele.Kg)
                domain.K[Rows[i],Cols[i]] += ele.Kg[i]
            end
        end
    end

    function assembleMassMatrix(domain::Domain;coupled::Bool=false)

        domain.M = zeros(domain.nDof,domain.nDof)
        for ele in values(domain.Eles)
            getMassMatrix(ele)
            # getDofMatrix(ele)
            Rows = ele.dofMatrix
            Cols = (ele.dofMatrix).'
            for i in 1:length(ele.Mg)
                domain.M[Rows[i],Cols[i]] += ele.Mg[i]
            end
        end
    end

    function assembleDampMatrix(domain::Domain;xi::Float64=0.05,order::Int64=2)
        domain.C = domain.alpha*domain.K + domain.beta*domain.M
    end

    function applyLoad(domain::Domain)
        # 施加节点荷载

        domain.P = zeros(domain.nDof)
        domain.dP = zeros(domain.nDof)
        for load in values(domain.Loads)
            dofOrder = domain.dofIndex[(load.node.tag,load.dofTag)]
            dofValue = getValue(load,domain.cstep)
            domain.P[dofOrder] += dofValue[1]
            domain.dP[dofOrder] += dofValue[2]
        end
    end

    function applyCons(domain::Domain)
        # '''施加支座约束'''

        for cons in values(domain.Cons)
            dofOrder = domain.dofIndex[[cons.node.tag,cons.dofTag]]
            dofValue = getValue(cons,domain.cstep)
            domain.U[dofOrder] += dofValue[1]
            domain.dU[dofOrder] += dofValue[2]
        end

        nActiveDof = domain.nActiveDof
        domain.linalgsys.K = sparse(domain.K[1:nActiveDof,1:nActiveDof])
        domain.linalgsys.M = sparse(domain.M[1:nActiveDof,1:nActiveDof])

        dPa = domain.dP[1:nActiveDof]
        dUc = domain.dU[nActiveDof+1:end]

        if norm(dUc) == 0.0
            domain.linalgsys.dP = dPa
        else
            Kac = domain.K[1:nActiveDof,nActiveDof+1:end]
            domain.linalgsys.dP = dPa - Kac*dUc
        end
    end

    function getSolution(domain::Domain)

        dUa = domain.linalgsys.dU[1:domain.nActiveDof]
        domain.dU[1:domain.nActiveDof] = dUa
        if domain.cstep <=1
            domain.U[1:domain.nActiveDof] = domain.dU[1:domain.nActiveDof]
        else
            domain.U[1:domain.nActiveDof] = domain.U[1:domain.nActiveDof]+domain.dU[1:domain.nActiveDof]
        end

        U = domain.U
        for node in values(domain.Nodes)
            setDisp(node,U)
        end

    end

    function getForce(domain::Domain)

        # 计算节点不平衡力向量(不更新刚度矩阵)

        domain.UF = zeros(domain.nDof)
        for ele in values(domain.Eles)
            setDeform(ele)
            getForce(ele)

            nDof = ele.nDof
            for i in 1:nDof
                dofOrder = ele.dofVector[i]
                domain.UF[dofOrder] += ele.nodeForces[i]
            end    
        end
        domain.linalgsys.UF = domain.UF[1:domain.nActiveDof]

        for i in domain.nActiveDof+1:domain.nDof
            nodeTag,dofTag = domain.IdofIndex[i]
            node = domain.Nodes[nodeTag]
            rf = -domain.UF[i]
            setReaction(node,dofTag,rf)
        end
    end

    function getForceStiff(domain::Domain)

        # 计算节点不平衡力向量同时更新刚度矩阵

        domain.K = zeros(domain.nDof,domain.nDof)
        domain.UF = zeros(domain.nDof)
        
        for ele in values(domain.ele)
            setDeform(ele)
            getForceStiff(ele)
            nDof = ele.nDof

            Rows = ele.dofMatrix
            Cols = (ele.dofMatrix).'

            for i in 1:length(ele.Kg)
                domain.K[Rows[i],Cols[i]] += ele.Kg[i]
            end

            for i in 1:nDof
                dofOrder = ele.dofVector[i]
                if dofOrder < domain.nActiveDof
                   domain.UF[dofOrder] += ele.nodeForces[i]
                end
            end
        end
        
        applyCons(domain)

    end

    function updateStiffHistPara(domain::Domain)
        # 采用修正牛顿法计算时,迭代收敛后更新刚度矩阵及历史变量
        for ni in values(domain.node)
            getUnbalancedForce(ni,domain.UF)
        end

        domain.K = zeros(domain.nDof,domain.nDof)
        for ei in values(domain.ele)
            getStiff(ei)
            tmp = 0.0
            for i in 1:ei.dof*ei.dof
                tmp = ei.KVals[i]
                if tmp != 0.0
                    domain.K[ei.KRows[i],ei.KCols[i]] += tmp
                end
            end
            updateHistoryPara(ei)
        end
        applyCons(domain)
    end

    function updateHistoryPara(domain::Domain)
        #  迭代收敛后更新历史变量
        for ni in values(domain.node)
            getUnbalancedForce(ni,domain.UF)
        end   

        for ei in values(domain.ele)
            updateHistoryPara(ei)
        end   
    end

    buildModel() = buildModel(currentDomain)
    nodeDofNumbering() = nodeDofNumbering(currentDomain)
    assembleSystemMatrix() = assembleSystemMatrix(currentDomain)
    applyLoad() = applyLoad(currentDomain)
    applyCons() = applyCons(currentDomain)
    getSolution() = getSolution(currentDomain)
    getForce() = getForce(currentDomain)
    getForceStiff() = getForceStiff(currentDomain)

    function readSAP2k(domain::Domain,fn::String)
        
        sf = open(fn,"r")
        ls = readlines(sf)
        titles = [
        "TABLE:  \"MATERIAL PROPERTIES 02 - BASIC MECHANICAL PROPERTIES\"\r\n",
        "TABLE:  \"FRAME SECTION PROPERTIES 01 - GENERAL\"\r\n",
        "TABLE:  \"JOINT COORDINATES\"\r\n",
        "TABLE:  \"CONNECTIVITY - FRAME\"\r\n",
        "TABLE:  \"FRAME SECTION ASSIGNMENTS\"\r\n",
        "TABLE:  \"ASSEMBLED JOINT MASSES\"\r\n",
        ]

        line = findin(ls,titles)

        matline, secline, jcline, cfline, fsline, jmline = line

        reg = r"\w+\=\-?\d+(\.\d+)?"

        matdict = Dict{String,Int64}()
        secdict = Dict{String,Int64}()
        
        i = matline+1
        lc = ls[i]
        while ismatch(reg,lc)
            execExpr(reg,lc)
            temp = match(r"Material=\w+",lc).match
            eval(parse(replace(temp,"=","=\"")*"\""))
            tag = i - matline
            matdict[Material] = tag
            matType = GeneralMaterial
            rho = UnitMass
            E0 = E1
            v = U12
            addMat(domain,matType,tag,rho,E0,v)
            i += 1
            lc = ls[i]
        end

        i = secline+1
        lc = ls[i]
        while ismatch(reg,lc)
            execExpr(reg,lc)
            temp = match(r"SectionName=\w+",lc).match
            eval(parse(replace(temp,"=","=\"")*"\""))
            temp = match(r"Material=\w+",lc).match
            eval(parse(replace(temp,"=","=\"")*"\""))
            tag = i - secline
            secdict[SectionName] = tag
            secType = GeneralSection
            A = Area
            I22 = I22
            I33 = I33
            mat = domain.Mats[matdict[Material]]
            addSec(domain,secType,tag,A,I33,I22,mat)
            i += 3
            lc = ls[i]
        end

        i = jcline+1
        lc = ls[i]
        mi = jmline+1
        mlc = ls[mi]
        while ismatch(reg,lc)
            execExpr(reg,lc)
            tag = Joint
            coord = float([XorR,Y,Z])
            execExpr(reg,mlc)
            mass = float([U1,U2,U3,R1,R2,R3])
            addNode(domain,tag,coord,mass)
            # println(tag,coord,mass)
            i += 1
            lc = ls[i]
            mi += 1
            mlc = ls[mi]
        end

        i = cfline+1
        lc = ls[i]
        si = fsline+1
        slc = ls[si]
        while ismatch(reg,lc)
            execExpr(reg,lc)
            tag = Frame
            nodeTags = [JointI,JointJ]
            eleType = GeneralLine
            temp = match(r"AnalSect=\w+",slc).match
            eval(parse(replace(temp,"=","=\"")*"\""))
            sec = domain.Secs[secdict[AnalSect]]
            addEle(domain,eleType,tag,nodeTags,sec)
            i += 1
            lc = ls[i]
            si += 1
            slc = ls[si]
        end
        nothing
    end

    function execExpr(reg::Regex,lc::String)

        expr = map(parse,matchall(reg,lc))
        exec = map(eval,expr)

        nothing
    end

    function saveModel(domain::Domain,fmt::String="gmsh")

        if fmt == "gmsh"
            writeGmsh(domain)
        else
            nothing
        end

    end

    saveModel(fmt::String="gmsh") = saveModel(currentDomain,fmt)

    function writeGmsh(domain::Domain)

        file = open("$(domain.name).msh","w")
        gmshHead = "\$MeshFormat\r\n2.8 0 8\r\n\$EndMeshFormat\r\n"
        write(file,gmshHead)
        write(file,"\$Nodes\r\n")
        write(file,"$(domain.nNode)\r\n")
        for node in values(domain.Nodes)
            i = node.tag
            x = node.coord[1]
            y = node.coord[2]
            z = 0.0
            if node.dim == 3
                z = node.coord[3]
            end
            write(file,"$i $x $y $z\r\n")
        end
        write(file,"\$EndNodes\r\n")
        
        nLine = 0
        lines = ""
        write(file,"\$Elements\r\n")
        for ele in values(domain.Eles)
            i = ele.tag
            if isa(ele,Line)
                nLine += 1
                n1 = ele.NodeI.tag
                n2 = ele.NodeJ.tag
                lines *= "$i 1 2 1 1 $n1 $n2\r\n"
            end
        end
        write(file,"$nLine\r\n")
        write(file,lines)
        write(file,"\$EndElements\r\n")
        
        close(file)

    end

    function setRayleighDampingFactor(domain::Domain,f1::Float64,f2::Float64,xi::Float64)

        alpha = 4.0*pi*xi*f1*f2/(f1+f2)
        beta = xi/pi/(f1+f2)
        domain.alpha = alpha
        domain.beta = beta
        nothing

    end

end