Merge pull request #59 from JuliaGNI/rigid_body

Rigid body

.github/workflows/Documenter.yml

@@ -15,7 +15,7 @@ jobs:
     steps:
       - uses: actions/checkout@v3
       - uses: julia-actions/setup-julia@latest
-      - name: Configure doc environment
+      - name: Install imagemagick, poppler-utils and texlive
         run: |
           sudo apt-get install imagemagick
           sudo apt-get install poppler-utils

docs/Makefile

@@ -20,7 +20,7 @@ documenter:
 images:
 	$(MAKE) all   -C src/images
 	$(MAKE) clean -C src/images
-
+	
 clean:
 	$(MAKE) empty -C src/images
 	rm -Rf build

docs/Project.toml

@@ -1,4 +1,5 @@
 [deps]
+CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 DocumenterCitations = "daee34ce-89f3-4625-b898-19384cb65244"
 GeometricEquations = "c85262ba-a08a-430a-b926-d29770767bf2"

docs/make.jl

@@ -34,6 +34,7 @@ makedocs(;
              "Point Vortices"              => "point_vortices.md",
              "Inner Solar System"          => "inner_solar_system.md",
              "Outer Solar System"          => "outer_solar_system.md",
+             "Rigid body"                  => "rigid_body.md",
              "Toda Lattice"                => "toda_lattice.md",
             ]
 )

docs/src/GeometricProblems.bib

@@ -25,4 +25,14 @@ @book{hairer2006geometric
   year={2006},
   publisher={Springer},
   address={Berlin}
-}
+}
+
+@article{bajars2023locally,
+  title={Locally-symplectic neural networks for learning volume-preserving dynamics},
+  author={Baj{\=a}rs, J{\=a}nis},
+  journal={Journal of Computational Physics},
+  volume={476},
+  pages={111911},
+  year={2023},
+  publisher={Elsevier}
+}

docs/src/images/Makefile

@@ -10,6 +10,18 @@ png:
 	pdftocairo  -png -r 250 -transp -singlefile  coupled_harmonic_oscillator.pdf        coupled_harmonic_oscillator              
 	pdftocairo  -png -r 250 -transp -singlefile  coupled_harmonic_oscillator_dark.pdf   coupled_harmonic_oscillator_dark       
 
