Merge pull request #266 from marinlauber/RigidBodies

Add `RigidMap` for fluid-(rigid)body interaction

Author: weymouth

Project.toml

@@ -4,6 +4,7 @@ authors = ["Gabriel Weymouth <gabriel.weymouth@gmail.com>"]
 version = "1.5.2"
 
 [deps]
+ConstructionBase = "187b0558-2788-49d3-abe0-74a17ed4e7c9"
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
 EllipsisNotation = "da5c29d0-fa7d-589e-88eb-ea29b0a81949"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
@@ -38,6 +39,7 @@ WaterLilyReadVTKExt = "ReadVTK"
 WaterLilyWriteVTKExt = "WriteVTK"
 
 [compat]
+ConstructionBase = "1.6.0"
 DocStringExtensions = "0.9"
 EllipsisNotation = "1.8"
 ForwardDiff = "0.10.18, 1"

src/AutoBody.jl

@@ -1,28 +1,22 @@
 """
-    AutoBody(sdf,map=(x,t)->x; compose=true) <: AbstractBody
+    AutoBody(sdf,map=(x,t)->x) <: AbstractBody
 
   - `sdf(x::AbstractVector,t::Real)::Real`: signed distance function
   - `map(x::AbstractVector,t::Real)::AbstractVector`: coordinate mapping function
-  - `compose::Bool=true`: Flag for composing `sdf=sdf∘map`
 
 Implicitly define a geometry by its `sdf` and optional coordinate `map`. Note: the `map`
-is composed automatically if `compose=true`, i.e. `sdf(x,t) = sdf(map(x,t),t)`.
-Both parameters remain independent otherwise. It can be particularly heplful to set
-`compose=false` when adding mulitple bodies together to create a more complex one.
+is composed automatically i.e. `sdf(body::AutoBody,x,t) = body.sdf(body.map(x,t),t)`.
 """
 struct AutoBody{F1<:Function,F2<:Function} <: AbstractBody
     sdf::F1
     map::F2
-    function AutoBody(sdf, map=(x,t)->x; compose=true)
-        comp(x,t) = compose ? sdf(map(x,t),t) : sdf(x,t)
-        new{typeof(comp),typeof(map)}(comp, map)
-    end
 end
+AutoBody(sdf, map=(x,t)->x) = AutoBody(sdf, map)
 
 """
-    d = sdf(body::AutoBody,x,t) = body.sdf(x,t)
+    d = sdf(body::AutoBody,x,t) = body.sdf(body.map(x,t),t)
 """
-sdf(body::AutoBody,x,t=0;kwargs...) = body.sdf(x,t)
+@inline sdf(body::AutoBody,x,t=0;kwargs...) = body.sdf(body.map(x,t),t)
 
 using ForwardDiff
 """
@@ -35,9 +29,9 @@ Skips the `n,V` calculation when `d²>fastd²`.
 """
 function measure(body::AutoBody,x,t;fastd²=Inf)
     # eval d=f(x,t), and n̂ = ∇f
-    d = body.sdf(x,t)
+    d = sdf(body,x,t)
     d^2>fastd² && return (d,zero(x),zero(x)) # skip n,V
-    n = ForwardDiff.gradient(x->body.sdf(x,t), x)
+    n = ForwardDiff.gradient(x->sdf(body,x,t), x)
     any(isnan.(n)) && return (d,zero(x),zero(x))
 
     # correct general implicit fnc f(x₀)=0 to be a pseudo-sdf

src/RigidMap.jl

