substructuring, fixed tests, possibility to save to integration points

src/abaqus_reader_old.jl

@@ -44,7 +44,8 @@ function parse_element_section(model, header, data)
     eldims = Dict(
         "C3D10" => 10,
         "C3D4" => 4,
-        "S3" => 3)
+        "S3" => 3,
+        "STRI65" => 6)
     eltype = header["options"]["TYPE"]
     if !(eltype in keys(eldims))
         throw("Element $eltype dimension information missing")

src/assembly.jl

@@ -3,6 +3,16 @@
 
 # Functions to handle global assembly of problem
 
+type CAssembly
+    interior_dofs :: Vector{Int}
+    boundary_dofs :: Vector{Int}
+    F :: Factorization
+    Kc :: SparseMatrixCSC
+    fc :: SparseMatrixCSC
+    Ki :: SparseMatrixCSC
+    fi :: SparseMatrixCSC
+end
+
 function assemble!(assembly::Assembly, problem::AllProblems, time::Float64, empty_assembly::Bool=true)
     if empty_assembly
         empty!(assembly)
@@ -20,6 +30,73 @@ function assemble(problem::AllProblems, time::Float64)
     return assembly
 end
 
+""" Return condensed system. """
+function assemble(problem::FieldProblem, time::Float64, boundary_dofs::Vector{Int})
+    assembly = Assembly()
+    for element in get_elements(problem)
+        assemble!(assembly, problem, element, time)
+    end
+    return condensate(assembly, boundary_dofs)
+end
+
+function condensate(assembly::Assembly, boundary_dofs_::Vector{Int})
+    K = sparse(assembly.stiffness_matrix)
+    all_dofs = unique(assembly.stiffness_matrix.I)
+    boundary_dofs = intersect(all_dofs, boundary_dofs_)
+    interior_dofs = setdiff(all_dofs, boundary_dofs_)
+    dim = size(K, 1)
+    f = sparse(assembly.force_vector, dim, 1)
+
+    # check that matrix is symmetric
+    asdf = maximum(abs(1/2*(K + K') - K))
+    if asdf > 1.0e-6
+        info(full(K))
+        error("asdf $asdf > 1.0e-6")
+    end
+    
+    K = 1/2*(K + K')
+    F::Factorization = cholfact(K[interior_dofs, interior_dofs])
+
+#   info("condensation: all dofs: ", all_dofs)
+#   info("condensation: interior dofs: ", interior_dofs)
+#   info("condensation: boundary dofs: ", boundary_dofs)
+#   info("manually condensated")
+#   Kman = K[boundary_dofs, boundary_dofs] - K[boundary_dofs,interior_dofs] * inv(full(K[interior_dofs, interior_dofs])) * K[interior_dofs, boundary_dofs]
+#   info("\n$(full(Kman))")
+    #info("K = \n$(full(K))")
+    #Ki = K[interior_dofs, boundary_dofs]
+    Ki = K[interior_dofs, boundary_dofs]
+    fi = f[interior_dofs]
+#   info("condensated using factorization")
+#   LL = K[boundary_dofs, boundary_dofs] - K[boundary_dofs, interior_dofs] * (K[interior_dofs, interior_dofs] \ K[interior_dofs, boundary_dofs])
+#   info(LL)
+
+    Ks = F \ Ki
+    Fs = F \ fi
+
+    dim = size(K, 1)
+    Kc = spzeros(dim, dim)
+    fc = spzeros(dim, 1)
+    Kc[boundary_dofs, boundary_dofs] = K[boundary_dofs, boundary_dofs] - Ki' * Ks
+    fc[boundary_dofs] = f[boundary_dofs] - Ki' * Fs
+
+
+    return CAssembly(interior_dofs, boundary_dofs, F, Kc, fc, Ki, fi)
+end
+
+function reconstruct!(ca::CAssembly, x::SparseMatrixCSC)
+#   info("size of la = ", size(la))
+#   info("size of ca.Ki = ", size(ca.Ki))
+#   info("size of ca.fi = ", size(ca.fi))
+#   info("size of la[ca.interior_dofs] = ", size(la[ca.interior_dofs]))
+#   info("interior dofs: $(ca.interior_dofs)")
+#   info("boundary dofs: $(ca.boundary_dofs)")
+#   info("ca.fi = $(ca.fi')")
+#   info("sol1 = ", full(ca.F \ ca.fi)')
+#   info("sol2 = ", full(ca.F \ (ca.Ki*x[ca.boundary_dofs]))')
+    x[ca.interior_dofs] += ca.F \ (ca.fi - ca.Ki*x[ca.boundary_dofs])
+end
+
 function Base.(:+)(ass1::Assembly, ass2::Assembly)
     mass_matrix = ass1.mass_matrix + ass2.mass_matrix
     stiffness_matrix = ass1.stiffness_matrix + ass2.stiffness_matrix

src/directsolver.jl

@@ -99,48 +99,93 @@ function call(solver::DirectSolver, time::Number=0.0)
 
         mapper = solver.parallel ? pmap : map
 
+        info("Assembling problems.")
         # assemble boundary problems
         tic(timing, "boundary assembly")
         boundary_assembly = sum(mapper((p)->assemble(p, time), solver.boundary_problems))
         boundary_dofs = unique(boundary_assembly.stiffness_matrix.I)
+#       boundary_dofs = collect(range(1, 12))
+        info("# of interface dofs: $(length(boundary_dofs))")
+        #info("dofs = $boundary_dofs")
         toc(timing, "boundary assembly")
 
-        # assemble field problems
-        # in principle if we want to static condensation we need to pass boundary dofs
-        # to field problems in order to know which dofs are interior dofs and can be
-        # condensated.
-        tic(timing, "field assembly")
-        field_assembly = sum(mapper((p)->assemble(p, time), solver.field_problems))
-        field_dofs = unique(field_assembly.stiffness_matrix.I)
-        info("# of dofs: $(length(field_dofs)), # of interface dofs: $(length(boundary_dofs))")
-        toc(timing, "field assembly")
+        #static_condensation = false
 
+        # assemble field problems
+        dim = 0
+        assemblies = []
+        for (i, problem) in enumerate(solver.field_problems)
+            tic(timing, "field assembly")
+            field_assembly = assemble(problem, time)
+            #info("full assembly body $i")
+            #info(round(full(field_assembly.stiffness_matrix), 3))
+            toc(timing, "field assembly")
+            field_dofs = unique(field_assembly.stiffness_matrix.I)
+            #info("# of dofs in problem $i: $(length(field_dofs))")
+            dim = maximum([dim, maximum(field_dofs)])
+            #tic(timing, "condensate")
+            #cfield_assembly = condensate(field_assembly, boundary_dofs)
+            #toc(timing, "condensate")
+            #push!(assemblies, cfield_assembly)
+            push!(assemblies, field_assembly)
+        end
+
+        #info("dim = $dim")
+
+        #info("assembly done")
         tic(timing, "create sparse matrices")
         # create sparse matrices and saddle point problem
-        K = sparse(field_assembly.stiffness_matrix)
-        dim = size(K, 1)
-        r = sparse(field_assembly.force_vector, dim, 1)
+        #K = sparse(field_assembly.stiffness_matrix)
+        #dim = size(K, 1)
+        #r = sparse(field_assembly.force_vector, dim, 1)
+
+        K = spzeros(dim, dim)
+        r = spzeros(dim, 1)
+
+        for (i, assembly) in enumerate(assemblies)
+            #info("body $i")
+            #info(round(full(assembly.Kc), 3))
+            #resize!(assembly.stiffness_matrix, dim, dim)
+            #resize!(assembly.force_vector, dim, 1)
+            K += sparse(assembly.stiffness_matrix, dim, dim)
+            r += sparse(assembly.force_vector, dim, 1)
+        end
+        #info(round(full(K), 3))
+
         C = sparse(boundary_assembly.stiffness_matrix, dim, dim)
         g = sparse(boundary_assembly.force_vector, dim, 1)
         A = [K C'; C spzeros(dim, dim)]
         b = [r; g]
         toc(timing, "create sparse matrices")
+        #info("problem size = ", size(A))
 
+        info("Solving system")
         tic(timing, "solution of system")
         # solve increment for linearized problem
         nz = unique(rowvals(A))  # take only non-zero rows
         sol = zeros(b)
-        sol[nz] = lufact(A[nz,nz]) \ full(b[nz])
-        info("solved. length of solution vector = $(length(sol))")
-        toc(timing, "solution of system")
-        #info(full(sol[nz]))
+        sol[nz] = A[nz,nz] \ full(b[nz])
+        #info("solution vector before reconstruction")
+        #info(full(sol)')
+        la = sol[dim+1:end]
+        #for assembly in assemblies
+        #    reconstruct!(assembly, sol)
+        #end
+        la = vec(full(la))
+        sol = vec(full(sol))
+        #info("la = ", la')
+        #info("sol = ", sol')
 
+        info("solved. solution norm: $(norm(sol[1:dim]))")
+        toc(timing, "solution of system")
+
+        info("Updating element data")
         tic(timing, "update element data")
         # update elements in field problems
         for field_problem in solver.field_problems
             for element in get_elements(field_problem)
                 gdofs = get_gdofs(element, field_dim)
-                local_sol = vec(full(sol[gdofs]))  # incremental data for element
+                local_sol = sol[gdofs]  # incremental data for element
                 local_sol = reshape(local_sol, field_dim, length(element))
                 local_sol = Vector{Float64}[local_sol[:,i] for i=1:length(element)]
                 last(element[field_name]).data += local_sol  # <-- added
@@ -150,8 +195,8 @@ function call(solver::DirectSolver, time::Number=0.0)
         # update elements in boundary problems
         for boundary_problem in solver.boundary_problems
             for element in get_elements(boundary_problem)
-                gdofs = get_gdofs(element, field_dim) + dim
-                local_sol = vec(full(sol[gdofs]))
+                gdofs = get_gdofs(element, field_dim)
+                local_sol = la[gdofs]
                 local_sol = reshape(local_sol, field_dim, length(element))
                 local_sol = Vector{Float64}[local_sol[:,i] for i=1:length(element)]
                 last(element["reaction force"]).data = local_sol  # <-- replaced
@@ -164,6 +209,7 @@ function call(solver::DirectSolver, time::Number=0.0)
             info("timing info for non-linear iteration:")
             info("boundary assembly      : ", time_elapsed(timing, "boundary assembly"))
             info("field assembly         : ", time_elapsed(timing, "field assembly"))
+#           info("condensate             : ", time_elapsed(timing, "condensate"))
             info("create sparse matrices : ", time_elapsed(timing, "create sparse matrices"))
             info("solution of system     : ", time_elapsed(timing, "solution of system"))
             info("update element data    : ", time_elapsed(timing, "update element data"))

src/elasticity.jl

@@ -76,6 +76,12 @@ function get_residual_vector{P<:ElasticityProblem}(problem::Problem{P}, element:
         end
         E = 1/2*(F'*F - I)  # strain
         S = lambda*trace(E)*I + 2*mu*E
+
+        J = det(element, ip, time)
+        T = J^-1*F*S*F'
+        #ip["cauchy stress"] = T
+        ip["gl strain"] = E
+
         r += F*S*dbasis
     end
 

src/elements.jl

@@ -5,7 +5,6 @@ abstract AbstractElement
 
 type Element{E<:AbstractElement}
     connectivity :: Vector{Int}
-#   integration_points :: Vector{IntegrationPoint}
     fields :: Dict{ASCIIString, Field}
 end
 
@@ -15,10 +14,23 @@ function convert{E}(::Type{Element{E}}, connectivity::Vector{Int})
 end
 
 function get_integration_points{E}(element::Element{E})
-#   return element.integration_points
     return get_integration_points(E)
 end
 
+function update_gauss_fields!(element::Element, data::Vector{IntegrationPoint}, time::Real)
+    if haskey(element, "integration points")
+        # push or update
+        if !isapprox(last(element["integration points"]).time, time)
+            push!(element["integration points"], time => data)
+        else
+            last(element["integration points"]).data = data
+        end
+    else
+        # create
+        element["integration points"] = Field(time => data)
+    end
+end
+
 """
 Test routine for element. If this passes, element interface is properly
 defined.
@@ -46,7 +58,7 @@ function test_element(element_type)
 
     info("Initializing element")
     try
-        element = element_type(collect(1:n))
+        element = Element{element_type}(collect(1:n))
     catch
         error("""
         Unable to create element with default constructor define function
@@ -55,22 +67,20 @@ function test_element(element_type)
     end
 
     # try to interpolate some scalar field
-    element["field1"] = Field(collect(1:n))
+    element["field1"] = range(1, n)
     # TODO: how to parametrize this?
-    element["geometry"] = Field(Vector[[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]])
+    element["geometry"] = Vector{Float64}[[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]]
 
     # evaluate basis functions at middle point of element
-    basis = get_basis(element)
-    dbasis = grad(basis)
     mid = zeros(dim)
-    val1 = basis(mid, 0.0)
+    val1 = element(mid, 0.0)
     info("basis at $mid: $val1")
-    val2 = basis("field1", mid, 0.0)
+    val2 = element("field1", mid, 0.0)
     info("field val at $mid: $val2")
-    val3 = dbasis(mid, 0.0)
+    val3 = element(mid, 0.0, Val{:grad})
     info("derivative of basis at $mid:\n$val3")
-    val4 = dbasis("field1", mid, 0.0)
-    info("field val at $mid: $val4")
+    #val4 = element("field1", mid, Val{:grad})
+    #info("field val at $mid: $val4")
 
     info("Element $element_type passed tests.")
 end
@@ -122,11 +132,16 @@ function call(element::Element, field_name::ASCIIString, xi::VecOrIP, time::Numb
 end
 
 function call(element::Element, field_name::ASCIIString, xi::VecOrIP)
-    return element.basis(element[field_name], xi)
+    field = element[field_name]
+    basis = get_basis(element)
+    return basis(element[field_name], xi)
 end
 
 function call(element::Element, field_name::ASCIIString, xi::VecOrIP, ::Type{Val{:grad}})
-    return element.basis(element["geometry"], element[field_name], xi, Val{:grad})
+    field = element[field_name]
+    geom = element["geometry"]
+    basis = get_basis(element)
+    return basis(geom, field, xi, Val{:grad})
 end
 
 function call(element::Element, field_name::ASCIIString, time::Number)
@@ -137,6 +152,10 @@ function get_basis{E}(element::Element{E}, ip::IntegrationPoint)
     return get_basis(E, ip.xi)
 end
 
+function get_basis{E}(::Type{Element{E}}, xi::Vector{Float64})
+    return get_basis(E, xi)
+end
+
 function get_basis{E}(element::Element{E}, xi::Vector{Float64})
     return get_basis(E, xi)
 end
@@ -152,9 +171,15 @@ function get_basis{E}(element::Element{E})
     return basis
 end
 
+function call{E}(element::Element{E}, xi::VecOrIP, ::Type{Val{:grad}})
+    basis = get_basis(element)
+    geom = element["geometry"]
+    return basis(geom, xi, Val{:grad})
+end
+
 function call{E}(element::Element{E}, xi::VecOrIP, time::Float64, ::Type{Val{:grad}})
     basis = get_basis(element)
-    return basis(element["geometry"], xi, Val{:grad})
+    return basis(element["geometry"](time), xi, Val{:grad})
 end
 
 function call(element::Element, field_name::ASCIIString)

src/equations.jl

@@ -95,7 +95,7 @@ function assemble!(assembly::Assembly, problem::Problem, element::Element, time:
     end
 
     # 2. energy form -- user has defined potential energy W -> min!
-    if has_potential_energy(problem, element)
+    if has_potential_energy(problem, element) && haskey(element, unknown_field_name)
         field = element[unknown_field_name](time)
 
         """ Wrapper for potential energy for ForwardDiff. """
@@ -122,7 +122,7 @@ function assemble!(assembly::Assembly, problem::Problem, element::Element, time:
     end
 
     # 3. virtual work -- user has defined some residual r = p - f = 0
-    if has_residual_vector(problem, element)
+    if has_residual_vector(problem, element) && haskey(element, unknown_field_name)
 
         field = DVTI(last(element[unknown_field_name]).data)
 
@@ -130,17 +130,24 @@ function assemble!(assembly::Assembly, problem::Problem, element::Element, time:
         function calc_R(data::Vector)
             R = zeros(length(data))
             df = similar(field, data)
+            gauss_fields = IntegrationPoint[]
             # integrate residual vector
             for ip in get_integration_points(element)
                 s = ip.weight*det(element, ip, time)
                 dr = get_residual_vector(problem, element, ip, time; variation=df)
                 R += s*dr
+                if ip.changed
+                    push!(gauss_fields, ip)
+                end
             end
             # external loads -- if any nodal loads is defined, decrease from residual
             if haskey(element, "$unknown_field_name nodal load")
                 R -= vec(element["$unknown_field_name nodal load"](time))
             end
             #info("return = $R")
+            if length(gauss_fields) != 0
+                update_gauss_fields!(element, gauss_fields, time)
+            end
             return R
         end
 

src/problems.jl

@@ -15,6 +15,8 @@ type BoundaryProblem{T<:AbstractProblem}
     elements :: Vector{Element}
 end
 
+typealias FieldProblem Problem
+
 typealias AllProblems Union{Problem, BoundaryProblem}
 
 function get_elements(problem::AllProblems)

src/types.jl

@@ -21,34 +21,48 @@ fields
 type IntegrationPoint
     xi :: Vector
     weight :: Float64
-    fields :: FieldSet
+    fields :: Dict{ASCIIString, Field}
+    changed :: Bool
 end
 
 function IntegrationPoint(xi, weight)
-    return IntegrationPoint(xi, weight, FieldSet())
+    return IntegrationPoint(xi, weight, FieldSet(), false)
 end
 
-function Base.convert(::Type{Number}, ip::IntegrationPoint)
+function setindex!{T<:ForwardDiff.ForwardDiffNumber}(ip::IntegrationPoint, data::Array{T,2}, field_name::ASCIIString)
+    data = ForwardDiff.get_value(data)
+    setindex!(ip, data, field_name)
+end
+function setindex!(ip::IntegrationPoint, data, field_name)
+    ip.fields[field_name] = Field(data)
+    ip.changed = true
+end
+
+function getindex(ip::IntegrationPoint, field_name::ASCIIString)
+    ip.fields[field_name]
+end
+
+function convert(::Type{Number}, ip::IntegrationPoint)
     return ip.xi
 end
 
-function Base.call(field::CVTI, ip::IntegrationPoint)
+function call(field::CVTI, ip::IntegrationPoint)
     return call(field, ip.xi)
 end
 
-function Base.call(basis::CVTI, field::DCTI, ip::IntegrationPoint)
+function call(basis::CVTI, field::DCTI, ip::IntegrationPoint)
     call(basis, field, ip.xi)
 end
 
-function Base.call(basis::CVTI, field::DVTI, ip::IntegrationPoint, ::Type{Val{:grad}})
+function call(basis::CVTI, field::DVTI, ip::IntegrationPoint, ::Type{Val{:grad}})
     call(basis, field, ip.xi, Val{:grad})
 end
 
-function Base.call(basis::CVTI, field::DVTI, ip::IntegrationPoint)
+function call(basis::CVTI, field::DVTI, ip::IntegrationPoint)
     call(basis, field, ip.xi)
 end
 
-function Base.call(basis::CVTI, geometry::DVTI, field::DVTI, ip::IntegrationPoint, ::Type{Val{:grad}})
+function call(basis::CVTI, geometry::DVTI, field::Union{DCTI, DVTI}, ip::IntegrationPoint, ::Type{Val{:grad}})
     call(basis, geometry, field, ip.xi, Val{:grad})
 end