In [1]:
using JuliaFEM

using JuliaFEM.Core: Node, Element, Seg2, Tri3, Quad4, Hex8
using JuliaFEM.Core: Problem, FieldProblem, BoundaryProblem, Dirichlet, Elasticity, Mortar
using JuliaFEM.Core: Solver, SparseMatrixCOO
using JuliaFEM.Core: get_elements, update!, calculate_normal_tangential_coordinates!,
get_connectivity, get_field_assembly, get_boundary_problems, handle_overconstraint_error!
using JuliaFEM.Preprocess: parse_aster_med_file, aster_create_elements
import JuliaFEM.Core: solve_linear_system

using PyPlot

ENV["COLUMNS"] = 300


Out[1]:
300

In [2]:
function create_ironing_problem(meshfile="/geometry/3d_ironing/MESH_SPARSE.med")

    mesh = parse_aster_med_file(Pkg.dir("JuliaFEM")*meshfile)
    
    body1 = Problem(Elasticity, "SLAB", 3)
    body1_elements = aster_create_elements(mesh, :SLAB, :HE8)
    update!(body1_elements, "youngs modulus", 288.0)
    update!(body1_elements, "poissons ratio", 1/3)
    push!(body1, body1_elements...)

    body2 = Problem(Elasticity, "DIE", 3)
    body2_elements = aster_create_elements(mesh, :DIE, :HE8)
    update!(body2_elements, "youngs modulus", 288000.0)
    update!(body2_elements, "poissons ratio", 1/3)
    push!(body2, body2_elements...)

    # boundary conditions
    bc1 = Problem(Dirichlet, "bottom surface of slab", 3, "displacement")
    bc1_elements = aster_create_elements(mesh, :GROUND, :QU4)
    update!(bc1_elements, "displacement 1", 0.0)
    update!(bc1_elements, "displacement 2", 0.0)
    update!(bc1_elements, "displacement 3", 0.0)
    push!(bc1, bc1_elements...)

    bc2 = Problem(Dirichlet, "top surface of die", 3, "displacement")
    bc2_elements = aster_create_elements(mesh, :TOP, :QU4)
    #d = [0.0, 0.0, -1.0]
    #update!(bc2_elements, "displacement", Vector{Float64}[d, d, d, d])
    update!(bc2_elements, "displacement 1",  0.0)
    update!(bc2_elements, "displacement 2",  0.0)
    update!(bc2_elements, "displacement 3",  -1.0)
    push!(bc2, bc2_elements...)

    # contact
    bc3 = Problem(Mortar, "contact between slab and die", 3, "displacement")
    bc3_slave_elements = aster_create_elements(mesh, :SLAB_TO_DIE, :QU4)
    bc3_master_elements = aster_create_elements(mesh, :DIE_TO_SLAB, :QU4)
    update!(bc3_slave_elements, "master elements", bc3_master_elements)
    info("normal tangential for first slave element")
    calculate_normal_tangential_coordinates!(bc3_slave_elements, 0.0)
    Q = bc3_slave_elements[1]("normal-tangential coordinates", [0.0, 0.0], 0.0)
    dump(Q)
    push!(bc3, bc3_slave_elements...)
    push!(bc3, bc3_master_elements...)
    info("# of master elements: $(length(bc3_master_elements))")
    info("# of slave elements: $(length(bc3_slave_elements))")

    return body1, body2, bc1, bc2, bc3

end

body1, body2, bc1, bc2, bc3 = create_ironing_problem();


INFO: Found 6 element sets: TOP, DIE, GROUND, SLAB, DIE_TO_SLAB, SLAB_TO_DIE
INFO: normal tangential for first slave element
Array(Float64,(3,3)) 3x3 Array{Float64,2}
INFO: # of master elements: 120
INFO: # of slave elements: 80

In [3]:
body1, body2, bc1, bc2, bc3 = create_ironing_problem()
bc3.properties.inequality_constraints = true
#bc3.properties.normal_condition = :Contact
#bc3.properties.tangential_condition = :Slip
#bc3.properties.minimum_distance = 1
solver = Solver("solve ironing problem.")
push!(solver, body1, body2, bc1, bc2, bc3)
call(solver)


