elasticity solver example updated

src/abaqus_reader_old.jl

@@ -0,0 +1,140 @@
+# This file is a part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
+
+global handlers = Dict()
+
+"""
+Register new handler for parser
+"""
+function add_handler(section, function_name)
+  handlers[section] = function_name
+end
+
+function create_or_get(model, key)
+    if !(key in keys(model))
+        model[key] = Dict()
+    end
+    return model[key]
+end
+
+function parse_header(header_line)
+    args = map(s -> strip(s), split(header_line, ","))
+    args[1] = strip(args[1], '*')
+    d = Dict("section" => args[1], "options" => Dict())
+    options = d["options"]
+    for k in args[2:end]
+        args2 = split(k, "=")
+        options[args2[1]] = args2[2]
+    end
+    return d
+end
+
+function parse_node_section(model, header, data)
+    nodes = create_or_get(model, "nodes")
+    for line in split(data, "\n")
+        m = matchall(r"[-0-9.]+", line)
+        id = parse(Int, m[1])
+        coords = float(m[2:end])
+        nodes[id] = coords
+    end
+end
+
+function parse_element_section(model, header, data)
+    info("Parsing elements")
+    eldims = Dict(
+        "C3D10" => 10,
+        "C3D4" => 4,
+        "S3" => 3)
+    eltype = header["options"]["TYPE"]
+    if !(eltype in keys(eldims))
+        throw("Element $eltype dimension information missing")
+    end
+    eldim = eldims[eltype]
+    test_match = matchall(r"[0-9]+", "234, 242")
+    m = matchall(r"[0-9]+", data)   
+    m = map((s) -> parse(Int, s), m)
+    elements = create_or_get(model, "elements")
+    m = reshape(m, eldim+1, round(Int, length(m)/(eldim+1)))
+    nel = size(m)[2]
+    info("$nel elements found")
+    for i=1:nel
+        elements[m[1,i]] = m[2:end,i]
+    end
+    if "ELSET" in keys(header["options"])
+        elsets = create_or_get(model, "elsets")
+        elset_name = header["options"]["ELSET"]
+        info("Creating ELSET $elset_name")
+        elsets[elset_name] = Int64[]
+        for i=1:nel
+            push!(elsets[elset_name], m[1,i])
+        end
+    end
+end
+
+function parse_elset_section(model, header, data)
+    elset_name = header["options"]["ELSET"]
+    info("Creating element set $elset_name")
+    m = matchall(r"[0-9]+", data)
+    element_ids = map((s) -> parse(Int, s), m)
+    elsets = create_or_get(model, "elsets")
+    elsets[elset_name] = Int64[]
+    for j in element_ids
+        push!(elsets[elset_name], j)
+    end
+end
+
+function parse_nodeset_section(model, header, data)
+    nset_name = header["options"]["NSET"]
+    info("Creating node set $nset_name")
+    m = matchall(r"[0-9]+", data)
+    node_ids = map((s) -> parse(Int, s), m)
+    nsets = create_or_get(model, "nsets")
+    nsets[nset_name] = Int64[]
+    for j in node_ids
+        push!(nsets[nset_name], j)
+    end
+end
+
+function parse_abaqus(fid::IOStream)
+    model = Dict()
+    section = nothing
+    header = nothing
+    data = ASCIIString[]
+    info("Registered handlers: $(keys(handlers))")
+
+    function process_section(section)
+        if section == nothing
+            return
+        end
+        if !(section in keys(handlers))
+            info("Don't know what to do with data in section $section")
+            info("Skipping $(length(data)) bytes of unknown data")
+            return
+        end
+        joined = join(data, "")
+        handlers[section](model, header, strip(joined))
+        empty!(data)
+    end
+
+    line_idx = 0
+    for line in eachline(fid)
+        if startswith(line, "**")
+            continue
+        end
+        if startswith(line, "*")
+            process_section(section)
+            header = parse_header(line)
+            section = header["section"]
+            continue
+        end
+        push!(data, line)
+    end 
+    process_section(section)
+    return model
+end
+
+# add handlers
+add_handler("NODE", parse_node_section)
+add_handler("ELEMENT", parse_element_section)
+add_handler("NSET", parse_nodeset_section)
+add_handler("ELSET", parse_elset_section)

