Merge pull request #27 from axla-io/optimize-performance-element-cs

Optimize performance element cs

Project.toml

@@ -4,30 +4,13 @@ author = ["Axel Larsson"]
 version = "0.2.0"
 
 [deps]
+StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6"
+StaticGraphs = "4c8beaf5-199b-59a0-a7f2-21d17de635b6"
 DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
 DiffEqCallbacks = "459566f4-90b8-5000-8ac3-15dfb0a30def"
-DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
-FiniteDiff = "6a86dc24-6348-571c-b903-95158fe2bd41"
-Flux = "587475ba-b771-5e3f-ad9e-33799f191a9c"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
-GraphRecipes = "bd48cda9-67a9-57be-86fa-5b3c104eda73"
-Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6"
-LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
-ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
-NBodySimulator = "0e6f8da7-a7fc-5c8b-a220-74e902c310f9"
-Optim = "429524aa-4258-5aef-a3af-852621145aeb"
-Optimization = "7f7a1694-90dd-40f0-9382-eb1efda571ba"
-OptimizationNLopt = "4e6fcdb7-1186-4e1f-a706-475e75c168bb"
-OptimizationOptimJL = "36348300-93cb-4f02-beb5-3c3902f8871e"
-OptimizationOptimisers = "42dfb2eb-d2b4-4451-abcd-913932933ac1"
-OptimizationPolyalgorithms = "500b13db-7e66-49ce-bda4-eed966be6282"
-OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
-Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
-SciMLSensitivity = "1ed8b502-d754-442c-8d5d-10ac956f44a1"
-StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
-StaticGraphs = "4c8beaf5-199b-59a0-a7f2-21d17de635b6"
-SteadyStateDiffEq = "9672c7b4-1e72-59bd-8a11-6ac3964bc41f"
-Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
 
 [extras]
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

examples/beam_steady_state.jl

@@ -5,9 +5,8 @@ using DynamicRelaxation
 using Graphs
 using StaticGraphs
 using DiffEqCallbacks
-using DifferentialEquations
+using NBodySimulator
 using SteadyStateDiffEq
-
 # Define a simple graph system
 n_elem = 17
 n_pt = n_elem + 1
@@ -18,31 +17,36 @@ system = default_system(graph, Node6DOF, :catenary, n_pt)
 ext_f = uniform_load(Pz(-10_000, system), system)
 
 # Set parameters
-maxiters = 300
+maxiters = 500
 dt = 0.01
 tspan = (0.0, 10.0)
 
 # Create problem
-simulation = RodSimulation{StructuralGraphSystem{Node6DOF},Float64,eltype(ext_f)}(system, tspan, dt, ext_f)
-prob = ODEProblem(simulation)
+simulation = RodSimulation(system, tspan, dt)
+prob = ODEProblem(simulation, ext_f)
+
 ssprob = SteadyStateProblem(prob)
 
-# Create decay callback
+# Create callback for damping
 c = 0.7
 (u0, v0, n, u_len, v_len) = get_u0(simulation)
-(dx_ids, dr_ids, v_ids, ω_ids) = get_vel_ids(u_len, v_len)
+(dx_ids, dr_ids, v_ids, ω_ids) = get_vel_ids(u_len, v_len, system)
 v_decay!(integrator) = velocitydecay!(integrator, vcat(v_ids, ω_ids), c)
-cb = PeriodicCallback(v_decay!, 3 * dt; initial_affect=true)
+cb1 = PeriodicCallback(v_decay!, 2 * dt; initial_affect = true)
+
+# Create callback for termination
+abstol = reltol = 1e-2
+ktc(integrator, abstol, reltol, min_t) = ke_termination_cond(integrator, abstol, reltol, min_t, v_ids)
+cb2 = TerminateSteadyState(abstol, reltol, ktc, min_t = 20*dt)
 
-# Set termination condition
-cond = NLSolveTerminationCondition(NLSolveTerminationMode.AbsSafe; abstol = 1e-1)
+cb = CallbackSet(cb1,cb2)
 
 # Set algorithm for solver
 alg = RK4()
 
 # Solve problem
-@time sol = solve(ssprob, DynamicSS(alg, termination_condition = cond),  maxiters=maxiters, callback = cb);
+@time sol = solve(ssprob, DynamicSS(alg), custom_termination_cond = true, maxiters=maxiters, callback = cb);
 
 # Plot final state
-u_final = get_state(sol.u, u_len)
+u_final = get_state(sol.u, u_len, simulation)
 plot(u_final[1, :], u_final[3, :])

examples/cantileverbeam.jl

@@ -0,0 +1,63 @@
+using DynamicRelaxation
+using OrdinaryDiffEq
+using Graphs
+using StaticGraphs
+using DiffEqCallbacks
+using Plots
+using LinearAlgebra
+
+n_elem = 6
+n_pt = n_elem + 1
+graph = StaticGraph(path_graph(n_pt))
+system = default_system(graph, Node6DOF, :cantilever, n_pt)
+
+# Set loads
+ext_f = point_loads([Pz(-1_000, system)], [n_pt], system)
+# Set parameters
+maxiters = 5000
+dt = 10
+tspan = (0.0, Inf)
+
+# Create problem
+simulation = RodSimulation(system, tspan, dt)
+prob = ODEProblem(simulation, ext_f)
+
+# Create callback
+c = 0.991
+(u0, v0, n, u_len, v_len) = get_u0(simulation)
+(dx_ids, dr_ids, v_ids, ω_ids) = get_vel_ids(u_len, v_len, system)
+v_decay!(integrator) = velocitydecay!(integrator, vcat(v_ids, ω_ids), c)
+cb1 = PeriodicCallback(v_decay!, 2 * dt; initial_affect = true)
+
+res = zero(prob.u0)
+prob.f(res, prob.u0, prob.p, 1.0)
+norm(res)
+# Set algorithm for solver
+alg = RK4()
+
+# Solve problem
+@time sol = solve(prob, alg, dt = simulation.dt, maxiters = maxiters, callback = cb1, 
+                  verbose = false);
+sol.retcode
+# Plot final state
+u_final = get_state(sol.u[end], u_len, simulation)
+
+# Select frames for animation
+itt = generate_range(100, 1, length(sol.u))
+u_red = sol.u[itt]
+
+gr()
+# Loop over the time values and create a plot for each frame
+anim = @animate for i in axes(u_red, 1)
+    u_final = get_state(u_red[i], u_len, simulation)
+    plot(u_final[1, :], u_final[3, :], xlims = (0,10), ylims = (-5, 1), label = "",c = :black, lw = 3.0)
+    quiver!([u_final[1, end]], [u_final[3, end]], quiver=([0], [-1]), color=:blue, label="Arrow Tip", lw = 2.5)
+    scatter!(u_final[1, 2:end], u_final[3, 2:end], label = "",markercolor=:white, markerstrokecolor=:black, markersize=5, lw = 3.0)
+end
+
+res = zero(prob.u0)
+prob.f(res, sol.u[end], prob.p, 0.0)
+norm(res)
+
+gif(anim, "cantilever.gif")
+