:
 0.0  1.0  0.0
 0.0  0.0  1.0
 1.0  0.0  0.0
Array(Float64,(3,3)) 3x3 Array{Float64,2}:
 0.0  1.0  0.0
 0.0  0.0  1.0
 1.0  0.0  0.0
INFO: Found 6 element sets: TOP, DIE, GROUND, SLAB, DIE_TO_SLAB, SLAB_TO_DIE
INFO: normal tangential for first slave element
INFO: # of master elements: 120
INFO: # of slave elements: 80
INFO: solving linear system of 5 problems.
INFO: PDASS: Starting primal-dual active set strategy to determine active constraints
INFO: PDASS: contact nodes: [304,307,308,309,310,318,322,323,324,325,348,349,350,351,352,353,354,355,356,357,358,359,360,361,362,363,364,365,366,386,387,388,389,390,391,392,393,394,395,396,397,398,399,400,401,402,403,404,550,551,552,553,554,555,556,557,558,559,560,561,562,563,564,565,566,567,568,569,570,571,572,573,574,575,576,577,578,579,580,581,582,583,584,585,586,587,588,589,590,591,592,593,594,595,596,597,598,599,600,601,602,603,604,605,606]
INFO: PDASS: active nodes: Int64[]
INFO: PDASS: inactive nodes: [304,307,308,309,310,318,322,323,324,325,348,349,350,351,352,353,354,355,356,357,358,359,360,361,362,363,364,365,366,386,387,388,389,390,391,392,393,394,395,396,397,398,399,400,401,402,403,404,550,551,552,553,554,555,556,557,558,559,560,561,562,563,564,565,566,567,568,569,570,571,572,573,574,575,576,577,578,579,580,581,582,583,584,585,586,587,588,589,590,591,592,593,594,595,596,597,598,599,600,601,602,603,604,605,606]
INFO: UMFPACK: solved in 0.6064090728759766 seconds. norm = 17.32050807568924
INFO: solving linear system of 5 problems.
INFO: PDASS: Starting primal-dual active set strategy to determine active constraints
INFO: PDASS: contact nodes: [304,307,308,309,310,318,322,323,324,325,348,349,350,351,352,353,354,355,356,357,358,359,360,361,362,363,364,365,366,386,387,388,389,390,391,392,393,394,395,396,397,398,399,400,401,402,403,404,550,551,552,553,554,555,556,557,558,559,560,561,562,563,564,565,566,567,568,569,570,571,572,573,574,575,576,577,578,579,580,581,582,583,584,585,586,587,588,589,590,591,592,593,594,595,596,597,598,599,600,601,602,603,604,605,606]
INFO: PDASS: active nodes: [348,349,350,351,352,353,354,355,356,357,386,387,388,389,390,391,392,393,394,395,550,551,552,553,554,555,556,557,558,559,569,570,571,572,573,574,575,576,577,578,588,589,590,591,592,593,594,595,596,597]
INFO: PDASS: inactive nodes: [304,307,308,309,310,318,322,323,324,325,358,359,360,361,362,363,364,365,366,396,397,398,399,400,401,402,403,404,560,561,562,563,564,565,566,567,568,579,580,581,582,583,584,585,586,587,598,599,600,601,602,603,604,605,606]
INFO: UMFPACK: solved in 0.24304413795471191 seconds. norm = 18.1021717204183
INFO: solving linear system of 5 problems.
INFO: PDASS: Starting primal-dual active set strategy to determine active constraints
INFO: PDASS: contact nodes: [304,307,308,309,310,318,322,323,324,325,348,349,350,351,352,353,354,355,356,357,358,359,360,361,362,363,364,365,366,386,387,388,389,390,391,392,393,394,395,396,397,398,399,400,401,402,403,404,550,551,552,553,554,555,556,557,558,559,560,561,562,563,564,565,566,567,568,569,570,571,572,573,574,575,576,577,578,579,580,581,582,583,584,585,586,587,588,589,590,591,592,593,594,595,596,597,598,599,600,601,602,603,604,605,606]
INFO: PDASS: active nodes: [349,350,351,352,353,354,355,356,387,388,389,390,391,392,393,394,551,552,553,554,555,556,557,558,570,571,572,573,574,575,576,577,589,590,591,592,593,594,595,596]
INFO: PDASS: inactive nodes: [304,307,308,309,310,318,322,323,324,325,348,357,358,359,360,361,362,363,364,365,366,386,395,396,397,398,399,400,401,402,403,404,550,559,560,561,562,563,564,565,566,567,568,569,578,579,580,581,582,583,584,585,586,587,588,597,598,599,600,601,602,603,604,605,606]
INFO: UMFPACK: solved in 0.12660002708435059 seconds. norm = 18.174500081971626
INFO: solving linear system of 5 problems.
INFO: PDASS: Starting primal-dual active set strategy to determine active constraints
INFO: PDASS: contact nodes: [304,307,308,309,310,318,322,323,324,325,348,349,350,351,352,353,354,355,356,357,358,359,360,361,362,363,364,365,366,386,387,388,389,390,391,392,393,394,395,396,397,398,399,400,401,402,403,404,550,551,552,553,554,555,556,557,558,559,560,561,562,563,564,565,566,567,568,569,570,571,572,573,574,575,576,577,578,579,580,581,582,583,584,585,586,587,588,589,590,591,592,593,594,595,596,597,598,599,600,601,602,603,604,605,606]
INFO: PDASS: active nodes: [349,350,351,352,353,354,355,356,387,388,389,390,391,392,393,394,551,552,553,554,555,556,557,558,570,571,572,573,574,575,576,577,589,590,591,592,593,594,595,596]
INFO: PDASS: inactive nodes: [304,307,308,309,310,318,322,323,324,325,348,357,358,359,360,361,362,363,364,365,366,386,395,396,397,398,399,400,401,402,403,404,550,559,560,561,562,563,564,565,566,567,568,569,578,579,580,581,582,583,584,585,586,587,588,597,598,599,600,601,602,603,604,605,606]
INFO: UMFPACK: solved in 0.14113306999206543 seconds. norm = 18.174500081971626
INFO: Converged in 4 iterations.
Out[3]:
true

