start linear elasticity tutorial

docs/src/literate-tutorials/linear_elasticity.jl

@@ -41,6 +41,7 @@
 # \int_{\partial\Omega}
 #   \boldsymbol{t}^\ast \cdot \delta \boldsymbol{u}
 # \, \mathrm{d}A
+# +
 # \int_\Omega
 #   \boldsymbol{b} \cdot \delta \boldsymbol{u}
 # \, \mathrm{d}V
@@ -50,6 +51,42 @@
 # $\mathbb{T}$ are suitable trial and test function sets, respectively.
 #-
 
+# First we load Ferrite, and some other packages we need
+using Ferrite, FerriteGmsh, SparseArrays
+# Like for the Heat Equation example, we will use a unit square - but here we'll load the grid of the Ferrite logo! This is done by `FerriteGmsh` here. 
+grid = togrid("logo.geo") # TODO: where to store correctly and how to refer to this?
+# By default the grid lacks the facesets for the boundaries, so we add them by Ferrite here.
+# Note that approximate comparison to 0.0 doesn't work well, so we use a tolerance instead.
+addfaceset!(grid, "right", x->x[1] ≈ 1.0)
+addfaceset!(grid, "top", x->x[2] ≈ 1.0)
+addfaceset!(grid, "left", x->x[1] < 1e-6) 
+addfaceset!(grid, "bottom", x->x[2] < 1e-6)
+
+# ### Trial and test functions
+# We use linear Lagrange functions as test and trial functions. The grid is composed of triangular elements, thus we need the Lagrange functions defined on `RefTriangle`. All currently available interpolations can be found under [`Interpolation`](@ref).
+# Since the displacement field $\boldsymbol{u}$ is a vector valued field, we vectorize the interpolation by the spatial dimension `dim`.
+dim = 2
+order = 1 # linear interpolation
+ip = Lagrange{RefTriangle, order}()^dim
+qr = QuadratureRule{RefQuadrilateral}(2) # 2 quadrature points per spatial dimension
+cellvalues = CellValues(qr, ip)
+
+# ### Degrees of freedom
+# For distributing degrees of freedom, we define a `DofHandler`. The `DofHandler` knows that `u` has two degrees of freedom per node because we vectorized the interpolation above.
+dh = DofHandler(grid)
+add!(dh, :u, ip)
+close!(dh);
+
+# ### Boundary conditions
+# Now, we add boundary conditions. We simply support the bottom and the left side and prescribe a displacement upwards on the top edge.
+# The last argument to `Dirichlet` determines which components of the field should be constrained.
+ch = ConstraintHandler(dh)
+add!(ch, Dirichlet(:u, getfaceset(grid, "bottom"), (x, t) -> 0.0, 2))
+add!(ch, Dirichlet(:u, getfaceset(grid, "left"), (x, t) -> 0.0, 1))
+add!(ch, Dirichlet(:u, getfaceset(grid, "top"), (x, t) -> 0.1, 2))
+close!(ch)
+
+# ### Element routine
 function assemble_cell!(ke, re, cellvalues, cell, material, ue)
     fill!(ke, 0.0)
     fill!(re, 0.0)
@@ -64,15 +101,14 @@ function assemble_cell!(ke, re, cellvalues, cell, material, ue)
 
         dΩ = getdetJdV(cellvalues, q_point)
         for i in 1:n_basefuncs
-            Nᵢ∇ = shape_gradient(cellvalues, q_point, i)# shape_symmetric_gradient(cellvalues, q_point, i)
-            re[i] += σ ⊡ Nᵢ∇ * dΩ # add internal force to residual
-            for j in 1:i # loop only over lower half
+            ∇Nᵢ = shape_gradient(cellvalues, q_point, i)# shape_symmetric_gradient(cellvalues, q_point, i)
+            re[i] += σ ⊡ ∇Nᵢ * dΩ # add internal force to residual
+            for j in 1:n_basefuncs
                 ∇ˢʸᵐNⱼ = shape_symmetric_gradient(cellvalues, q_point, j)
-                ke[i, j] += Nᵢ∇ ⊡ ∂σ∂ε ⊡ ∇ˢʸᵐNⱼ * dΩ
+                ke[i, j] += ∇Nᵢ ⊡ ∂σ∂ε ⊡ ∇ˢʸᵐNⱼ * dΩ
             end
         end
     end
-    # symmetrize_lower!(ke) # needed? what does assembly in symmetric global matrix does?
 end
 
 

docs/src/literate-tutorials/logo.geo

@@ -0,0 +1,43 @@
+Point (1) = {1.000000000000,1.000000000000,0.000000000000};
+Point (2) = {1.000000000000,0.379757979263,-0.000000000000};
+Point (3) = {0.801534880751,0.454963545532,-0.000000000000};
+Point (4) = {0.656107331955,1.000000000000,-0.000000000000};
+Point (5) = {0.600672553037,0.775245709538,-0.000000000000};
+Point (6) = {0.000000000000,1.000000000000,0.000000000000};
+Point (7) = {0.392825178821,0.672136259831,-0.000000000000};
+Point (8) = {1.000000000000,0.000000000000,0.000000000000};
+Point (9) = {0.547800422194,-0.000000000000,-0.000000000000};
+Point (10) = {0.488710023938,0.224380304618,-0.000000000000};
+Point (11) = {0.000000000000,0.000000000000,0.000000000000};
+Point (12) = {-0.000000000000,0.324566579562,-0.000000000000};
+Point (13) = {0.172066723668,0.367888021869,-0.000000000000};
+Line (1) = {2,1};
+Line (2) = {1,4};
+Line (3) = {2,3};
+Line (4) = {3,5};
+Line (5) = {5,4};
+Line (6) = {4,6};
+Line (7) = {7,6};
+Line (8) = {5,7};
+Line (9) = {8,2};
+Line (10) = {9,8};
+Line (11) = {9,10};
+Line (12) = {10,3};
+Line (13) = {12,11};
+Line (14) = {11,9};
+Line (15) = {12,13};
+Line (16) = {13,10};
+Line (17) = {6,12};
+Line (18) = {7,13};
+Line Loop (1) = {-1,3,4,5,-2};
+Plane Surface (1) = {1}; Physical Surface (1) = {1};
+Line Loop (2) = {7,-6,-5,8};
+Plane Surface (2) = {2}; Physical Surface (2) = {2};
+Line Loop (3) = {-10,11,12,-3,-9};
+Plane Surface (3) = {3}; Physical Surface (3) = {3};
+Line Loop (4) = {-11,-14,-13,15,16};
+Plane Surface (4) = {4}; Physical Surface (4) = {4};
+Line Loop (5) = {-7,18,-15,-17};
+Plane Surface (5) = {5}; Physical Surface (5) = {5};
+Line Loop (6) = {-16,-18,-8,-4,-12};
+Plane Surface (6) = {6}; Physical Surface (6) = {6};