@@ -0,0 +1,51 @@
+"""
+    RigidMap(center, θ) <: AbstractBody
+
+  - `x₀::SVector{D}`: coordinate of the center of the body
+  - `θ::Union{Real, SVector{3}}`: rotation (single angle in 2D, and in 3D these are the rotation angle around
+                                  the x, y, and z axes respectively.)
+  - `V::SVector{D}=zero(center)`: linear velocity of the center
+  - `xₚ::SVector{D}=zero(center)`: offset of the pivot point compared to center
+  - `ω::Union{Real, SVector{3}}=zero(θ)`: angular velocity (scalar in 2D, vector in 3D)
+
+Define a `RigidMap` for any `AbstractBody` using rigid body motion parameters.
+
+RigidMap updates are computed externally via a set of ODEs and then updated in the
+simulation loop following:
+```julia
+using WaterLily,StaticArrays
+body = AutoBody((x,t)->sqrt(sum(abs2,x))-4,RigidMap(SA{Float32}[16,16],0.f0;ω=0.1f0))
+sim = Simulation((32,32),(1,0),8;body)
+for n in 1:10
+    # update body motion (example: constant angular velocity)
+    θ = sim.body.map.θ + sim.body.map.ω*sim.flow.Δt[end]
+    sim.body = setmap(sim.body; θ)
+    # remeasure and step
+    sim_step!(sim;remeasure=true)
+end
+```
+"""
+struct RigidMap{A<:AbstractVector,R} <: Function
+    x₀ :: A   # center of translation
+    θ  :: R   # rotation (angle in 2D, euler angles in 3D)
+    xₚ :: A   # rotation offset
+    V  :: A   # linear velocity of the center
+    ω  :: R   # angular velocity (scalar in 2D, vector in 3D)
+end
+RigidMap(x₀::SVector,θ;xₚ=zero(x₀),V=zero(x₀),ω=zero(θ)) = RigidMap(x₀, θ, xₚ, V, ω)
+
+# this is the function map(x,t) AND derivative(t->map(x,t),t)
+(map::RigidMap)(x::SVector,t=0) = rotation(map.θ)*(x-map.x₀-map.xₚ)+map.xₚ
+(map::RigidMap)(x::SVector,::ForwardDiff.Dual{Tag}) where Tag = Dual{Tag}.(map(x),-rotation(map.θ)*(map.V+map.ω×(x-map.x₀-map.xₚ)))
+
+# cross product in 2D and rotation matrix in 2D and 3D
+import WaterLily: ×
+×(a::Number,b::SVector{2,T}) where T = a*SA[-b[2],b[1]]
+rotation(θ::T) where T = SA{T}[cos(θ) sin(θ); -sin(θ) cos(θ)]
+rotation(θ::SVector{3,T}) where T = SA{T}[cos(θ[3])*cos(θ[2]) cos(θ[3])*sin(θ[2])*sin(θ[1])+sin(θ[3])*cos(θ[1]) -cos(θ[3])*sin(θ[2])*cos(θ[1])+sin(θ[3])*sin(θ[1]);
+                                         -sin(θ[3])*cos(θ[2]) -sin(θ[3])*sin(θ[2])*sin(θ[1])+cos(θ[3])*cos(θ[1]) sin(θ[3])*sin(θ[2])*cos(θ[1])+cos(θ[3])*sin(θ[1]);
+                                                sin(θ[2])                         -cos(θ[2])*sin(θ[1])                               cos(θ[2])*cos(θ[1])]
+import ConstructionBase: setproperties
+setmap(body::AbstractBody; kwargs...) = setproperties(body,map=setproperties(body.map; kwargs...))
+setmap(body::SetBody; kwargs...) = SetBody(body.op,setmap(body.a; kwargs...),setmap(body.b; kwargs...))
+setmap(body::NoBody; kwargs...) = NoBody()

src/WaterLily.jl

@@ -30,6 +30,10 @@ include("Metrics.jl")
 export MeanFlow,update!,uu!,uu
 
 abstract type AbstractSimulation end
