Merge 215ee13 into 4d337b7

REQUIRE

@@ -2,3 +2,5 @@ julia 1.0
 FEMBase
 ForwardDiff
 Tensors
+NLsolve
+Einsum

examples/continuum.jl

@@ -0,0 +1,138 @@
+# This file is a part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/Materials.jl/blob/master/LICENSE
+
+mutable struct Continuum3D <: FieldProblem
+    material_model :: Symbol
+end
+
+Continuum3D() = Continuum3D(:IdealPlastic)
+FEMBase.get_unknown_field_name(::Continuum3D) = "displacement"
+
+function FEMBase.assemble_elements!(problem::Problem{Continuum3D},
+                                    assembly::Assembly,
+                                    elements::Vector{Element{Hex8}},
+                                    time::Float64)
+
+
+    for element in elements
+        for ip in get_integration_points(element)
+            material = ip("material", time)
+            preprocess_increment!(material, element, ip, time)
+        end
+    end
+    bi = BasisInfo(Hex8)
+
+    dim = 3
+    nnodes = 8
+    ndofs = dim*nnodes
+
+    BL = zeros(6, ndofs)
+    Km = zeros(ndofs, ndofs)
+    f_int = zeros(ndofs)
+    f_ext = zeros(ndofs)
+    D = zeros(6, 6)
+    S = zeros(6)
+
+    dtime = 0.05
+    # super dirty hack
+    # data = first(elements).fields["displacement"].data
+    # if length(data) > 1
+    #     time0 = data[end-1].first
+    #     dtime = time - time0
+    # end
+
+    for element in elements
+
+        u = element("displacement", time)
+        fill!(Km, 0.0)
+        fill!(f_int, 0.0)
+        fill!(f_ext, 0.0)
+
+        for ip in get_integration_points(element)
+            J, detJ, N, dN = element_info!(bi, element, ip, time)
+            material = ip("material", time)
+            w = ip.weight*detJ
+
+            # Kinematic matrix, linear part
+
+            fill!(BL, 0.0)
+            for i=1:nnodes
+                BL[1, 3*(i-1)+1] = dN[1,i]
+                BL[2, 3*(i-1)+2] = dN[2,i]
+                BL[3, 3*(i-1)+3] = dN[3,i]
+                BL[4, 3*(i-1)+1] = dN[2,i]
+                BL[4, 3*(i-1)+2] = dN[1,i]
+                BL[5, 3*(i-1)+2] = dN[3,i]
+                BL[5, 3*(i-1)+3] = dN[2,i]
+                BL[6, 3*(i-1)+1] = dN[3,i]
+                BL[6, 3*(i-1)+3] = dN[1,i]
+            end
+
+            # Calculate stress response
+            integrate_material!(material)
+            D = material.jacobian
+            S = material.stress + material.dstress
+            #@info("material matrix", D)
+
+            # Material stiffness matrix
+            Km += w*BL'*D*BL
+
+            # Internal force vector
+            f_int += w*BL'*S
+
+            # External force vector
+            for i=1:dim
+                haskey(element, "displacement load $i") || continue
+                b = element("displacement load $i", ip, time)
+                f_ext[i:dim:end] += w*B*vec(N)
+            end
+
+        end
+
+        # add contributions to K, Kg, f
+        gdofs = get_gdofs(problem, element)
+        add!(assembly.K, gdofs, gdofs, Km)
+        add!(assembly.f, gdofs, f_ext - f_int)
+
+    end
+
+    return nothing
+end
+
+function FEMBase.assemble_elements!(problem::Problem{Continuum3D},
+                                    assembly::Assembly,
+                                    elements::Vector{Element{Quad4}},
+                                    time::Float64)
+
+    nnodes = 4
+    ndofs = 3
+    f = zeros(nnodes*ndofs)
+    bi = BasisInfo(Quad4)
+
+    for element in elements
+
+        fill!(f, 0.0)
+
+        for ip in get_integration_points(element)
+
+            J, detJ, N, dN = element_info!(bi, element, ip, time)
+            w = ip.weight*detJ
+
+            if haskey(element, "surface pressure")
+                J = element(ip, time, Val{:Jacobian})'
+                n = cross(J[:,1], J[:,2])
+                n /= norm(n)
+                # sign convention, positive pressure is towards surface
+                p = element("surface pressure", ip, time)
+                f += w*p*vec(n*N)
+            end
+        end
+
+        gdofs = get_gdofs(problem, element)
+        add!(assembly.f, gdofs, f)
+
+    end
+
+    return nothing
+
+end

examples/one_elem_disp_chaboche.jl

