Merge pull request #108 from bcube-project/dev_rotation

New affine transformation `Rotation`

src/Bcube.jl

@@ -24,7 +24,7 @@ import .LazyOperators:
 include("utils.jl")
 
 include("./mesh/transformation.jl")
-export Translation
+export Translation, Rotation
 
 include("./mesh/boundary_condition.jl")
 export BoundaryCondition, PeriodicBCType

src/mesh/domain.jl

@@ -115,8 +115,8 @@ function _compute_periodicity(mesh, labels1, labels2, A, tol = 1e-9)
     nbnd2 = length(bndfaces2)
 
     # Allocate result
-    bnd_f2n1 = [f2n[iface] for iface in bndfaces2]
-    bnd_f2n2 = [f2n[iface] for iface in bndfaces2]
+    bnd_f2n1 = [zero(f2n[iface]) for iface in bndfaces1]
+    bnd_f2n2 = [zero(f2n[iface]) for iface in bndfaces2]
     bnd_f2c = zeros(Int, nbnd2, 2) # Adjacent cells for each bnd face
     bnd_ftypes = Array{AbstractEntityType}(undef, nbnd2)
 

src/mesh/gmsh_utils.jl

@@ -1240,12 +1240,15 @@ function gen_cylinder_shell_mesh(
     # Mesh settings
     if transfinite
         _nθ = round(Int, nθ / 3)
-        map(
+        foreach(
             line -> gmsh.model.geo.mesh.setTransfiniteCurve(line, _nθ),
             (AOB1, BOC1, COA1, AOB2, BOC2, COA2),
         )
-        map(line -> gmsh.model.geo.mesh.setTransfiniteCurve(line, nz), (A1A2, B1B2, C1C2))
-        map(surf -> gmsh.model.geo.mesh.setTransfiniteSurface(surf), surfs)
+        foreach(
+            line -> gmsh.model.geo.mesh.setTransfiniteCurve(line, nz),
+            (A1A2, B1B2, C1C2),
+        )
+        foreach(gmsh.model.geo.mesh.setTransfiniteSurface, surfs)
     end
     recombine && map(surf -> gmsh.model.geo.mesh.setRecombine(2, surf), surfs)
 

src/mesh/transformation.jl

@@ -7,3 +7,32 @@ struct Translation{T} <: AbstractAffineTransformation
 end
 (t::Translation)(x) = x .+ t.vec
 Base.inv(t::Translation) = Translation(-t.vec)
+
+struct Rotation{T} <: AbstractAffineTransformation
+    A::T
+end
+
+"""
+    Rotation(u::AbstractVector, θ::Number)
+
+Return a `Rotation` operator build from a given
+rotation axis `u` and an angle of rotation `θ`.
+
+# Example :
+```julia
+julia> rot = Rotation([0,0,1], π/4);
+julia> x = [3.0,0,0];
+julia> x_rot_ref = 3√(2)/2 .*[1,1,0] # expected result
+julia> all(rot(x) .≈ x_rot_ref)
+true
+```
+"""
+function Rotation(u::AbstractVector, θ::Number)
+    Rotation((normalize(u), sincos(θ)))
+end
+
+function (r::Rotation)(x)
+    u, (sinθ, cosθ) = r.A
+    rx = cosθ * x + (1 - cosθ) * (u ⋅ x) * u + sinθ * cross(u, x)
+    return rx
+end

test/mesh/test_transformation.jl

@@ -0,0 +1,15 @@
+@testset "Transformation" begin
+    @testset "Translation" begin
+        u = [2, 0]
+        T = Translation(u)
+        x = [3, 1]
+        @test all(T(x) .== [5, 1])
+    end
+
+    @testset "Rotation" begin
+        rot = Rotation([0, 0, 1], π / 4)
+        x = [3.0, 0, 0]
+        x_rot_ref = 3√(2) / 2 .* [1, 1, 0]
+        @test all(rot(x) .≈ x_rot_ref)
+    end
+end

test/runtests.jl

@@ -140,6 +140,7 @@ fname2sum = Dict(r[2] => r[1] for r in eachrow(f))
     custom_include("./test_utils.jl")
     custom_include("./mesh/test_entity.jl")
     custom_include("./mesh/test_connectivity.jl")
+    custom_include("./mesh/test_transformation.jl")
     custom_include("./mesh/test_mesh.jl")
     custom_include("./mesh/test_mesh_generator.jl")
     custom_include("./mesh/test_gmsh.jl")