src/directsolver.jl

@@ -25,12 +25,28 @@ function DirectSolver()
     DirectSolver([], [], false, true, 10, 1.0e-6)
 end
 
+function tic(timing, what::ASCIIString)
+    timing[what * " start"] = time()
+end
+
+function toc(timing, what::ASCIIString)
+    timing[what * " finish"] = time()
+end
+
+function time_elapsed(timing, what::ASCIIString)
+    return timing[what * " finish"] - timing[what * " start"]
+end
+
 """ Call solver to solve a set of problems. """
 function call(solver::DirectSolver, time::Number=0.0)
     #@assert length(solver.field_problems) == 1
     info("# of field problems: $(length(solver.field_problems))")
     info("# of boundary problems: $(length(solver.boundary_problems))")
     @assert solver.nonlinear_problem == true
+    
+    timing = Dict{ASCIIString, Float64}()
+    tic(timing, "solver")
+    tic(timing, "initialization")
 
     # check that all problems are "same kind"
     field_name = get_unknown_field_name(solver.field_problems[1])
@@ -73,26 +89,33 @@ function call(solver::DirectSolver, time::Number=0.0)
         end
     end
 
+    toc(timing, "initialization")
+
     dim = 0
 
     for iter=1:solver.max_iterations
-        tic()
         info("Starting iteration $iter")
+        tic(timing, "non-linear iteration")
 
         mapper = solver.parallel ? pmap : map
 
         # 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)
+        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")
 
+        tic(timing, "create sparse matrices")
         # create sparse matrices and saddle point problem
         K = sparse(field_assembly.stiffness_matrix)
         dim = size(K, 1)
@@ -101,14 +124,18 @@ function call(solver::DirectSolver, time::Number=0.0)
         g = sparse(boundary_assembly.force_vector, dim, 1)
         A = [K C'; C spzeros(dim, dim)]
         b = [r; g]
+        toc(timing, "create sparse matrices")
 
+        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]))
 
+        tic(timing, "update element data")
         # update elements in field problems
         for field_problem in solver.field_problems
             for element in get_elements(field_problem)
@@ -130,10 +157,22 @@ function call(solver::DirectSolver, time::Number=0.0)
                 last(element["reaction force"]).data = local_sol  # <-- replaced
             end
         end
+        toc(timing, "update element data")
+        toc(timing, "non-linear iteration")
 
-        info("Non-linear iteration took $(toq()) seconds")
+        if true
+            info("timing info for non-linear iteration:")
+            info("boundary assembly      : ", time_elapsed(timing, "boundary assembly"))
+            info("field assembly         : ", time_elapsed(timing, "field assembly"))
+            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"))
+            info("non-linear iteration   : ", time_elapsed(timing, "non-linear iteration"))
+        end
 
         if norm(sol[1:dim]) < solver.tol
+            toc(timing, "solver")
+            info("solver finished in ", time_elapsed(timing, "solver"), " seconds.")
             return (iter, true)
         end
 

src/elasticity.jl

@@ -5,12 +5,6 @@
 
 abstract ElasticityProblem <: AbstractProblem
 
-abstract PlaneStressElasticityProblem <: ElasticityProblem
-
-function PlaneStressElasticityProblem(dim::Int=2, elements=[])
-    return Problem{PlaneStressElasticityProblem}(dim, elements)
-end
-
 function get_unknown_field_name{P<:ElasticityProblem}(::Type{P})
     return "displacement"
 end
@@ -19,6 +13,16 @@ function get_unknown_field_type{P<:ElasticityProblem}(::Type{P})
     return Vector{Float64}
 end
 