@@ -0,0 +1,177 @@
+# This file is a part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/Materials.jl/blob/master/LICENSE
+
+using JuliaFEM, FEMBase, LinearAlgebra
+
+include("continuum.jl")
+include("chaboche.jl")
+
+X = Dict(
+    1 => [0.0, 0.0, 0.0],
+    2 => [1.0, 0.0, 0.0],
+    3 => [1.0, 1.0, 0.0],
+    4 => [0.0, 1.0, 0.0],
+    5 => [0.0, 0.0, 1.0],
+    6 => [1.0, 0.0, 1.0],
+    7 => [1.0, 1.0, 1.0],
+    8 => [0.0, 1.0, 1.0])
+
+body_element = Element(Hex8, (1, 2, 3, 4, 5, 6, 7, 8))
+body_elements = [body_element]
+update!(body_elements, "geometry", X)
+update!(body_elements, "youngs modulus", 200.0e3)
+update!(body_elements, "poissons ratio", 0.3)
+update!(body_elements, "K_n", 100.0)
+update!(body_elements, "n_n", 10.0)
+update!(body_elements, "C_1", 10000.0)
+update!(body_elements, "D_1", 100.0)
+update!(body_elements, "C_2", 50000.0)
+update!(body_elements, "D_2", 1000.0)
+update!(body_elements, "Q", 50.0)
+update!(body_elements, "b", 0.1)
+update!(body_elements, "yield stress", 100.0)
+
+bc_element_1 = Element(Poi1, (1,))
+bc_element_2 = Element(Poi1, (2,))
+bc_element_3 = Element(Poi1, (3,))
+bc_element_4 = Element(Poi1, (4,))
+bc_element_5 = Element(Poi1, (5,))
+bc_element_6 = Element(Poi1, (6,))
+bc_element_7 = Element(Poi1, (7,))
+bc_element_8 = Element(Poi1, (8,))
+
+bc_elements = [bc_element_1, bc_element_2, bc_element_3, bc_element_4,
+               bc_element_5, bc_element_6, bc_element_7, bc_element_8]
+update!(bc_elements, "geometry", X)
+for element in (bc_element_1, bc_element_2, bc_element_3, bc_element_4)
+    update!(element, "displacement 3", 0.0)
+end
+for element in (bc_element_5, bc_element_6, bc_element_7, bc_element_8)
+    update!(element, "displacement 3", 0.0 => 0.0)
+    update!(element, "displacement 3", 1.0 => 5.0e-3)
+    update!(element, "displacement 3", 3.0 => -5.0e-3)
+    update!(element, "displacement 3", 5.0 => 5.0e-3)
+    update!(element, "displacement 3", 7.0 => -5.0e-3)
+    update!(element, "displacement 3", 9.0 => 5.0e-3)
+    update!(element, "displacement 3", 10.0 => 0.0)
+end
+
+update!(bc_element_1, "displacement 1", 0.0)
+update!(bc_element_1, "displacement 2", 0.0)
+update!(bc_element_2, "displacement 2", 0.0)
+update!(bc_element_4, "displacement 1", 0.0)
+
+update!(bc_element_5, "displacement 1", 0.0)
+update!(bc_element_5, "displacement 2", 0.0)
+update!(bc_element_6, "displacement 2", 0.0)
+update!(bc_element_8, "displacement 1", 0.0)
+
+#update!(bc_element_5, "displacement 1", 0.0)
+#update!(bc_element_5, "displacement 2", 0.0)
+#update!(bc_element_5, "displacement 3", 0.0 => 0.0)
+#update!(bc_element_5, "displacement 3", 1.0 => 1.0e-3)
+
+# Initialize material model to integration points
+for ip in get_integration_points(body_element)
+    update!(ip, "yield stress", 0.0 => body_element("yield stress", ip, 0.0))
+    update!(ip, "plastic strain", 0.0 => zeros(6))
+    update!(ip, "stress", 0.0 => zeros(6))
+    update!(ip, "strain", 0.0 => zeros(6))
+    update!(ip, "backstress 1", 0.0 => zeros(6))
+    update!(ip, "backstress 2", 0.0 => zeros(6))
+    update!(ip, "cumulative equivalent plastic strain", 0.0 => 0.0)
+    ip.fields["material"] = field(Chaboche(body_element, ip, 0.0))
+end
+
+body = Problem(Continuum3D, "1 element problem", 3)
+bc = Problem(Dirichlet, "fix displacement", 3, "displacement")
+add_elements!(body, body_elements)
+add_elements!(bc, bc_elements)
+analysis = Analysis(Nonlinear, "solve problem")
+# xdmf = Xdmf("results"; overwrite=true)
+# add_results_writer!(analysis, xdmf)
+add_problems!(analysis, body, bc)
+# time_end = 1.0
+time_end = 10.0
+dtime = 0.05
+
+for problem in get_problems(analysis)
+    initialize!(problem, analysis.properties.time)
+end
+
+while analysis.properties.time < time_end
+    analysis.properties.time += dtime
+    update!(body_element, "displacement", analysis.properties.time => Dict(j => zeros(3) for j in 1:8))
+    @info("time = $(analysis.properties.time)")
+    run!(analysis)
+    for ip in get_integration_points(body_element)
+        material = ip("material", analysis.properties.time)
+        material_matrix = zeros(6,6)
+        stress_vector = zeros(6)
+        calculate_stress!(material, body_element, ip,
+            analysis.properties.time, dtime, material_matrix, stress_vector)
+        update!(ip, "stress", analysis.properties.time => stress_vector)
+        # update!(ip, "stress", analysis.properties.time => copy(material.stress))
+        update!(ip, "plastic strain", analysis.properties.time => copy(material.plastic_strain))
+        update!(ip, "cumulative equivalent plastic strain", analysis.properties.time => copy(material.cumulative_equivalent_plastic_strain))
+        update!(ip, "backstress 1", analysis.properties.time => copy(material.backstress1))
+        update!(ip, "backstress 2", analysis.properties.time => copy(material.backstress2))
+        update!(ip, "yield stress", analysis.properties.time => copy(material.yield_stress))
+    end
+    # update material internal parameters
+end
+
+# close(xdmf)
+
+using Plots
+if true
+    ip1 = first(get_integration_points(body_element))
+    t = range(0, stop=time_end, length=Int(time_end/dtime)+1)
+    s11(t) = ip1("stress", t)[1]
+    s22(t) = ip1("stress", t)[2]
+    s33(t) = ip1("stress", t)[3]
+    s12(t) = ip1("stress", t)[4]
+    s23(t) = ip1("stress", t)[5]
+    s31(t) = ip1("stress", t)[6]
+    e11(t) = ip1("strain", t)[1]
+    e22(t) = ip1("strain", t)[2]
+    e33(t) = ip1("strain", t)[3]
+    s(t) = ip1("stress", t)
+    function vmis(t)
+        s11, s22, s33, s12, s23, s31 = ip1("stress", t)
+        return sqrt(1/2*((s11-s22)^2 + (s22-s33)^2 + (s33-s11)^2 + 6*(s12^2+s23^2+s31^2)))
+    end
+    path = joinpath("one_elem_disp_chaboche", "unitelement_results.rpt")
+    data = readdlm(path, Float64; skipstart=4)
+    t_ = data[:,1]
+    s11_ = data[:,2]
+    s12_ = data[:,3]
+    s13_ = data[:,4]
+    s22_ = data[:,5]
+    s23_ = data[:,6]
+    s33_ = data[:,7]
+    e11_ = data[:,8]
+    e12_ = data[:,9]
+    e13_ = data[:,10]
+    e22_ = data[:,11]
+    e23_ = data[:,12]
+    e33_ = data[:,13]
+    plot(e11.(t), s11.(t), label="\$\\sigma_{11}\$", legend=:topleft,
+        fg_legend=:transparent, bg_legend=:transparent)
+    plot!(e22.(t), s22.(t), label="\$\\sigma_{22}\$")
+    plot!(e33.(t), s33.(t), linecolor=:red, label="\$\\sigma_{33}\$")
+    plot!(e11_, s11_, ls=:dash, label="\$\\sigma_{11} \\quad \\mathrm{Commercial}\$")
+    plot!(e22_, s22_, ls=:dash, label="\$\\sigma_{22} \\quad \\mathrm{Commercial}\$")
+    plot!(e33_, s33_, linecolor=:black, lw=1, ls=:dash,
+        label="\$\\sigma_{33} \\quad \\mathrm{Commercial}\$")
+    title!("Chaboche plasticity model\nOne element model with uniaxial stress")
+    # xlabel!("\$\\varepsilon\$")
+    # ylabel!("\$\\sigma\$")
+    # labels = ["s11" "s22" "s33" "s12" "s23" "s31"]
+    # plot(t, s11, title="stress at integration point 1", label="s11")
+    # plot!(t, s22, label="s22")
+    # plot!(t, s33, label="s33")
+    # plot!(t, s12, label="s12")
+    # plot!(t, s23, label="s23")
+    # plot!(t, s31, label="s31")
+end