multiple dirichlet boundary conditions for vector valued functions. direct solver design.

src/JuliaFEM.jl

@@ -10,6 +10,10 @@ module JuliaFEM
 #@Logging.configure(level=DEBUG)
 #using Lexicon
 
+import Base: +, -, /, *, push!, convert, getindex, length, similar, call, vec, endof
+
+#importall Base
+
 """
 A very simple debugging macro. It prints debug message if environment variable
 JULIAFEM_DEBUG is found.
@@ -78,6 +82,7 @@ include("elasticity.jl")
 ### ASSEMBLY + SOLVE ###
 include("assembly.jl")
 include("solvers.jl")
+include("directsolver.jl") # parallel sparse direct solver for non-lniear problems
 
 # PRE AND POSTPROCESS
 include("xdmf.jl")

src/directsolver.jl

@@ -0,0 +1,114 @@
+# This file is a part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
+
+## Direct solver
+
+type DirectSolver <: Solver
+    field_problems :: Vector{FieldProblem}
+    boundary_problems :: Vector{BoundaryProblem}
+    nonlinear_problem :: Bool
+    max_iterations :: Int64
+    tol :: Float64
+end
+
+function push!(solver::DirectSolver, problem::FieldProblem)
+    push!(solver.field_problems, problem)
+end
+
+function push!(solver::DirectSolver, problem::BoundaryProblem)
+    push!(solver.boundary_problems, problem)
+end
+
+""" Default initializer. """
+function DirectSolver()
+    DirectSolver([], [], true, 10, 1.0e-6)
+end
+
+""" Call solver to solve a set of problems. """
+function call(solver::DirectSolver, time::Number=0.0)
+    @assert length(solver.field_problems) == 1
+    @assert length(solver.boundary_problems) == 1
+    @assert solver.nonlinear_problem == true
+
+    problem1 = solver.field_problems[1]
+    problem2 = solver.boundary_problems[1]
+
+    x = zeros(3)
+    dx = zeros(3)
+    dims = nothing
+
+    for iter=1:solver.max_iterations
+        tic()
+        info("Starting iteration $iter")
+        assembly1 = Assembly()
+        assemble!(assembly1, problem1, time)
+        assembly2 = Assembly()
+        assemble!(assembly2, problem2, time)
+
+        A1 = sparse(assembly1.stiffness_matrix)
+        dims = size(A1)
+        b1 = sparse(assembly1.force_vector, dims[1], 1)
+        A2 = sparse(assembly2.stiffness_matrix, dims[1], dims[2])
+        b2 = sparse(assembly2.force_vector, dims[1], 1)
+    
+        # create a saddle point problem
+        A = [A1 A2; A2' zeros(A2)]
+        b = [b1; b2]
+
+        if length(b) != length(x)
+            info("iter $iter: resizing solution vector")
+            resize!(x, length(b))
+            resize!(dx, length(b))
+            fill!(x, 0.0)
+            fill!(dx, 0.0)
+        end
+
+        # solve problem, update solution vector
+        nz = unique(rowvals(A))  # take only non-zero rows
+        dx[nz] = lufact(A[nz,nz]) \ full(b[nz])
+        x += dx
+
+        # get "problem-wise" solution vectors
+        x1 = x[1:dims[1]]
+        x2 = x[dims[1]+1:end]
+
+        # update field for elements in problem 1
+        for equation in get_equations(problem1)
+            element = get_element(equation)
+            field_name = get_unknown_field_name(problem1)
+            gdofs = get_gdofs(problem1, equation)
+            local_sol = vec(full(x1[gdofs]))
+            eqsize = size(equation)
+            if eqsize[1] != 1
+                local_sol = reshape(local_sol, eqsize)
+            end
+            #info("problem1: pushing to $field_name")
+            push!(element[field_name], time => local_sol)
+        end
+
+        # update field for elements in problem 2 (Dirichlet boundary)
+        for equation in get_equations(problem2)
+            element = get_element(equation)
+            field_name = "reaction force" #get_unknown_field_name(problem2)
+            gdofs = get_gdofs(problem2, equation)
+            local_sol = vec(full(x1[gdofs]))
+            eqsize = size(equation)
+            if eqsize[1] != 1
+                local_sol = reshape(local_sol, eqsize)
+            end
+            #info("problem2: pushing to $field_name")
+            push!(element[field_name], time => local_sol)
+        end
+        
+        if norm(dx[1:dims[1]]) < solver.tol
+            return (iter, true)
+        end
+
+        info("Iteration took $(toq()) seconds")
+    end
+
+    info("Warning: did not coverge in $(solver.max_iterations) iterations!")
+    return (solver.max_iterations, false)
+
+end
+

src/dirichlet.jl

@@ -3,11 +3,7 @@
 
 # Dirichlet boundary conditions in weak form
 
-abstract DirichletEquation <: Equation
-
-function get_unknown_field_name(equation::DirichletEquation)
-    return "reaction force"
-end
+abstract DirichletEquation <: BoundaryEquation
 
 ### Dirichlet problem + equations
 
@@ -15,8 +11,6 @@ type DirichletProblem <: BoundaryProblem
     unknown_field_name :: ASCIIString
     unknown_field_dimension :: Int
     equations :: Vector{DirichletEquation}
-#   element_mapping :: Dict{DataType, DataType}
-    field_value :: Function
 end
 
 """ Initialize new Dirichlet boundary condition.
