julia_help.jl

using Gridap
using Gridap.Geometry
const E = 70.0e9
const ν = 0.33
const λ = (E * ν) / ((1 + ν) * (1 - 2 * ν))
const μ = E / (2 * (1 + ν))
const density = 2710 # kg/m^3
σ(ε) = λ * tr(ε) * one(ε) + 2 * μ * ε

# for making the FEA grid. 
L = 7
W = 2
H = 3
mult = 3

# ODE Solver config
Δt = 0.01
θ = 0.5
time_start = 0.0
time_end = 3.0

function make_model()
    model = CartesianDiscreteModel((0, W, 0, L, 0, H), (W * mult, L * mult, H * mult))
    labels = get_face_labeling(model)
    add_tag_from_tags!(labels, "fixed", ["tag_24",
        "tag_03",
        "tag_04",
        "tag_07",
        "tag_08",
        "tag_10",
        "tag_12",
        "tag_19",
        "tag_20"
    ])
    add_tag_from_tags!(labels, "load", ["tag_23"])
    writevtk(model, "model_with_fixed")
    model
end

function make_U0_solution(model; order=1)
    reffe = ReferenceFE(lagrangian, VectorValue{3,Float64}, order)
    V0 = TestFESpace(model, reffe;
        conformity=:H1,
        dirichlet_tags=["fixed", "load"],
        dirichlet_masks=[(true, true, true), (true, false, false)])

    g2_disp(x) = VectorValue(0.05, 0.0, 0.0)
    g1_fixed(x) = VectorValue(0.0, 0.0, 0.0)

    # From functions `g1` and `g2`, we define the trial space as follows:
    U = TrialFESpace(V0, [g1_fixed, g2_disp])

    degree = 2 * order
    Ω = Triangulation(model)
    dΩ = Measure(Ω, degree)
    a(u, v) = ∫(ε(v) ⊙ (σ ∘ ε(u))) * dΩ
    l(v) = 0
    op = AffineFEOperator(a, l, U, V0)

    uh = solve(op)
    writevtk(Ω, "u_0_for_transcient", cellfields=["uh" => uh, "epsi" => ε(uh), "sigma" => σ ∘ ε(uh)])
    uh
end

function solve_transcient(model, U0; order=1)
    reffe = ReferenceFE(lagrangian, VectorValue{3,Float64}, order)
    V = TestFESpace(model, reffe;
        conformity=:H1,
        dirichlet_tags=["fixed"],
        dirichlet_masks=[(true, true, true)]
    )
    g1_fixed(x, t) = VectorValue(0.0, 0.0, 0.0)
    g1_fixed(t) = g1_fixed(0, 0)

    U = TransientTrialFESpace(V, g1_fixed)

    degree = 2 * order
    Ω = Triangulation(model)
    dΩ = Measure(Ω, degree)

    # try 5
    m(t, u, v) = ∫(density * outer(v, u))dΩ
    c(t, u, v) = ∫(0.0 * u ⋅ v)dΩ
    a(t, u, v) = ∫(ε(v) ⊙ (σ ∘ ε(u)))dΩ
    b(t, u, v) = ∫(0.0 * u ⋅ v)dΩ
    b(t, v) = b(t, 0, v)
    op = TransientAffineFEOperator(m, c, a, b, U, V)
    

    linear_solver = LUSolver()
    ode_solver = Newmark(linear_solver,Δt,0.5,0.25) # ode_solver = ThetaMethod(linear_solver, Δt, θ)

    V0 =  interpolate_everywhere( VectorValue(0.0,0.0,0.0),U(0.0))
    A0 =  interpolate_everywhere( VectorValue(0.0,0.0,0.0),U(0.0))
    u_transcient = solve(ode_solver, op, (U0, V0, A0), time_start, time_end)

    # cut out saving to .vtu file to help with debugging. 
    function run_solver_t(uₕₜ)
        for (uₕ, t) in uₕₜ
            println(t)
        end
    end

    run_solver_t(u_transcient)
end

model = make_model()
U0 = make_U0_solution(model)
U_trans = solve_transcient(model, U0)