+function ElasticityProblem(dim::Int=3, elements=[])
+    return Problem{PlaneStressElasticityProblem}(dim, elements)
+end
+
+abstract PlaneStressElasticityProblem <: ElasticityProblem
+
+function PlaneStressElasticityProblem(dim::Int=2, elements=[])
+    return Problem{PlaneStressElasticityProblem}(dim, elements)
+end
+
 """ Elasticity equations.
 
 Formulation
@@ -49,25 +53,31 @@ https://en.wikipedia.org/wiki/Plane_stress
 https://en.wikipedia.org/wiki/Hooke's_law
 
 """
-function get_residual_vector{EL<:CG}(problem::Problem{PlaneStressElasticityProblem}, element::Element{EL}, ip::IntegrationPoint, time::Number; variation=nothing)
+function get_residual_vector{P<:ElasticityProblem}(problem::Problem{P}, element::Element, ip::IntegrationPoint, time::Number; variation=nothing)
 
     basis = element(ip, time)
-    dbasis = element(ip, time, Val{:grad})
 
     u = element("displacement", ip, time, variation)
-    gradu = element("displacement", ip, time, Val{:grad}, variation)
-    F = I + gradu # deformation gradient
+
+    r = zeros(Float64, problem.dim, length(element))
 
     # internal forces