+MYDIR := .
+convert_with_sips: $(MYDIR)/*.pdf
+	for file in $^ ; do \
+		sips --setProperty format png --resampleHeightWidthMax 600 $${file} --out $${file%.*}.png ; \
+	done
+
+# this is converting pdfs to pngs using pdftocairo (linux version)
+convert_with_pdftocairo: $(MYDIR)/*.pdf
+	for file in $^ ; do \
+		pdftocairo -png -r 250 -transp -singlefile $${file} $${file%.*} ; \
+	done
+
 clean:
 	rm -f *.aux
 	rm -f *.fdb_latexmk

docs/src/index.md

@@ -38,6 +38,7 @@ Pages = [
 "point_vortices.md",
 "inner_solar_system.md",
 "outer_solar_system.md",
+"rigid_body.md",
 "toda_lattice.md"
 ]
 Depth = 1

docs/src/rigid_body.md

@@ -0,0 +1,61 @@
+# The rigid body 
+
+```@example 
+using GeometricProblems.RigidBody: odeensemble
+using GeometricIntegrators: integrate, ImplicitMidpoint
+using GeometricEquations: EnsembleProblem
+using GeometricSolutions: GeometricSolution
+using CairoMakie
+
+ics = [
+        [sin(1.1), 0., cos(1.1)],
+        [sin(2.1), 0., cos(2.1)],
+        [sin(2.2), 0., cos(2.2)],
+        [0., sin(1.1), cos(1.1)],
+        [0., sin(1.5), cos(1.5)], 
+        [0., sin(1.6), cos(1.6)]
+        ]
+
+ensemble_problem = odeensemble(ics)
+ensemble_solution = integrate(ensemble_problem, ImplicitMidpoint())
+
+function plot_geometric_solution!(p, solution::GeometricSolution; kwargs...)
+    lines!(p, solution.q[:, 1].parent, solution.q[:, 2].parent, solution.q[:, 3].parent; kwargs...)
+end
+
+function sphere(r, C)   # r: radius; C: center [cx,cy,cz]
+           n = 100
+           u = range(-π, π; length = n)
+           v = range(0, π; length = n)
+           x = C[1] .+ r * cos.(u) * sin.(v)'
+           y = C[2] .+ r * sin.(u) * sin.(v)'
+           z = C[3] .+ r * ones(n) * cos.(v)'
+           return x, y, z
+       end
+
+fig, ax, plt = surface(sphere(1., [0., 0., 0.])..., alpha = .6)
+
+for (i, solution) in zip(1:length(ensemble_solution), ensemble_solution.s)
+    plot_geometric_solution!(ax, solution; label = "trajectory "*string(i), linewidth=2)
+end
+
+fig
+```
+
+
+## Library functions
+
+```@docs 
+GeometricProblems.RigidBody
+```
+
+```@autodocs
+Modules = [GeometricProblems.RigidBody]
+Order = [:constant, :type, :macro, :function]
+```
+
+```@bibliography
+Pages = []
+
+bajars2023locally
+```

src/GeometricProblems.jl

@@ -31,6 +31,7 @@ module GeometricProblems
     include("pendulum.jl")
     include("point_vortices.jl")
     include("point_vortices_linear.jl")
+    include("rigid_body.jl")
     include("toda_lattice.jl")
 
 end

src/rigid_body.jl

@@ -0,0 +1,57 @@
+@doc raw"""
+# Rigid body
+
+```math
+\begin{aligned}
+    \dot{x} = Ayz \\
+    \dot{y} = Bxz \\ 
+    \dot{z} = Cxy
+\end{aligned},
+```
+
+where ``A = (I_2 - I_3)/(I_2I_3)``, ``B = (I_3 - I_1)/(I_3I_1)`` and ``C = (I_1 - I_2)/(I_1I_2)``; ``I_{\cdot}`` are the *principal components of inertia*. 
+
+The initial condition and the default parameters are taken from [bajars2023locally](@cite).
+"""
+module RigidBody
+
+    using GeometricEquations 
+    using GeometricSolutions
+    using Parameters 
+
+    export odeproblem, odeensemble
+
+    const tspan = (0.0, 100.0)
+    const tstep = 0.1
+
+    const default_parameters = (
+        I₁ = 2.,
+        I₂ = 1.,
+        I₃ = 2. / 3.
+    )
+
+    const q₀ = [cos(1.1), 0., sin(1.1)]
+    const q₁ = [cos(2.1), 0., sin(2.1)]
+    const q₂ = [cos(2.2), 0., sin(2.2)]
+
+    function rigid_body_v(v, t, q, params)
+        @unpack I₁, I₂, I₃ = params
+        A = (I₂ - I₃) / (I₂ * I₃)
+        B = (I₃ - I₁) / (I₃ * I₁)
+        C = (I₁ - I₂) / (I₁ * I₂)
+        v[1] = A * q[2] * q[3]
+        v[2] = B * q[1] * q[3]
+        v[3] = C * q[1] * q[2]
+
+        nothing
+    end
+
+    function odeproblem(q₀ = q₀; tspan = tspan, tstep = tstep, parameters = default_parameters)
+        ODEProblem(rigid_body_v, tspan, tstep, q₀; parameters = parameters)
+    end
+
+    function odeensemble(samples = [q₀, q₁, q₂]; parameters = default_parameters, tspan = tspan, tstep = tstep)
+        ODEEnsemble(rigid_body_v, tspan, tstep, samples; parameters = parameters)
+    end
+
+end

test/rigid_body_test.jl

@@ -0,0 +1,5 @@
+using GeometricProblems.RigidBody
+using Test 
+
+@test_nowarn odeproblem()
+@test_nowarn odeensemble()

test/runtests.jl

@@ -14,3 +14,4 @@ using SafeTestsets
 @safetestset "Lotka-Volterra 4D (Lagrangian)                                                  " begin include("lotka_volterra_4d_lagrangian_tests.jl") end
 @safetestset "Point Vortices                                                                  " begin include("point_vortices_tests.jl") end
 @safetestset "Point Vortices (linear)                                                         " begin include("point_vortices_linear_tests.jl") end
+@safetestset "Rigid Body                                                                      " begin include("rigid_body_test.jl") end