Merge pull request #141 from GalerkinToolkit/docs

More examples (Stokes and linear elasticity)

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"

docs/src/examples/problem_types.jl

@@ -3,11 +3,16 @@
 
 import GalerkinToolkit as GT
 import PartitionedSolvers as PS
-import ForwardDiff
 import GLMakie as Makie
+import ForwardDiff
+import StaticArrays
+import Tensors
 using LinearAlgebra
 import FileIO # hide
 
+
+
+
 # ## Poisson
 #
 # Solve the following Poisson equation on the unit square,
@@ -147,10 +152,104 @@ nothing # hide
 #     (except for periodic boundary conditions)
 #     Help wanted.
 #
-# ## Stokes equation
+
+# ## Linear elasticity
 #
 # !!! warning
-#     TODO Lid cavity problem. This will illustrate how to use vector-valued spaces and multifield.
-#     The key ingredient missing is to develop a strategy to impose a zero mean constrain on the pressure.
+#     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
+#     TODO
+#     * Problem statement
+#     * Allow g1 and g1 to be defined on the boundary
+#     * Build a pressure field with zero mean
 #
+
+domain = (0,1,0,1)
+cells = (20,20)
+D = length(cells)
+mesh = GT.cartesian_mesh(domain,cells)
+Ω = GT.interior(mesh)
+Γ1 = GT.boundary(mesh;physical_names=["1-face-2"])
+Γ2 = GT.boundary(mesh;physical_names=["1-face-1","1-face-3","1-face-4"])
+#TODO
+#g1 = GT.analytical_field(x->SVector(1,0),Γ1)
+#g2 = GT.analytical_field(x->SVector(0,0),Γ2)
+#g = GT.piecewiese_field(g1,g2)
+#Γ = GT.domain(g)
+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
+V = GT.lagrange_space(Ω,order;space=:Q,shape=(D,),dirichlet_boundary=Γ)
+Q = GT.lagrange_space(Ω,order-1;space=:P,dirichlet_boundary=GT.last_dof())
+VxQ = V × Q
+u_field, p_field = 1,2
+uhph_dirichlet = GT.dirichlet_field(Float64,VxQ)
+GT.interpolate_dirichlet!(g,uhph_dirichlet,u_field)
+dΩ = GT.measure(Ω,2*order)
+∇ = ForwardDiff.jacobian
+div(u,x) = tr(∇(u,x))
+a((u,p),(v,q)) = GT.∫( x-> ∇(v,x)⋅∇(u,x) - div(v,x)*p(x) + q(x)*div(u,x), dΩ)
+l((v,q)) = 0
+p = GT.linear_problem(uhph_dirichlet,a,l)
+s = PS.LinearAlgebra_lu(p)
+s = PS.solve(s)
+#TODO
+#uh = GT.solution_field(uhph_dirichlet,s,u_field)
+#ph = GT.solution_field(uhph_dirichlet,s,p_field;zeromean=true)
+uh,ph = GT.solution_field(uhph_dirichlet,s)
+Makie.plot(Ω,color=ph)
+Makie.arrows!(uh;color=x->norm(uh(x)),lengthscale=0.1)
+FileIO.save(joinpath(@__DIR__,"fig_pt_stokes.png"),Makie.current_figure()) # hide
+
+# ![](fig_pt_stokes.png)
+
+
+
+
+