From 30f20f396ccf6b7b97ba7c80614b15cd8201f425 Mon Sep 17 00:00:00 2001
From: Francesc Verdugo <f.verdugo.rojano@vu.nl>
Date: Wed, 13 Nov 2024 13:32:47 +0100
Subject: [PATCH] Linear elasticity example

---
 docs/Project.toml                  |  1 +
 docs/src/examples/problem_types.jl | 58 ++++++++++++++++++++++++++----
 ext/GalerkinToolkitMakieExt.jl     |  2 +-
 src/interpolation.jl               |  4 +--
 4 files changed, 56 insertions(+), 9 deletions(-)

diff --git a/docs/Project.toml b/docs/Project.toml
index 28213c2..4803729 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -10,3 +10,4 @@ Metis = "2679e427-3c69-5b7f-982b-ece356f1e94b"
 PartitionedArrays = "5a9dfac6-5c52-46f7-8278-5e2210713be9"
 PartitionedSolvers = "11b65f7f-80ac-401b-9ef2-3db765482d62"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
+Tensors = "48a634ad-e948-5137-8d70-aa71f2a747f4"
diff --git a/docs/src/examples/problem_types.jl b/docs/src/examples/problem_types.jl
index ebc4f16..731975d 100644
--- a/docs/src/examples/problem_types.jl
+++ b/docs/src/examples/problem_types.jl
@@ -3,17 +3,16 @@
 
 import GalerkinToolkit as GT
 import PartitionedSolvers as PS
-import ForwardDiff
 import GLMakie as Makie
+import ForwardDiff
+import StaticArrays
+import Tensors
 using LinearAlgebra
-using StaticArrays
 import FileIO # hide
 
 
 
 
-
-
 # ## Poisson
 #
 # Solve the following Poisson equation on the unit square,
@@ -153,6 +152,53 @@ nothing # hide
 #     (except for periodic boundary conditions)
 #     Help wanted.
 #
+
+# ## Linear elasticity
+#
+# !!! warning
+#     TODO
+#     * Problem statement
+#     * Hide call
+#     * Use Voigt notation instead of using Tensors.jl?
+#
+
+domain = (0,1,0,0.2)
+cells = (20,5)
+D = length(cells)
+mesh = GT.cartesian_mesh(domain,cells)
+Ω = GT.interior(mesh)
+Γ = GT.boundary(mesh;physical_names=["1-face-3"])
+f = GT.analytical_field(x->StaticArrays.SVector(0,-1),Ω)
+const E = 1
+const ν = 0.33
+const λ = (E*ν)/((1+ν)*(1-2*ν))
+const μ = E/(2*(1+ν))
+σ(ε) = λ*tr(ε)*one(ε) + 2*μ*ε
+order = 2
+V = GT.lagrange_space(Ω,order;shape=(D,),dirichlet_boundary=Γ)
+uhd = GT.dirichlet_field(Float64,V)
+dΩ = GT.measure(Ω,2*order)
+∇ = ForwardDiff.jacobian
+#TODO this function should be better in Tensors.jl
+function symmetrize(m::StaticArrays.SMatrix{2,2})
+    T = eltype(m)
+    Tensors.SymmetricTensor{2,2,T}((m[1,1],0.5*(m[2,1]+m[1,2]),m[2,2]))
+end
+#TODO hide GT.call
+ε(u,x) = GT.call(symmetrize,∇(u,x))
+a(u,v) = GT.∫(x-> GT.call(Tensors.:⊡,ε(v,x),GT.call(σ,ε(u,x))), dΩ)
+l(v) = GT.∫(x-> v(x)⋅f(x), dΩ)
+p = GT.linear_problem(uhd,a,l)
+s = PS.LinearAlgebra_lu(p)
+s = PS.solve(s)
+@show maximum(abs,PS.solution(s))
+uh = GT.solution_field(uhd,s)
+Makie.plot(Ω;color=x->norm(uh(x)),warp_by_vector=uh,warp_scale=0.002)
+Makie.plot!(Ω;color=nothing,strokecolor=:black)
+FileIO.save(joinpath(@__DIR__,"fig_pt_linelast.png"),Makie.current_figure()) # hide
+
+# ![](fig_pt_linelast.png)
+
 # ## Stokes equation
 #
 # !!! warning
@@ -174,8 +220,8 @@ mesh = GT.cartesian_mesh(domain,cells)
 #g2 = GT.analytical_field(x->SVector(0,0),Γ2)
 #g = GT.piecewiese_field(g1,g2)
 #Γ = GT.domain(g)
-g1 = GT.analytical_field(x->SVector(1,0),Ω)
-g2 = GT.analytical_field(x->SVector(0,0),Ω)
+g1 = GT.analytical_field(x->StaticArrays.SVector(1,0),Ω)
+g2 = GT.analytical_field(x->StaticArrays.SVector(0,0),Ω)
 g = GT.piecewiese_field(g1,g2)
 Γ = GT.piecewiese_domain(Γ1,Γ2)
 order = 2
diff --git a/ext/GalerkinToolkitMakieExt.jl b/ext/GalerkinToolkitMakieExt.jl
index eb6fe59..0ee9135 100644
--- a/ext/GalerkinToolkitMakieExt.jl
+++ b/ext/GalerkinToolkitMakieExt.jl
@@ -118,7 +118,7 @@ function Makie.plot!(sc::MakiePlot{<:Tuple{<:GT.AbstractDomain}})
     warp_by_vector = valid_attributes[:warp_by_vector]
     color = valid_attributes[:color]
     args = Makie.lift(dom,color) do dom,color
-        if isa(color,GT.AbstractQuantity)
+        if isa(color,GT.AbstractQuantity) || isa(color,Function)
             label = string(gensym())
             plt = GT.plot(dom)
             GT.plot!(plt,color;label)
diff --git a/src/interpolation.jl b/src/interpolation.jl
index 2d90c96..7d97c1a 100644
--- a/src/interpolation.jl
+++ b/src/interpolation.jl
@@ -183,7 +183,7 @@ function face_dofs_from_nodes(fe::AbstractLagrangeFE,face_to_nodes)
             else
                 error("Not Implemented")
             end
-            reduce(vcat,nested)
+            reduce(vcat,nested;init=Int32[])
         end
     end
 end
@@ -289,7 +289,7 @@ function face_own_dof_permutations(fe::AbstractLagrangeFE,d)
                 else
                     error("Not Implemented")
                 end
-                reduce(vcat,nested)
+                reduce(vcat,nested;init=Tv[])
             end
         end
     end