-    young = element("youngs modulus", ip, time)
-    poisson = element("poissons ratio", ip, time)
-    mu = young/(2*(1+poisson))
-    lambda = young*poisson/((1+poisson)*(1-2*poisson))
-    lambda = 2*lambda*mu/(lambda + 2*mu)  # <- correction for 2d
-    E = 1/2*(F'*F - I)  # strain
-    S = lambda*trace(E)*I + 2*mu*E
+    if haskey(element, "youngs modulus") && haskey(element, "poissons ratio")
+        dbasis = element(ip, time, Val{:grad})
+        gradu = element("displacement", ip, time, Val{:grad}, variation)
+        F = I + gradu # deformation gradient
 
-    r = F*S*dbasis
+        young = element("youngs modulus", ip, time)
+        poisson = element("poissons ratio", ip, time)
+        mu = young/(2*(1+poisson))
+        lambda = young*poisson/((1+poisson)*(1-2*poisson))
+        if P == PlaneStressElasticityProblem
+            lambda = 2*lambda*mu/(lambda + 2*mu)  # <- correction for 2d problems
+        end
+        E = 1/2*(F'*F - I)  # strain
+        S = lambda*trace(E)*I + 2*mu*E
+        r += F*S*dbasis
+    end
 
     # external forces - volume load
     if haskey(element, "displacement load")
@@ -75,92 +85,12 @@ function get_residual_vector{EL<:CG}(problem::Problem{PlaneStressElasticityProbl
         r -= b*basis
     end
 
-    return vec(r)
-end
-
-""" Surface load for plane stress model. """
-function get_residual_vector(problem::Problem{PlaneStressElasticityProblem}, element::Element{Seg2}, ip::IntegrationPoint, time::Number; variation=nothing)
-
-    u = element("displacement", ip, time, variation)
-    r = zeros(problem.dim, length(element))
-
+    # external forces - surface traction force
     if haskey(element, "displacement traction force")
         T = element("displacement traction force", ip, time)
-        r -= T*element(ip, time)
+        r -= T*basis
     end
 
     return vec(r)
 end
 
-
-#=
-
-### 3d continuum elasticity ###
-
-type ContinuumElasticityProblem <: ElasticityProblem
-    unknown_field_name :: ASCIIString
-    unknown_field_dimension :: Int
-    equations :: Vector{ElasticityEquation}
-end
-
-function ContinuumElasticityProblem(equations=[])
-    return PlaneStressElasticityProblem("displacement", 3, equations)
-end
-
-### Equations ###
-
-""" 4-node plane stress element. """
-type C3D10 <: ContinuumElasticityEquation
-    element :: Quad4
-    integration_points :: Vector{IntegrationPoint}
-end
-
-function Base.size(equation::CPS4)
-    return (2, 4)
-end
-
-function Base.convert(::Type{PlaneStressElasticityEquation}, element::Quad4)
-    integration_points = get_integration_points(element)
-    if !haskey(element, "displacement")
-        element["displacement"] = 0.0 => [zeros(2) for i=1:4]
-    end
-    CPS4(element, integration_points)
-end
-
-""" Boundary element for plane stress problem for surface loads. """
-type CPS2 <: PlaneStressElasticityEquation
-    element :: Seg2
-    integration_points :: Vector{IntegrationPoint}
-end
-
-function Base.size(equation::CPS2)
-    return (2, 2)
-end
-
-function Base.convert(::Type{PlaneStressElasticityEquation}, element::Seg2)
-    integration_points = get_integration_points(element)
-    if !haskey(element, "displacement")
-        element["displacement"] = 0.0 => [zeros(2) for i=1:2]
-    end
-    CPS2(element, integration_points)
-end
-
-function get_residual_vector(equation::CPS2, ip::IntegrationPoint, time::Number; variation=nothing)
-
-    element = get_element(equation)
-    basis = get_basis(element)
-
-    u = basis("displacement", ip, time, variation)
-    r = zeros(size(equation))
-
-    if haskey(element, "displacement traction force")
-        T = basis("displacement traction force", ip, time)
-        r -= T*basis(ip, time)
-    end
-
-    return vec(r)
-end
-
-=#
-
-

src/integrate.jl

@@ -2,17 +2,9 @@
 # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
 
 # Let's drop here all integration schemes and some defaults for different element types
+# maybe parse from txt file ..?
 
-function get_integration_points(::Type{Quad4})
-    [
-        IntegrationPoint(1.0/sqrt(3.0)*[-1, -1], 1.0),
-        IntegrationPoint(1.0/sqrt(3.0)*[ 1, -1], 1.0),
-        IntegrationPoint(1.0/sqrt(3.0)*[ 1,  1], 1.0),
-        IntegrationPoint(1.0/sqrt(3.0)*[-1,  1], 1.0)
-    ]
-end
-
-typealias LineElement Union{Seg2, Seg3}
+### 1d elements
 
 function get_integration_points(::Type{Seg2}, ::Type{Val{1}})
     [
@@ -27,7 +19,7 @@ function get_integration_points(::Type{Seg2}, ::Type{Val{2}})
     ]
 end
 
-function get_integration_points(element::LineElement, ::Type{Val{3}})
+function get_integration_points(::Type{Seg3}, ::Type{Val{3}})
     [
         IntegrationPoint([0.0], 8/9),
         IntegrationPoint([-sqrt(3/5)], 5/9),
@@ -35,7 +27,7 @@ function get_integration_points(element::LineElement, ::Type{Val{3}})
     ]
 end
 
-function get_integration_points(element::LineElement, ::Type{Val{4}})
+function get_integration_points(::Type{Seg3}, ::Type{Val{4}})
     [
         IntegrationPoint([+sqrt(3/7 - 2/7*sqrt(6/5))], (18+sqrt(30))/36)
         IntegrationPoint([-sqrt(3/7 - 2/7*sqrt(6/5))], (18+sqrt(30))/36)
@@ -44,7 +36,7 @@ function get_integration_points(element::LineElement, ::Type{Val{4}})
     ]
 end
 
-function get_integration_points(element::LineElement, ::Type{Val{5}})
+function get_integration_points(::Type{Seg3}, ::Type{Val{5}})
     [
         IntegrationPoint([-1/3*sqrt(5 + 2*sqrt(10/7))], (322-13*sqrt(70))/900),
         IntegrationPoint([-1/3*sqrt(5 - 2*sqrt(10/7))], (322+13*sqrt(70))/900),
@@ -58,24 +50,56 @@ function get_integration_points(::Type{Seg2})
     return get_integration_points(Seg2, Val{2})
 end
 
-function get_integration_points(element::Seg3)
-    return get_integration_points(element, Val{3})
+function get_integration_points(::Type{Seg3})
+    return get_integration_points(Seg3, Val{3})
 end
 
-### 3D elements
+### 2d elements
 
-function get_integration_points(element::Tet10, ::Type{Val{4}})
+function get_integration_points(::Type{Tri3})
+    # http://libmesh.github.io/doxygen/quadrature__gauss__2D_8C_source.html
+    [
+        IntegrationPoint([1.0/3.0, 1.0/3.0], 0.5)
+    ]
+end
+
+function get_integration_points(::Type{Tri6})
+    # http://libmesh.github.io/doxygen/quadrature__gauss__2D_8C_source.html
+    [
+        IntegrationPoint([2.0/3.0, 1.0/6.0], 1.0/6.0)
+        IntegrationPoint([1.0/6.0, 2.0/3.0], 1.0/6.0)
+        IntegrationPoint([1.0/6.0, 1.0/6.0], 1.0/6.0)
+    ]
+end
+
+function get_integration_points(::Type{Quad4})
+    [
+        IntegrationPoint(1.0/sqrt(3.0)*[-1, -1], 1.0),
+        IntegrationPoint(1.0/sqrt(3.0)*[ 1, -1], 1.0),
+        IntegrationPoint(1.0/sqrt(3.0)*[ 1,  1], 1.0),
+        IntegrationPoint(1.0/sqrt(3.0)*[-1,  1], 1.0)
+    ]
+end
+
+### 3d elements
+
+function get_integration_points(::Type{Tet4})
+    # http://libmesh.github.io/doxygen/quadrature__gauss__3D_8C_source.html
+    [
+        IntegrationPoint([0.25, 0.25, 0.25], 1.0/6.0)
+    ]
+end
+
+function get_integration_points(::Type{Tet10})
+    # http://libmesh.github.io/doxygen/quadrature__gauss__3D_8C_source.html
     a = .585410196624969
     b = .138196601125011
     w = .041666666666667
-    integration_points = [
+    [
         IntegrationPoint([a, b, b], w),
         IntegrationPoint([b, a, b], w),
         IntegrationPoint([b, b, a], w),
-        IntegrationPoint([b, b, b], w)]
-end
-
-function get_integration_points(element::Tet10)
-    return get_integration_points(element, Val{4})
+        IntegrationPoint([b, b, b], w)
+    ]
 end
 

src/lagrange.jl

@@ -64,6 +64,11 @@ end
      0.0 0.0 1.0],
     (xi) -> [1.0, xi[1], xi[2]])
 
+@create_lagrange_element(Tri6, "6 node quadratic triangle element",
+    [0.0 1.0 0.0 0.5 0.5 0.0
+     0.0 0.0 1.0 0.0 0.5 0.5],
+    (xi) -> [1.0, xi[1], xi[2], xi[1]^2, xi[2]^2, xi[1]*xi[2]])
+
 @create_lagrange_element(Quad4, "4 node bilinear quadrangle element",
     [-1.0  1.0 1.0 -1.0
      -1.0 -1.0 1.0  1.0],
@@ -71,6 +76,12 @@ end
 
 # 3d Lagrange elements
 
+@create_lagrange_element(Tet4, "4 node tetrahedron",
+    [0.0 1.0 0.0 0.0
+     0.0 0.0 1.0 0.0
+     0.0 0.0 0.0 1.0],
+    (xi) -> [1.0, xi[1], xi[2], xi[3]])
+
 @create_lagrange_element(Tet10, "10 node quadratic tetrahedron",
     [0.0 1.0 0.0 0.0 0.5 0.5 0.0 0.0 0.5 0.0
      0.0 0.0 1.0 0.0 0.0 0.5 0.5 0.0 0.0 0.5