@@ -24,27 +18,14 @@ end
 Parameters
 ----------
 dimension
-    dimension of unknown field
-field_value
-    boundary function
+    dimension of unknown field (scalar, vector, ...)
 
 Examples
 --------
 
-Create u(X) = 0.0 boundary condition for three-dimensional elasticity problem:
-
->>> u(X) = [0.0, 0.0, 0.0]
->>> bc = DirichletProblem(3, u)
 """
-function DirichletProblem(dimension::Int=1, field_value::Function=(X)->[0.0,0.0,0.0], equations=[])
-#    element_mapping = nothing
-#    if dimension == 1
-#        element_mapping = Dict(
-#            Seg2 => DBC2D2
-#            )
-#    end
-    DirichletProblem("reaction force", dimension, equations, field_value)
-#    DirichletProblem("reaction force", dimension, [], element_mapping, field_value)
+function DirichletProblem(unknown_field_name::ASCIIString, dimension::Int=1)
+    DirichletProblem(unknown_field_name, dimension, [])
 end
 
 """ Dirichlet boundary condition element for 2 node line segment """
@@ -57,31 +38,45 @@ function Base.size(equation::DBC2D2)
     return (1, 2)
 end
 
-#function DBC2D2(element::Seg2)
 function Base.convert(::Type{DirichletEquation}, element::Seg2)
     integration_points = line3()
-    haskey(element, "reaction force") || (element["reaction force"] = zeros(1, 2))
+    haskey(element, "reaction force") || (element["reaction force"] = 0.0 => zeros(2))
     DBC2D2(element, integration_points)
 end
 
 
 function assemble!(assembly::Assembly, equation::DirichletEquation, time::Number=0.0, problem=nothing)
-    gdofs = get_gdofs(equation)
-#   info("gdofs = $gdofs")
+    isa(problem, Void) && error("Dicihlet boundary condition needs problem defined")
+    field_dim = problem.unknown_field_dimension
+    field_name = problem.unknown_field_name
     element = get_element(equation)
+    gdofs = get_gdofs(element, field_dim)
     basis = get_basis(element)
     detJ = det(basis)
     for ip in get_integration_points(equation)
         w = ip.weight * detJ(ip)
         N = basis(ip, time)
-        add!(assembly.stiffness_matrix, gdofs, gdofs, w*N'*N)
-#       info("added $(w*N'*N)")
-#        if !isa(problem, Void)
-#            X = basis("geometry", ip, time)
-#            u = problem.field_value(X)[1:length(gdofs)]
-#            add!(assembly.force_vector, gdofs, w*N'*u)
-#        end
+        A = w*N'*N
+
+        if haskey(element, field_name)
+            # add all dimensions at once
+            for i=1:field_dim
+                g = element(field_name, ip, time)
+                ldofs = gdofs[i:field_dim:end]
+                add!(assembly.stiffness_matrix, ldofs, ldofs, A)
+                add!(assembly.force_vector, ldofs, w*g*N')
+            end
+        end
+  
+        for i=1:field_dim
+            if haskey(element, field_name*" $i")
+                # add single component
+                g = element(field_name*" $i", ip, time)
+                ldofs = gdofs[i:field_dim:end]
+                add!(assembly.stiffness_matrix, ldofs, ldofs, A)
+                add!(assembly.force_vector, ldofs, w*g*N')
+            end
+        end
     end
-#   info("assembly done")
 end
 

src/elasticity.jl

@@ -3,7 +3,7 @@
 
 # Elasticity problems
 
-abstract ElasticityProblem <: Problem
+abstract ElasticityProblem <: FieldProblem
 abstract ElasticityEquation <: Equation
 
 function get_unknown_field_name(equation::ElasticityEquation)

src/elements.jl

@@ -67,6 +67,10 @@ function Base.getindex(element::Element, field_name)
     return element.fields[field_name]
 end
 
+function get_integration_points(element)
+    return get_default_integration_points(element)
+end
+
 """Add new Field to element.
 
 Examples