In [4]:
using JuliaFEM.Postprocess: xdmf_new_model, xdmf_new_temporal_collection, xdmf_new_grid,
                            xdmf_new_mesh!, xdmf_new_nodal_field!, xdmf_save_model
using JuliaFEM.Core: Element

function xdmf_dump(all_elements, filename="/tmp/xdmf_result.xmf")
    info("$(length(all_elements)) elements.")
    xdoc, xmodel = xdmf_new_model()
    coll = xdmf_new_temporal_collection(xmodel)
    grid = xdmf_new_grid(coll; time=0.0)

    Xg = Dict{Int64, Vector{Float64}}()
    ug = Dict{Int64, Vector{Float64}}()
    nids = Dict{Int64, Int64}()
    for element in all_elements
        conn = get_connectivity(element)
        for (i, c) in enumerate(conn)
            nids[c] = c
        end
        X = element("geometry", 0.0)
        for (i, c) in enumerate(conn)
            Xg[c] = X[i]
        end
        haskey(element, "displacement") || continue
        u = element("displacement", 0.0)
        for (i, c) in enumerate(conn)
            ug[c] = u[i]
        end
    end
    perm = sort(collect(keys(Xg)))
    nodes = Vector{Float64}[Xg[i] for i in perm]
    disp = Vector{Float64}[ug[i] for i in perm]
    nids = Int[nids[i] for i in perm]
    inids = Dict{Int64, Int64}()
    for (i, nid) in enumerate(nids)
        inids[nid] = i
    end
    elements = []
    for element in all_elements
        isa(element, Element{Hex8}) || continue
        conn = get_connectivity(element)
        nconn = [inids[i] for i in conn]
        push!(elements, (:Hex8, nconn))
    end

    xdmf_new_mesh!(grid, nodes, elements)
    xdmf_new_nodal_field!(grid, "displacement", disp)
    xdmf_save_model(xdoc, filename)
    info("model dumped to $filename")
end

xdmf_dump([body1.elements; body2.elements], "/tmp/ironing.xmf")


INFO: 360 elements.
INFO: XDFM: ndim = 2160
INFO: model dumped to /tmp/ironing.xmf

In [ ]: