JuliaFEM.jl/src/problems.jl

# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md

abstract AbstractProblem
abstract FieldProblem <: AbstractProblem
abstract BoundaryProblem <: AbstractProblem
abstract MixedProblem <: AbstractProblem

"""
General linearized problem to solve
    K*u +   C1'*la  = f
    C2*u +  D*la    = g
"""
type Assembly
    # for field assembly
    M :: SparseMatrixCOO  # mass matrix
    K :: SparseMatrixCOO  # stiffness matrix
    f :: SparseMatrixCOO  # force vector
    # for boundary assembly
    C1 :: SparseMatrixCOO
    C2 :: SparseMatrixCOO
    D :: SparseMatrixCOO 
    g :: SparseMatrixCOO

    u :: Vector{Float64}  # solution vector u
    u_prev :: Vector{Float64}  # previous solution vector u
    u_norm_change :: Real  # change of norm in u

    la :: Vector{Float64}  # solution vector la
    la_prev :: Vector{Float64}  # previous solution vector u
    la_norm_change :: Real # change of norm in la

    prehooks :: Vector{Tuple{Symbol,Any,Any}}  # assign possible prehooks before assembly
    posthooks :: Vector{Tuple{Symbol,Any,Any}}  # assign possible posthooks after assembly
    changed :: Bool  # flag to control is reassembly needed
end

function Assembly()
    return Assembly(
        SparseMatrixCOO(),
        SparseMatrixCOO(),
        SparseMatrixCOO(),
        SparseMatrixCOO(),
        SparseMatrixCOO(),
        SparseMatrixCOO(),
        SparseMatrixCOO(),
        [], [], Inf,
        [], [], Inf,
        [], [], true)
end

function Base.empty!(assembly::Assembly)
    empty!(assembly.M)
    empty!(assembly.K)
    empty!(assembly.f)
    empty!(assembly.C1)
    empty!(assembly.C2)
    empty!(assembly.D)
    empty!(assembly.g)
    assembly.changed = true
end

type Problem{P<:AbstractProblem}
    name :: ASCIIString              # descriptive name for problem
    dimension :: Int                 # degrees of freedom per node
    parent_field_name :: ASCIIString # (optional) name of parent field e.g. "displacement"
    elements :: Vector{Element}
    assembly :: Assembly
    properties :: P
end

""" Construct a new field problem.

Examples
--------
Create vector-valued (dim=3) elasticity problem:

julia> prob = Problem(Elasticity, "this is my problem", 3)

"""
function Problem{P<:FieldProblem}(::Type{P}, name, dimension, elements=[])
    Problem{P}(name, dimension, "none", elements, Assembly(), P())
end

""" Construct a new boundary problem.

Examples
--------
Create Dirichlet boundary problem for vector-valued (dim=3) elasticity problem.

julia> bc1 = Problem(Dirichlet, "support", 3, "displacement")

"""
function Problem{P<:BoundaryProblem}(::Type{P}, name, dimension, parent_field_name, elements=[])
    Problem{P}(name, dimension, parent_field_name, elements, Assembly(), P())
end

function get_formulation_type{P<:FieldProblem}(problem::Problem{P})
    return :total
end

function get_formulation_type{P<:BoundaryProblem}(problem::Problem{P})
    return :total
end

function get_assembly(problem)
    return problem.assembly
end

""" Initialize unknown field ready for nonlinear iterations, i.e.,
    take last known value and set it as a initial quess for next
    time increment.
"""
function initialize!(problem::Problem, time::Real)
    field_name = get_unknown_field_name(problem)
    field_dim = get_unknown_field_dimension(problem)
    for element in get_elements(problem)
        gdofs = get_gdofs(element, problem)
        if haskey(element, field_name)
            # if field is found, copy last known solution to new time as initial guess
            if !isapprox(last(element[field_name]).time, time)
                last_data = copy(last(element[field_name]).data)
                push!(element[field_name], time => last_data)
            end
        else # if field not found at all, initialize new zero field.
            data = Vector{Float64}[zeros(field_dim) for i in 1:length(element)]
            element[field_name] = (time => data)
        end
    end
end

""" Update problem solution vector for assembly. """
function update_assembly!(problem, u, la)

    assembly = get_assembly(problem)

    # resize & fill with zeros vectors if length mismatch with current solution
    if length(u) != length(assembly.u)
        resize!(assembly.u, length(u))
        fill!(assembly.u, 0.0)
    end
    if length(la) != length(assembly.la)
        resize!(assembly.la, length(la))
        fill!(assembly.la, 0.0)
    end

    # copy current solutions to previous ones and add/replace new solution
    assembly.u_prev = copy(assembly.u)
    assembly.la_prev = copy(assembly.la)
    if get_formulation_type(problem) == :incremental
        assembly.u += u
        assembly.la += la
    else
        assembly.u = u
        assembly.la = la
    end

    # calculate change of norm
    assembly.u_norm_change = norm(assembly.u - assembly.u_prev)
    assembly.la_norm_change = norm(assembly.la - assembly.la_prev)
    return assembly.u_norm_change, assembly.la_norm_change
end

""" Update solutions to elements.

Notes
-----
This assumes that element is properly initialized so that last known field data
is from current time. For boundary problems solution is updated from lambda vector
and for field problems from actual solution vector.
"""
function update_elements!(problem, u, la)
    field_name = get_unknown_field_name(problem)
    field_dim = get_unknown_field_dimension(problem)
    nnodes = round(Int, length(u)/field_dim)

    solution = nothing
    if is_field_problem(problem)
        solution = reshape(u, field_dim, nnodes)
    elseif is_boundary_problem(problem)
        solution = reshape(la, field_dim, nnodes)
    else
        error("update_elements!(): unknown problem type $(typeof(problem))")
    end

    for element in get_elements(problem)
        connectivity = get_connectivity(element) # node ids
        local_sol = Vector{Float64}[solution[:, node_id] for node_id in connectivity]
        last(element[field_name]).data = local_sol
    end
end

#=
function add_postprocessor!(problem::Union{FieldProblem, BoundaryProblem}, postprocessor_name::Symbol, args...; kwargs...)
    push!(problem.postprocessors, (postprocessor_name, args, kwargs))
end

function add_preprocessor!(problem::Union{FieldProblem, BoundaryProblem}, preprocessor_name::Symbol, args...; kwargs...)
    push!(problem.preprocessors, (preprocessor_name, args, kwargs))
end

=#

function get_elements(problem)
    return problem.elements
end

""" Return the dimension of the unknown field of this problem. """
function get_unknown_field_dimension(problem::Problem)
    return problem.dimension
end

""" Return the name of the unknown field of this problem. """
function get_unknown_field_name{P}(problem::Problem{P})
    return get_unknown_field_name(P)
end

""" Return the name of the parent field of this (boundary) problem. """
function get_parent_field_name{P<:BoundaryProblem}(problem::Problem{P})
    return problem.parent_field_name
end

function push!(problem::Problem, element)
    push!(problem.elements, element)
end