@@ -78,6 +82,9 @@ Examples
 function Base.setindex!(element::Element, data, name::ASCIIString)
     element.fields[name] = Field(data)
 end
+function Base.setindex!(element::Element, field::Field, name::ASCIIString)
+    element.fields[name] = field
+end
 function Base.setindex!(element::Element, data::Tuple, name::ASCIIString)
     element.fields[name] = Field(data...)
 end
@@ -143,6 +150,39 @@ function call(gradu::GradientFunctionSpace, field_name, xi::Union{Vector, Integr
     gradu.basis(geometry, field, xi, Val{:grad})
 end
 
+typealias VecOrIP Union{Vector, IntegrationPoint}
+
+function Base.call(element::Element, field_name::ASCIIString, xi::VecOrIP, time::Number)
+    return element.basis(element[field_name](time), xi)
+end
+
+function Base.call(element::Element, field_name::ASCIIString, xi::VecOrIP, time::Number, ::Type{Val{:grad}})
+    return element.basis(element["geometry"](time), element[field_name](time), xi, Val{:grad})
+end
+
+function Base.call(element::Element, field_name::ASCIIString, xi::VecOrIP)
+    return element.basis(element[field_name], xi)
+end
+
+function Base.call(element::Element, field_name::ASCIIString, xi::VecOrIP, ::Type{Val{:grad}})
+    return element.basis(element["geometry"], element[field_name], xi, Val{:grad})
+end
+
+function Base.call(element::Element, field_name::ASCIIString, time::Number)
+    return element[field_name](time)
+end
+
+function Base.call(element::Element, xi::VecOrIP)
+    element.basis(xi)
+end
+
+function Base.call(element::Element, xi::VecOrIP, ::Type{Val{:grad}})
+    element.basis(element["geometry"], xi, Val{:grad})
+end
+
+function Base.call(element::Element, field_name::ASCIIString)
+    return element[field_name]
+end
 
 # on-line functions to get api more easy to use, ip -> xi.ip
 #call(u::FunctionSpace, ip::IntegrationPoint, t::Number=Inf) = call(u, ip.xi, t)
@@ -176,12 +216,19 @@ function LinAlg.det(u::FunctionSpace, xi::Vector, time::Number=0.0)
     m, n = size(J)
     return m == n ? det(J) : norm(J)
 end
+
 function LinAlg.det(u::FunctionSpace, ip::IntegrationPoint, time::Number=0.0)
     LinAlg.det(u, ip.xi, time)
 end
+
 function LinAlg.det(u::FunctionSpace)
     return (args...) -> det(u, args...)
 end
+
+function LinAlg.det(element::Element)
+    return det(get_basis(element))
+end
+
 #Base.(:+)(u::FunctionSpace, v::FunctionSpace) = (args...) -> u(args...) + v(args...)
 #Base.(:-)(u::FunctionSpace, v::FunctionSpace) = (args...) -> u(args...) - v(args...)
 #Base.(:+)(u::GradientFunctionSpace, v::GradientFunctionSpace) = (args...) -> u(args...) + v(args...)

src/equations.jl

@@ -4,6 +4,8 @@
 # Functions to handle element level things -- integration, assembly, ...
 
 abstract Equation
+abstract FieldEquation <: Equation
+abstract BoundaryEquation <: Equation
 
 type Assembly
     mass_matrix :: SparseMatrixIJV

src/fields2.jl

@@ -12,6 +12,7 @@ abstract Variable <: AbstractField
 abstract TimeVariant <: AbstractField
 abstract TimeInvariant <: AbstractField
 
+
 type Field{A<:Union{Discrete,Continuous}, B<:Union{Constant,Variable}, C<:Union{TimeVariant,TimeInvariant}}
     data
 end
@@ -81,6 +82,10 @@ typealias TimeVariantField   Union{DCTV, DVTV, CCTV, CVTV}
 
 ### Convenient functions to create fields
 
+#function Base.convert(::Type{Field}, data)
+#    return Field(data)
+#end
+
 function Field(data)
     return DCTI(data)
 end
@@ -101,6 +106,20 @@ function Base.convert{T}(::Type{DCTV}, data::Pair{Float64, Vector{T}}...)
     return DCTV([Increment{Vector{T}}(d[1], d[2]) for d in data])
 end
 
+function Field(func::Function)
+    if method_exists(func, Tuple{})
+        return CCTI(func)
+    elseif method_exists(func, Tuple{Float64})
+        return CCTV(func)
+    elseif method_exists(func, Tuple{Vector})
+        return CVTI(func)
+    elseif method_exists(func, Tuple{Vector, Number})
+        return CVTV(func)
+    else
+        error("no proper definition found for function: check methods.")
+    end
+end
+
 function CVTI(basis::Function, dbasis::Function)
     return CVTI(Basis(basis, dbasis))
 end
@@ -203,28 +222,55 @@ function Base.convert(::Type{Basis}, field::CVTI)
     return field.data
 end
 
+function Base.call(field::CCTV, time::Number)
+    return field.data(time)
+end
+
 ### Interpolation 
 
+""" Interpolate time-invariant field in time direction. """
 function Base.call(field::DVTI, time::Float64)
-    # interpolating time-invariant field in time direction -> no effect
     return field
 end
-
 function Base.call(field::DCTI, time::Float64)
-    # interpolating time-invariant field in time direction -> no effect
     return field
 end
+function Base.call(field::CVTI, time::Float64)
+    return field.data()
+end
+function Base.call(field::CCTI, time::Float64)
+    return field.data()
+end
 
+""" Interpolate time-variant field in time direction. """
+function Base.call(field::DCTV, time::Float64)
+    for i=reverse(1:length(field))
+        if isapprox(field[i].time, time)
+            return DCTI(field[i].data)
+        end
+    end
+    info(field.data)
+    info(time)
+    error("interpolate DCTV: not implemented yet")
+end
+
+function Base.call(field::DVTV, time::Float64, time_extrapolation::Symbol=:linear)
+    for i=reverse(1:length(field))
+        if isapprox(field[i].time, time)
+            return DVTI(field[i].data)
+        end
+    end
+    info(field.data)
+    info(time)
+    error("interpolate DVTV: not implemented yet")
+end
+
+""" Interpolate constant field in spatial dimension. """
 function Base.call(basis::CVTI, field::DCTI, xi::Vector)
-    # try to interpolate constant value -> no effect
-    return field
+    return field.data
 end
 
-#function Base.call(basis::Basis, field::DCTI, xi::Vector)
-    # calling constant field with basis -> no effect
-#    return field
-#end
-
+""" Interpolate variable field in spatial dimension. """
 function Base.call(basis::CVTI, values::DVTI, xi::Vector)
     N = basis(xi)
     return sum([N[i]*values[i] for i=1:length(N)])
@@ -244,27 +290,8 @@ function Base.call(basis::CVTI, geometry::DVTI, values::DVTI, xi::Vector, ::Type
     return length(gradf) == 1 ? gradf[1] : gradf
 end
 
-function Base.call(field::DCTV, time::Float64)
-    for i in length(field)
-        if isapprox(field[i].time, time)
-            return DCTI(field[i].data)
-        end
-    end
-    error("interpolate DCTV: not implemented yet")
-end
-
-function Base.call(field::DVTV, time::Float64, time_extrapolation::Symbol=:linear)
-#    info("length of field DVTV: $(length(field))")
-    for i=reverse(1:length(field))
-        res = isapprox(field[i].time, time)
-        #info("isapprox $(field[i].time) to $time ? $res")
-        if isapprox(field[i].time, time)
-            return DVTI(field[i].data)
-        end
-    end
-    info(field.data)
-    info(time)
-    error("interpolate DVTV: not implemented yet")
+function Base.call(basis::CVTI, xi::Vector, time::Number)
+    call(basis, xi)
 end
 
 ### FIELDSET ###

src/solvers.jl

@@ -9,7 +9,7 @@ abstract Solver
 Solve field equations for single element with some dofs fixed. This can be used
 to test nonlinear element formulations.
 """
-function solve!(equation::Equation, free_dofs::Vector{Int}, time::Number; max_iterations::Int=10, tolerance::Float64=1.0e-12, dump_matrices::Bool=false)
+function solve!(equation::Equation, free_dofs::Vector{Int}, time::Number; max_iterations::Int=10, tolerance::Float64=1.0e-12, dump_matrices::Bool=false, callback=nothing)
     unknown_field_name = get_unknown_field_name(equation)
     element = get_element(equation)
     x0 = element[unknown_field_name](0.0)
@@ -31,6 +31,9 @@ function solve!(equation::Equation, free_dofs::Vector{Int}, time::Number; max_it
         data = eqsize[1] != 1 ? reshape(x, eqsize) : x
         push!(element[unknown_field_name], time => data)
         norm(dx) < tolerance && return
+        if !isa(callback, Void)
+            callback(x)
+        end
     end
     error("Did not converge in $max_iterations iterations")
 end
@@ -41,7 +44,7 @@ to test nonlinear element formulations. Dirichlet boundary is assumed to be homo
 and degrees of freedom are eliminated. So if boundary condition is known in nodal
 points and everything is zero this should be quite good.
 """
-function solve!(problem::Problem, free_dofs::Vector{Int}, time::Float64; max_iterations::Int=10, tolerance::Float64=1.0e-12, dump_matrices::Bool=false)
+function solve!(problem::Problem, free_dofs::Vector{Int}, time::Float64; max_iterations::Int=10, tolerance::Float64=1.0e-12, dump_matrices::Bool=false, callback=nothing)
     info("start solver")
     assembly = Assembly()
 #   x = zeros(ga.ndofs)
@@ -66,6 +69,9 @@ function solve!(problem::Problem, free_dofs::Vector{Int}, time::Float64; max_ite
         dx[free_dofs] = lufact(A[free_dofs,free_dofs]) \ full(b)[free_dofs]
         info("Difference in solution norm: $(norm(dx))")
         x += dx
+        if !(isa(callback, Void))
+            callback(x)
+        end
         for equation in get_equations(problem)
             element = get_element(equation)
             gdofs = get_gdofs(equation)
@@ -81,11 +87,6 @@ function solve!(problem::Problem, free_dofs::Vector{Int}, time::Float64; max_ite
     error("Did not converge in $max_iterations iterations")
 end
 
-""" Add new problem to solver. """
-function add_problem!(solver::Solver, problem::Problem)
-    push!(solver.problems, problem)
-end
-
 function Base.push!(solver::Solver, problem::Problem)
     push!(solver.problems, problem)
 end
@@ -126,9 +127,10 @@ function call(solver::SimpleSolver, time::Number=0.0)
 
 #   info("Creating sparse matrices")
     A1 = sparse(assembly1.stiffness_matrix)
-    b1 = sparse(assembly1.force_vector, size(A1, 1), 1)
-    A2 = sparse(assembly2.stiffness_matrix)
-    b2 = sparse(assembly2.force_vector, size(A2, 1), 1)
+    dims = size(A1)
+    b1 = sparse(assembly1.force_vector, dims[1], 1)
+    A2 = sparse(assembly2.stiffness_matrix, dims[1], dims[2])
+    b2 = sparse(assembly2.force_vector, dims[1], 1)
     
     # create a saddle point problem
     A = [A1 A2; A2' zeros(A2)]
@@ -149,16 +151,28 @@ function call(solver::SimpleSolver, time::Number=0.0)
         field_name = get_unknown_field_name(problem1)
         gdofs = get_gdofs(problem1, equation)
         local_sol = vec(full(x1[gdofs]))
+        eqsize = size(equation)
+        if eqsize[1] != 1
+            local_sol = reshape(local_sol, eqsize)
+        end
+        #info("problem1: pushing to $field_name")
         push!(element[field_name], time => local_sol)
     end
 
     # update field for elements in problem 2 (Dirichlet boundary)
     for equation in get_equations(problem2)
         element = get_element(equation)
-        field_name = get_unknown_field_name(problem2)
+        field_name = "reaction force" #get_unknown_field_name(problem2)
         gdofs = get_gdofs(problem2, equation)
         local_sol = vec(full(x1[gdofs]))
+        eqsize = size(equation)
+        if eqsize[1] != 1
+            local_sol = reshape(local_sol, eqsize)
+        end
+        #info("problem2: pushing to $field_name")
         push!(element[field_name], time => local_sol)
     end
+
+    return norm(x1)
 end