+
+include("RigidMap.jl")
+export RigidMap,setmap
+
 """
     Simulation(dims::NTuple, uBC::Union{NTuple,Function}, L::Number;
                U=norm2(Uλ), Δt=0.25, ν=0., ϵ=1, g=nothing,

test/maintests.jl

@@ -257,7 +257,7 @@ end
     body = AutoBody(circ)
     for i in 2:20
         body += AutoBody(circ,(x,t)->x-rand(2))
-        @test sizeof(body) == i
+        @test sizeof(body) ≤ i
     end
 
     # test curvature, 2D and 3D
@@ -596,4 +596,58 @@ end
         @test all(meanflow1.t .== meanflow2.t)
     end
     @test_nowarn rm(test_dir, recursive=true)
+end
+@testset "RigidMap.jl" begin
+    for T ∈ (Float32,Float64)
+        # initialize a rigid body
+        sdf(x,t) = sqrt(sum(abs2,x))-1
+        body = AutoBody(sdf, RigidMap(SA{T}[0,0],T(0)))
+        # check sdf
+        @test all(measure(body,SA{T}[1.5,0],0) .≈ (1/2,SA{T}[1,0],SA{T}[0,0]))
+        # rotate and add linear velocity
+        body = setmap(body;θ=T(π/4),V=SA{T}[1.0,0])
+        # check sdf and velocity
+        @test all(measure(body,SA{T}[1.5,0],0) .≈ (1/2,SA{T}[1,0],SA{T}[1,0]))
+        # add angular velocity
+        body = setmap(body;ω=T(0.1))
+        @test all(measure(body,SA{T}[1.5,0],0) .≈ (1/2,SA{T}[1,0],SA{T}[1,1.5*0.1]))
+        # 3D rigid body
+        body3D = AutoBody(sdf, RigidMap(SA{T}[0,0,0],SA{T}[0,0,0];xₚ=SA{T}[-.5,0,0]))
+        @test all(measure(body3D,SA{T}[1.5,0,0],0) .≈ (1/2,SA{T}[1,0,0],SA{T}[0,0,0]))
+        # test rotations about x, y, and z
+        # rotate by 180 degrees about x-axis, should not change
+        body3D = setmap(body3D;θ=SA{T}[π,0,0])
+        @test all(measure(body3D,SA{T}[1.5,0,0],0) .≈ (1/2,SA{T}[1,0,0],SA{T}[0,0,0]))
+        # now rotate by 180 around y=axis, should invert z-component of normal
+        body3D = setmap(body3D;θ=SA{T}[0,π,0],V=SA{T}[1.0,0,0])
+        @test all(measure(body3D,SA{T}[1.5,0,0],0) .≈ (1.5,SA{T}[1,0,0],SA{T}[1,0,0]))
+        body3D = setmap(body3D;θ=SA{T}[0,0,π],V=SA{T}[1.0,0,0])
+        @test all(measure(body3D,SA{T}[1.5,0,0],0) .≈ (1.5,SA{T}[1,0,0],SA{T}[1,0,0]))
+        # 3D rigid body with linear and angular velocity
+        body3D = setmap(body3D;θ=SA{T}[0,0,0],V=SA{T}[1.0,0,0],ω=SA{T}[0,0,0.1])
+        @test all(measure(body3D,SA{T}[1.5,0,0],0) .≈ (1/2,SA{T}[1,0,0],SA{T}[1,0.2,0]))
+        @test all(measure(body3D,SA{T}[0,1.5,0],0) .≈ (1/2,SA{T}[0,1,0],SA{T}[0.85,0.05,0]))
+        @test all(measure(body3D,SA{T}[1.5,1.5,1.5],0) .≈ (√(3*(1.5^2))-1,SA{T}[√(1/3),√(1/3),√(1/3)],SA{T}[.85,0.2,0]))
+        # three 3D rotations
+        body3D = setmap(body3D;V=SA{T}[1.0,0,0],ω=SA{T}[0,-0.1,0.1])
+        @test all(measure(body3D,SA{T}[1.5,0,0],0) .≈ (1/2,SA{T}[1,0,0],SA{T}[1,0.2,0.2]))
+        @test all(measure(body3D,SA{T}[0,1.5,1.5],0) .≈ (√(2*(1.5^2))-1,SA{T}[0,√(1/2),√(1/2)],SA{T}[0.7,0.05,0.05]))
+        # test for a SetMap
+        body = AutoBody(sdf, RigidMap(SA{T}[0,0],T(0))) +AutoBody(sdf, RigidMap(SA{T}[1,1],T(0)))
+        body = setmap(body;θ=T(π/4),V=SA{T}[1.0,0])
+        @test all(body.a.map.θ == body.b.map.θ  == T(π/4))
+        @test all(body.a.map.V .≈ body.b.map.V  .≈ [1,0])
+        # try measure in the sim using different backends
+        for array in arrays
+            body = AutoBody((x,t)->sqrt(sum(abs2,x))-4,RigidMap(SA{T}[16,16,16],SA{T}[0,0,0];
+                             V=SA{T}[0,0,0],ω=SA{T}[0,-0.1,0.1]))
+            sim = Simulation((32,32,32),(1,0,0),8;body,T,mem=array)
+            @test GPUArrays.@allowscalar all(extrema(sim.flow.V) .≈ (-0.9,0.9))
+            sim.body = setmap(sim.body;x₀=SA{T}[16,16,12])
+            @test GPUArrays.@allowscalar all(sim.flow.μ₀[17,17,17,:] .≈ 0)
+        end
+    end
+    rmap = RigidMap(SA[0.,0.],π/4)
+    body = AutoBody((x,t)->√(x'x)-1,rmap)-AutoBody((x,t)->√(x'x)-0.5,rmap) # annulus
+    @test all(measure(setmap(body,ω=1.),SA[0.25,0.],0) .≈ (0.25,SA[-1,0],SA[0,0.25]))
 end