Linear Elasticity
Figure 1: Linear elastically deformed 1mm $\times$ 1mm Ferrite logo.
The full explanation for the underlying FEM theory in this example can be found in the Linear Elasticity tutorial of the Ferrite.jl documentation.
Implementation
The following code is based on the Linear Elasticity tutorial from the Ferrite.jl documentation, with some comments removed for brevity. There are two main modifications:
- Second-order
Lagrangeshape functions are used for field approximation:ip = Lagrange{RefTriangle,2}()^2. - Four quadrature points are used to accommodate the second-order shape functions:
qr = QuadratureRule{RefTriangle}(4).
using Ferrite, FerriteGmsh, SparseArrays
using Downloads: download
Emod = 200.0e3 # Young's modulus [MPa]
ν = 0.3 # Poisson's ratio [-]
Gmod = Emod / (2(1 + ν)) # Shear modulus
Kmod = Emod / (3(1 - 2ν)) # Bulk modulus
C = gradient(ϵ -> 2 * Gmod * dev(ϵ) + 3 * Kmod * vol(ϵ), zero(SymmetricTensor{2,2}))
function assemble_external_forces!(f_ext, dh, facetset, facetvalues, prescribed_traction)
# Create a temporary array for the facet's local contributions to the external force vector
fe_ext = zeros(getnbasefunctions(facetvalues))
for facet in FacetIterator(dh, facetset)
# Update the facetvalues to the correct facet number
reinit!(facetvalues, facet)
# Reset the temporary array for the next facet
fill!(fe_ext, 0.0)
# Access the cell's coordinates
cell_coordinates = getcoordinates(facet)
for qp in 1:getnquadpoints(facetvalues)
# Calculate the global coordinate of the quadrature point.
x = spatial_coordinate(facetvalues, qp, cell_coordinates)
tₚ = prescribed_traction(x)
# Get the integration weight for the current quadrature point.
dΓ = getdetJdV(facetvalues, qp)
for i in 1:getnbasefunctions(facetvalues)
Nᵢ = shape_value(facetvalues, qp, i)
fe_ext[i] += tₚ ⋅ Nᵢ * dΓ
end
end
# Add the local contributions to the correct indices in the global external force vector
assemble!(f_ext, celldofs(facet), fe_ext)
end
return f_ext
end
function assemble_cell!(ke, cellvalues, C)
for q_point in 1:getnquadpoints(cellvalues)
# Get the integration weight for the quadrature point
dΩ = getdetJdV(cellvalues, q_point)
for i in 1:getnbasefunctions(cellvalues)
# Gradient of the test function
∇Nᵢ = shape_gradient(cellvalues, q_point, i)
for j in 1:getnbasefunctions(cellvalues)
# Symmetric gradient of the trial function
∇ˢʸᵐNⱼ = shape_symmetric_gradient(cellvalues, q_point, j)
ke[i, j] += (∇Nᵢ ⊡ C ⊡ ∇ˢʸᵐNⱼ) * dΩ
end
end
end
return ke
end
function assemble_global!(K, dh, cellvalues, C)
# Allocate the element stiffness matrix
n_basefuncs = getnbasefunctions(cellvalues)
ke = zeros(n_basefuncs, n_basefuncs)
# Create an assembler
assembler = start_assemble(K)
# Loop over all cells
for cell in CellIterator(dh)
# Update the shape function gradients based on the cell coordinates
reinit!(cellvalues, cell)
# Reset the element stiffness matrix
fill!(ke, 0.0)
# Compute element contribution
assemble_cell!(ke, cellvalues, C)
# Assemble ke into K
assemble!(assembler, celldofs(cell), ke)
end
return K
end
function linear_elasticity_2d(C)
logo_mesh = "logo.geo"
asset_url = "https://raw.githubusercontent.com/Ferrite-FEM/Ferrite.jl/gh-pages/assets/"
isfile(logo_mesh) || download(string(asset_url, logo_mesh), logo_mesh)
grid = togrid(logo_mesh)
addfacetset!(grid, "top", x -> x[2] ≈ 1.0) # facets for which x[2] ≈ 1.0 for all nodes
addfacetset!(grid, "left", x -> abs(x[1]) < 1.0e-6)
addfacetset!(grid, "bottom", x -> abs(x[2]) < 1.0e-6)
dim = 2
order = 2 # quadratic interpolation
ip = Lagrange{RefTriangle,order}()^dim # vector valued interpolation
qr = QuadratureRule{RefTriangle}(4) # 4 quadrature point
qr_face = FacetQuadratureRule{RefTriangle}(1)
cellvalues = CellValues(qr, ip)
facetvalues = FacetValues(qr_face, ip)
dh = DofHandler(grid)
add!(dh, :u, ip)
close!(dh)
ch = ConstraintHandler(dh)
add!(ch, Dirichlet(:u, getfacetset(grid, "bottom"), (x, t) -> 0.0, 2))
add!(ch, Dirichlet(:u, getfacetset(grid, "left"), (x, t) -> 0.0, 1))
close!(ch)
traction(x) = Vec(0.0, 20.0e3 * x[1])
A = allocate_matrix(dh)
assemble_global!(A, dh, cellvalues, C)
b = zeros(ndofs(dh))
assemble_external_forces!(b, dh, getfacetset(grid, "top"), facetvalues, traction)
apply!(A, b, ch)
return A, b, dh, cellvalues, ch
endlinear_elasticity_2d (generic function with 1 method)Near Null Space (NNS)
In multigrid methods for problems with vector-valued unknowns, such as linear elasticity, the near null space represents the low energy mode or the smooth error that needs to be captured in the coarser grid when using SA-AMG (Smoothed Aggregation Algebraic Multigrid), more on the topic can be found in Schroder [1].
For 2D linear elasticity problems, the rigid body modes are:
- Translation in the x-direction,
- Translation in the y-direction,
- Rotation about the z-axis (i.e., $x_3$): each point (x, y) is mapped to (-y, x).
The function create_nns constructs the NNS matrix B ∈ ℝ^{n × 3}, where n is the number of degrees of freedom (DOFs) for the case of p = 1 (i.e., linear interpolation), because B is only relevant for AMG.
function create_nns(dh)
##Ndof = ndofs(dh)
grid = dh.grid
Ndof = 2 * (grid.nodes |> length) # nns at p = 1 for AMG
B = zeros(Float64, Ndof, 3)
B[1:2:end, 1] .= 1 # x - translation
B[2:2:end, 2] .= 1 # y - translation
# in-plane rotation (x,y) → (-y,x)
coords = reduce(hcat, grid.nodes .|> (n -> n.x |> collect))' # convert nodes to 2d array
y = coords[:, 2]
x = coords[:, 1]
B[1:2:end, 3] .= -y
B[2:2:end, 3] .= x
return B
endcreate_nns (generic function with 1 method)Setup the linear elasticity problem
Load FerriteMultigrid to access the p-multigrid solver.
using FerriteMultigridConstruct the linear elasticity problem with 2nd order polynomial shape functions.
A, b, dh, cellvalues, ch = linear_elasticity_2d(C);Info : Reading 'logo.geo'...
Info : Done reading 'logo.geo'
Info : Meshing 1D...
Info : [ 0%] Meshing curve 1 (Line)
Info : [ 10%] Meshing curve 2 (Line)
Info : [ 20%] Meshing curve 3 (Line)
Info : [ 20%] Meshing curve 4 (Line)
Info : [ 30%] Meshing curve 5 (Line)
Info : [ 30%] Meshing curve 6 (Line)
Info : [ 40%] Meshing curve 7 (Line)
Info : [ 40%] Meshing curve 8 (Line)
Info : [ 50%] Meshing curve 9 (Line)
Info : [ 60%] Meshing curve 10 (Line)
Info : [ 60%] Meshing curve 11 (Line)
Info : [ 70%] Meshing curve 12 (Line)
Info : [ 70%] Meshing curve 13 (Line)
Info : [ 80%] Meshing curve 14 (Line)
Info : [ 80%] Meshing curve 15 (Line)
Info : [ 90%] Meshing curve 16 (Line)
Info : [ 90%] Meshing curve 17 (Line)
Info : [100%] Meshing curve 18 (Line)
Info : Done meshing 1D (Wall 0.0010289s, CPU 0.001029s)
Info : Meshing 2D...
Info : [ 0%] Meshing surface 1 (Plane, Frontal-Delaunay)
Info : [ 20%] Meshing surface 2 (Plane, Frontal-Delaunay)
Info : [ 40%] Meshing surface 3 (Plane, Frontal-Delaunay)
Info : [ 60%] Meshing surface 4 (Plane, Frontal-Delaunay)
Info : [ 70%] Meshing surface 5 (Plane, Frontal-Delaunay)
Info : [ 90%] Meshing surface 6 (Plane, Frontal-Delaunay)
Info : Done meshing 2D (Wall 0.00191703s, CPU 0.001917s)
Info : 104 nodes 245 elementsConstruct the near null space (NNS) matrix
B = create_nns(dh)208×3 Matrix{Float64}:
1.0 0.0 -1.0
0.0 1.0 1.0
1.0 0.0 -0.379758
0.0 1.0 1.0
1.0 0.0 -0.454964
0.0 1.0 0.801535
1.0 0.0 -1.0
0.0 1.0 0.656107
1.0 0.0 -0.775246
0.0 1.0 0.600673
⋮
0.0 1.0 0.567386
1.0 0.0 -0.487217
0.0 1.0 0.365932
1.0 0.0 -0.477795
0.0 1.0 0.677392
1.0 0.0 -0.534479
0.0 1.0 0.52256
1.0 0.0 -0.640261
0.0 1.0 0.473388Since NNS matrix is only relevant for AMG, and it is not used in the p-multigrid solver, therefore, B has to provided using linear field approximation (i.e., p = 1) when using AMG as the coarse solver, otherwise (e.g., using Pinv as the coarse solver), then we don't have to provide it.
Construct the finite element space $\mathcal{V}_{h,p = 2}$
fe_space = FESpace(dh, cellvalues, ch)FESpace{Ferrite.DofHandler{2, Ferrite.Grid{2, Ferrite.Triangle, Float64}}, Ferrite.CellValues{Ferrite.FunctionValues{1, Ferrite.VectorizedInterpolation{2, Ferrite.RefTriangle, 2, Ferrite.Lagrange{Ferrite.RefTriangle, 2}}, Matrix{Tensors.Vec{2, Float64}}, Matrix{Tensors.Tensor{2, 2, Float64, 4}}, Matrix{Tensors.Tensor{2, 2, Float64, 4}}, Nothing, Nothing}, Ferrite.GeometryMapping{1, Ferrite.Lagrange{Ferrite.RefTriangle, 1}, Matrix{Float64}, Matrix{Tensors.Vec{2, Float64}}, Nothing}, Ferrite.QuadratureRule{Ferrite.RefTriangle, Vector{Float64}, Vector{Tensors.Vec{2, Float64}}}, Vector{Float64}}, Ferrite.ConstraintHandler{Ferrite.DofHandler{2, Ferrite.Grid{2, Ferrite.Triangle, Float64}}, Float64}}(Ferrite.DofHandler{2, Ferrite.Grid{2, Ferrite.Triangle, Float64}}(Ferrite.SubDofHandler{Ferrite.DofHandler{2, Ferrite.Grid{2, Ferrite.Triangle, Float64}}}[Ferrite.SubDofHandler{Ferrite.DofHandler{2, Ferrite.Grid{2, Ferrite.Triangle, Float64}}}(Ferrite.DofHandler{2, Ferrite.Grid{2, Ferrite.Triangle, Float64}}(#= circular reference @-3 =#), OrderedCollections.OrderedSet{Int64}([1, 2, 3, 4, 5, 6, 7, 8, 9, 10 … 165, 166, 167, 168, 169, 170, 171, 172, 173, 174]), [:u], Ferrite.Interpolation[Ferrite.VectorizedInterpolation{2, Ferrite.RefTriangle, 2, Ferrite.Lagrange{Ferrite.RefTriangle, 2}}(Ferrite.Lagrange{Ferrite.RefTriangle, 2}())], Int64[], 12)], [:u], [1, 2, 3, 4, 5, 6, 7, 8, 9, 10 … 241, 242, 681, 682, 705, 706, 755, 756, 685, 686], [1, 13, 25, 37, 49, 61, 73, 85, 97, 109 … 1969, 1981, 1993, 2005, 2017, 2029, 2041, 2053, 2065, 2077], [1, 1, 1, 1, 1, 1, 1, 1, 1, 1 … 1, 1, 1, 1, 1, 1, 1, 1, 1, 1], true, Ferrite.Grid{2, Ferrite.Triangle, Float64}(Ferrite.Triangle[Ferrite.Triangle((3, 20, 63)), Ferrite.Triangle((55, 22, 56)), Ferrite.Triangle((20, 57, 63)), Ferrite.Triangle((20, 14, 57)), Ferrite.Triangle((55, 56, 61)), Ferrite.Triangle((58, 55, 61)), Ferrite.Triangle((55, 58, 60)), Ferrite.Triangle((5, 22, 55)), Ferrite.Triangle((2, 14, 20)), Ferrite.Triangle((57, 56, 63)) … Ferrite.Triangle((97, 22, 100)), Ferrite.Triangle((21, 97, 102)), Ferrite.Triangle((7, 94, 104)), Ferrite.Triangle((94, 52, 101)), Ferrite.Triangle((31, 7, 104)), Ferrite.Triangle((95, 99, 101)), Ferrite.Triangle((99, 94, 101)), Ferrite.Triangle((99, 96, 103)), Ferrite.Triangle((94, 99, 103)), Ferrite.Triangle((100, 31, 104))], Ferrite.Node{2, Float64}[Ferrite.Node{2, Float64}([1.0, 1.0]), Ferrite.Node{2, Float64}([1.0, 0.379757979263]), Ferrite.Node{2, Float64}([0.801534880751, 0.454963545532]), Ferrite.Node{2, Float64}([0.656107331955, 1.0]), Ferrite.Node{2, Float64}([0.600672553037, 0.775245709538]), Ferrite.Node{2, Float64}([0.0, 1.0]), Ferrite.Node{2, Float64}([0.392825178821, 0.672136259831]), Ferrite.Node{2, Float64}([1.0, 0.0]), Ferrite.Node{2, Float64}([0.547800422194, -0.0]), Ferrite.Node{2, Float64}([0.488710023938, 0.224380304618]) … Ferrite.Node{2, Float64}([0.371815483395454, 0.387541210567491]), Ferrite.Node{2, Float64}([0.5769806418976112, 0.4322841106613628]), Ferrite.Node{2, Float64}([0.6165746595005666, 0.5544955814857399]), Ferrite.Node{2, Float64}([0.47736325015237235, 0.33970530481730987]), Ferrite.Node{2, Float64}([0.46406870923827526, 0.45073264476884783]), Ferrite.Node{2, Float64}([0.5673862832010109, 0.6502733177441277]), Ferrite.Node{2, Float64}([0.3659322703052846, 0.48721715390516906]), Ferrite.Node{2, Float64}([0.67739173818206, 0.47779490322723067]), Ferrite.Node{2, Float64}([0.5225603855864607, 0.5344788107726228]), Ferrite.Node{2, Float64}([0.473387537304908, 0.6402608858085923])], Dict{String, OrderedCollections.OrderedSet{Int64}}("4" => OrderedCollections.OrderedSet{Int64}([114, 105, 110, 111, 92, 99, 98, 89, 106, 96 … 102, 93, 113, 97, 88, 104, 91, 87, 100, 101]), "1" => OrderedCollections.OrderedSet{Int64}([5, 16, 20, 12, 24, 28, 8, 17, 30, 1 … 33, 4, 13, 15, 2, 10, 18, 21, 26, 27]), "5" => OrderedCollections.OrderedSet{Int64}([140, 123, 122, 133, 137, 136, 117, 135, 121, 115 … 139, 128, 124, 129, 131, 120, 127, 116, 132, 134]), "2" => OrderedCollections.OrderedSet{Int64}([35, 52, 37, 53, 47, 41, 43, 45, 36, 44 … 39, 51, 46, 40, 48, 34, 50, 54, 38, 42]), "6" => OrderedCollections.OrderedSet{Int64}([158, 169, 173, 141, 148, 160, 154, 171, 145, 146 … 168, 174, 163, 147, 165, 144, 142, 167, 157, 150]), "3" => OrderedCollections.OrderedSet{Int64}([78, 56, 55, 79, 81, 58, 60, 72, 75, 83 … 71, 66, 76, 57, 59, 70, 65, 63, 86, 62])), Dict{String, OrderedCollections.OrderedSet{Int64}}(), Dict{String, OrderedCollections.OrderedSet{Ferrite.FacetIndex}}("left" => OrderedCollections.OrderedSet{Ferrite.FacetIndex}([Ferrite.FacetIndex((96, 1)), Ferrite.FacetIndex((101, 1)), Ferrite.FacetIndex((108, 1)), Ferrite.FacetIndex((115, 2)), Ferrite.FacetIndex((118, 2)), Ferrite.FacetIndex((122, 1)), Ferrite.FacetIndex((125, 1)), Ferrite.FacetIndex((127, 1))]), "bottom" => OrderedCollections.OrderedSet{Ferrite.FacetIndex}([Ferrite.FacetIndex((69, 1)), Ferrite.FacetIndex((73, 1)), Ferrite.FacetIndex((78, 1)), Ferrite.FacetIndex((79, 1)), Ferrite.FacetIndex((90, 3)), Ferrite.FacetIndex((104, 1)), Ferrite.FacetIndex((106, 1)), Ferrite.FacetIndex((107, 1))]), "top" => OrderedCollections.OrderedSet{Ferrite.FacetIndex}([Ferrite.FacetIndex((21, 1)), Ferrite.FacetIndex((23, 1)), Ferrite.FacetIndex((28, 1)), Ferrite.FacetIndex((35, 1)), Ferrite.FacetIndex((36, 2)), Ferrite.FacetIndex((42, 3)), Ferrite.FacetIndex((43, 3)), Ferrite.FacetIndex((47, 3))])), Dict{String, OrderedCollections.OrderedSet{Ferrite.VertexIndex}}()), 762), Ferrite.CellValues{Ferrite.FunctionValues{1, Ferrite.VectorizedInterpolation{2, Ferrite.RefTriangle, 2, Ferrite.Lagrange{Ferrite.RefTriangle, 2}}, Matrix{Tensors.Vec{2, Float64}}, Matrix{Tensors.Tensor{2, 2, Float64, 4}}, Matrix{Tensors.Tensor{2, 2, Float64, 4}}, Nothing, Nothing}, Ferrite.GeometryMapping{1, Ferrite.Lagrange{Ferrite.RefTriangle, 1}, Matrix{Float64}, Matrix{Tensors.Vec{2, Float64}}, Nothing}, Ferrite.QuadratureRule{Ferrite.RefTriangle, Vector{Float64}, Vector{Tensors.Vec{2, Float64}}}, Vector{Float64}}(Ferrite.FunctionValues{1, Ferrite.VectorizedInterpolation{2, Ferrite.RefTriangle, 2, Ferrite.Lagrange{Ferrite.RefTriangle, 2}}, Matrix{Tensors.Vec{2, Float64}}, Matrix{Tensors.Tensor{2, 2, Float64, 4}}, Matrix{Tensors.Tensor{2, 2, Float64, 4}}, Nothing, Nothing}(Ferrite.VectorizedInterpolation{2, Ferrite.RefTriangle, 2, Ferrite.Lagrange{Ferrite.RefTriangle, 2}}(Ferrite.Lagrange{Ferrite.RefTriangle, 2}()), Tensors.Vec{2, Float64}[[-0.04820837781550805, 0.0] [-0.04820837781550805, 0.0] … [-0.07480380774819603, 0.0] [0.5176323419876757, 0.0]; [0.0, -0.04820837781550805] [0.0, -0.04820837781550805] … [0.0, -0.07480380774819603] [0.0, 0.5176323419876757]; … ; [0.1928335112620322, 0.0] [0.7954802262009061, 0.0] … [0.03354481152314795, 0.0] [0.2992152309927843, 0.0]; [0.0, 0.1928335112620322] [0.0, 0.7954802262009061] … [0.0, 0.03354481152314795] [0.0, 0.2992152309927843]], Tensors.Vec{2, Float64}[[-0.04820837781550805, 0.0] [-0.04820837781550805, 0.0] … [-0.07480380774819603, 0.0] [0.5176323419876757, 0.0]; [0.0, -0.04820837781550805] [0.0, -0.04820837781550805] … [0.0, -0.07480380774819603] [0.0, 0.5176323419876757]; … ; [0.1928335112620322, 0.0] [0.7954802262009061, 0.0] … [0.03354481152314795, 0.0] [0.2992152309927843, 0.0]; [0.0, 0.1928335112620322] [0.0, 0.7954802262009061] … [0.0, 0.03354481152314795] [0.0, 0.2992152309927843]], Tensors.Tensor{2, 2, Float64, 4}[[8.594688477838039 -2.4066055855083404; 0.0 0.0] [8.594688477838039 -2.4066055855083404; 0.0 0.0] … [-6.948780702511049 1.9457336347043792; 0.0 0.0] [24.863056007279276 -6.961924171447853; 0.0 0.0]; [0.0 0.0; 8.594688477838039 -2.4066055855083404] [0.0 0.0; 8.594688477838039 -2.4066055855083404] … [0.0 0.0; -6.948780702511049 1.9457336347043792] [0.0 0.0; 24.863056007279276 -6.961924171447853]; … ; [-12.471156596172827 -17.888642063253705; 0.0 0.0] [2.3474142350836686 -22.038001917739507; 0.0 0.0] … [0.4820451499822285 -4.5255378368976835; 0.0 0.0] [-27.512058462374867 -31.45953424974912; 0.0 0.0]; [0.0 0.0; -12.471156596172827 -17.888642063253705] [0.0 0.0; 2.3474142350836686 -22.038001917739507] … [0.0 0.0; 0.4820451499822285 -4.5255378368976835] [0.0 0.0; -27.512058462374867 -31.45953424974912]], Tensors.Tensor{2, 2, Float64, 4}[[0.78379396366388 0.0; 0.0 0.0] [0.78379396366388 0.0; 0.0 0.0] … [-0.6336951459609201 0.0; 0.0 0.0] [2.26739029192184 0.0; 0.0 0.0]; [0.0 0.0; 0.78379396366388 0.0] [0.0 0.0; 0.78379396366388 0.0] … [0.0 0.0; -0.6336951459609201 0.0] [0.0 0.0; 2.26739029192184 0.0]; … ; [-1.35138189099164 -1.78379396366388; 0.0 0.0] [-3.9968028886505635e-14 -1.78379396366388; 0.0 0.0] … [3.3306690738754696e-16 -0.36630485403908; 0.0 0.0] [-2.9010854378827595 -3.26739029192184; 0.0 0.0]; [0.0 0.0; -1.35138189099164 -1.78379396366388] [0.0 0.0; -3.9968028886505635e-14 -1.78379396366388] … [0.0 0.0; 3.3306690738754696e-16 -0.36630485403908] [0.0 0.0; -2.9010854378827595 -3.26739029192184]], nothing, nothing), Ferrite.GeometryMapping{1, Ferrite.Lagrange{Ferrite.RefTriangle, 1}, Matrix{Float64}, Matrix{Tensors.Vec{2, Float64}}, Nothing}(Ferrite.Lagrange{Ferrite.RefTriangle, 1}(), [0.44594849091597 0.44594849091597 … 0.09157621350977 0.81684757298046; 0.44594849091597 0.10810301816807 … 0.81684757298046 0.09157621350977; 0.10810301816806 0.44594849091596 … 0.09157621350977008 0.09157621350977005], Tensors.Vec{2, Float64}[[1.0, 0.0] [1.0, 0.0] … [1.0, 0.0] [1.0, 0.0]; [0.0, 1.0] [0.0, 1.0] … [0.0, 1.0] [0.0, 1.0]; [-1.0, -1.0] [-1.0, -1.0] … [-1.0, -1.0] [-1.0, -1.0]], nothing), Ferrite.QuadratureRule{Ferrite.RefTriangle, Vector{Float64}, Vector{Tensors.Vec{2, Float64}}}([0.111690794839005, 0.111690794839005, 0.111690794839005, 0.05497587182766, 0.05497587182766, 0.05497587182766], Tensors.Vec{2, Float64}[[0.44594849091597, 0.44594849091597], [0.44594849091597, 0.10810301816807], [0.10810301816807, 0.44594849091597], [0.09157621350977, 0.09157621350977], [0.09157621350977, 0.81684757298046], [0.81684757298046, 0.09157621350977]]), [0.0008497905835494427, 0.0008497905835494427, 0.0008497905835494427, 0.00041827957504382953, 0.00041827957504382953, 0.00041827957504382953]), Ferrite.ConstraintHandler{Ferrite.DofHandler{2, Ferrite.Grid{2, Ferrite.Triangle, Float64}}, Float64}(Ferrite.Dirichlet[Ferrite.Dirichlet(Main.var"#6#11"(), OrderedCollections.OrderedSet{Ferrite.FacetIndex}(Ferrite.FacetIndex[Ferrite.FacetIndex((69, 1)), Ferrite.FacetIndex((73, 1)), Ferrite.FacetIndex((78, 1)), Ferrite.FacetIndex((79, 1)), Ferrite.FacetIndex((90, 3)), Ferrite.FacetIndex((104, 1)), Ferrite.FacetIndex((106, 1)), Ferrite.FacetIndex((107, 1))]), :u, [2], [2, 4, 8, 4, 6, 10, 6, 2, 12], [1, 4, 7, 10]), Ferrite.Dirichlet(Main.var"#7#12"(), OrderedCollections.OrderedSet{Ferrite.FacetIndex}(Ferrite.FacetIndex[Ferrite.FacetIndex((96, 1)), Ferrite.FacetIndex((101, 1)), Ferrite.FacetIndex((108, 1)), Ferrite.FacetIndex((115, 2)), Ferrite.FacetIndex((118, 2)), Ferrite.FacetIndex((122, 1)), Ferrite.FacetIndex((125, 1)), Ferrite.FacetIndex((127, 1))]), :u, [1], [1, 3, 7, 3, 5, 9, 5, 1, 11], [1, 4, 7, 10])], Ferrite.ProjectedDirichlet[], [177, 354, 356, 358, 370, 374, 382, 392, 400, 404 … 539, 547, 551, 561, 563, 567, 591, 593, 609, 619], [1, 2, 3, 4, 5, 6, 7, 8, 9, 10 … 753, 754, 755, 756, 757, 758, 759, 760, 761, 762], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 … 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], Union{Nothing, Float64}[nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing … nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing], Union{Nothing, Vector{Pair{Int64, Float64}}}[nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing … nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing], Dict(479 => 14, 499 => 17, 392 => 8, 370 => 5, 514 => 20, 551 => 27, 177 => 1, 561 => 28, 563 => 29, 374 => 6…), Ferrite.BCValues{Float64}[Ferrite.BCValues{Float64}([1.0 0.0 0.5; 0.0 1.0 0.5; 0.0 0.0 0.0;;; 0.0 0.0 0.0; 1.0 0.0 0.5; 0.0 1.0 0.5;;; 0.0 1.0 0.5; 0.0 0.0 0.0; 1.0 0.0 0.5], [3, 3, 3], 1), Ferrite.BCValues{Float64}([1.0 0.0 0.5; 0.0 1.0 0.5; 0.0 0.0 0.0;;; 0.0 0.0 0.0; 1.0 0.0 0.5; 0.0 1.0 0.5;;; 0.0 1.0 0.5; 0.0 0.0 0.0; 1.0 0.0 0.5], [3, 3, 3], 1)], Ferrite.DofHandler{2, Ferrite.Grid{2, Ferrite.Triangle, Float64}}(Ferrite.SubDofHandler{Ferrite.DofHandler{2, Ferrite.Grid{2, Ferrite.Triangle, Float64}}}[Ferrite.SubDofHandler{Ferrite.DofHandler{2, Ferrite.Grid{2, Ferrite.Triangle, Float64}}}(Ferrite.DofHandler{2, Ferrite.Grid{2, Ferrite.Triangle, Float64}}(#= circular reference @-3 =#), OrderedCollections.OrderedSet{Int64}([1, 2, 3, 4, 5, 6, 7, 8, 9, 10 … 165, 166, 167, 168, 169, 170, 171, 172, 173, 174]), [:u], Ferrite.Interpolation[Ferrite.VectorizedInterpolation{2, Ferrite.RefTriangle, 2, Ferrite.Lagrange{Ferrite.RefTriangle, 2}}(Ferrite.Lagrange{Ferrite.RefTriangle, 2}())], Int64[], 12)], [:u], [1, 2, 3, 4, 5, 6, 7, 8, 9, 10 … 241, 242, 681, 682, 705, 706, 755, 756, 685, 686], [1, 13, 25, 37, 49, 61, 73, 85, 97, 109 … 1969, 1981, 1993, 2005, 2017, 2029, 2041, 2053, 2065, 2077], [1, 1, 1, 1, 1, 1, 1, 1, 1, 1 … 1, 1, 1, 1, 1, 1, 1, 1, 1, 1], true, Ferrite.Grid{2, Ferrite.Triangle, Float64}(Ferrite.Triangle[Ferrite.Triangle((3, 20, 63)), Ferrite.Triangle((55, 22, 56)), Ferrite.Triangle((20, 57, 63)), Ferrite.Triangle((20, 14, 57)), Ferrite.Triangle((55, 56, 61)), Ferrite.Triangle((58, 55, 61)), Ferrite.Triangle((55, 58, 60)), Ferrite.Triangle((5, 22, 55)), Ferrite.Triangle((2, 14, 20)), Ferrite.Triangle((57, 56, 63)) … Ferrite.Triangle((97, 22, 100)), Ferrite.Triangle((21, 97, 102)), Ferrite.Triangle((7, 94, 104)), Ferrite.Triangle((94, 52, 101)), Ferrite.Triangle((31, 7, 104)), Ferrite.Triangle((95, 99, 101)), Ferrite.Triangle((99, 94, 101)), Ferrite.Triangle((99, 96, 103)), Ferrite.Triangle((94, 99, 103)), Ferrite.Triangle((100, 31, 104))], Ferrite.Node{2, Float64}[Ferrite.Node{2, Float64}([1.0, 1.0]), Ferrite.Node{2, Float64}([1.0, 0.379757979263]), Ferrite.Node{2, Float64}([0.801534880751, 0.454963545532]), Ferrite.Node{2, Float64}([0.656107331955, 1.0]), Ferrite.Node{2, Float64}([0.600672553037, 0.775245709538]), Ferrite.Node{2, Float64}([0.0, 1.0]), Ferrite.Node{2, Float64}([0.392825178821, 0.672136259831]), Ferrite.Node{2, Float64}([1.0, 0.0]), Ferrite.Node{2, Float64}([0.547800422194, -0.0]), Ferrite.Node{2, Float64}([0.488710023938, 0.224380304618]) … Ferrite.Node{2, Float64}([0.371815483395454, 0.387541210567491]), Ferrite.Node{2, Float64}([0.5769806418976112, 0.4322841106613628]), Ferrite.Node{2, Float64}([0.6165746595005666, 0.5544955814857399]), Ferrite.Node{2, Float64}([0.47736325015237235, 0.33970530481730987]), Ferrite.Node{2, Float64}([0.46406870923827526, 0.45073264476884783]), Ferrite.Node{2, Float64}([0.5673862832010109, 0.6502733177441277]), Ferrite.Node{2, Float64}([0.3659322703052846, 0.48721715390516906]), Ferrite.Node{2, Float64}([0.67739173818206, 0.47779490322723067]), Ferrite.Node{2, Float64}([0.5225603855864607, 0.5344788107726228]), Ferrite.Node{2, Float64}([0.473387537304908, 0.6402608858085923])], Dict{String, OrderedCollections.OrderedSet{Int64}}("4" => OrderedCollections.OrderedSet{Int64}([114, 105, 110, 111, 92, 99, 98, 89, 106, 96 … 102, 93, 113, 97, 88, 104, 91, 87, 100, 101]), "1" => OrderedCollections.OrderedSet{Int64}([5, 16, 20, 12, 24, 28, 8, 17, 30, 1 … 33, 4, 13, 15, 2, 10, 18, 21, 26, 27]), "5" => OrderedCollections.OrderedSet{Int64}([140, 123, 122, 133, 137, 136, 117, 135, 121, 115 … 139, 128, 124, 129, 131, 120, 127, 116, 132, 134]), "2" => OrderedCollections.OrderedSet{Int64}([35, 52, 37, 53, 47, 41, 43, 45, 36, 44 … 39, 51, 46, 40, 48, 34, 50, 54, 38, 42]), "6" => OrderedCollections.OrderedSet{Int64}([158, 169, 173, 141, 148, 160, 154, 171, 145, 146 … 168, 174, 163, 147, 165, 144, 142, 167, 157, 150]), "3" => OrderedCollections.OrderedSet{Int64}([78, 56, 55, 79, 81, 58, 60, 72, 75, 83 … 71, 66, 76, 57, 59, 70, 65, 63, 86, 62])), Dict{String, OrderedCollections.OrderedSet{Int64}}(), Dict{String, OrderedCollections.OrderedSet{Ferrite.FacetIndex}}("left" => OrderedCollections.OrderedSet{Ferrite.FacetIndex}([Ferrite.FacetIndex((96, 1)), Ferrite.FacetIndex((101, 1)), Ferrite.FacetIndex((108, 1)), Ferrite.FacetIndex((115, 2)), Ferrite.FacetIndex((118, 2)), Ferrite.FacetIndex((122, 1)), Ferrite.FacetIndex((125, 1)), Ferrite.FacetIndex((127, 1))]), "bottom" => OrderedCollections.OrderedSet{Ferrite.FacetIndex}([Ferrite.FacetIndex((69, 1)), Ferrite.FacetIndex((73, 1)), Ferrite.FacetIndex((78, 1)), Ferrite.FacetIndex((79, 1)), Ferrite.FacetIndex((90, 3)), Ferrite.FacetIndex((104, 1)), Ferrite.FacetIndex((106, 1)), Ferrite.FacetIndex((107, 1))]), "top" => OrderedCollections.OrderedSet{Ferrite.FacetIndex}([Ferrite.FacetIndex((21, 1)), Ferrite.FacetIndex((23, 1)), Ferrite.FacetIndex((28, 1)), Ferrite.FacetIndex((35, 1)), Ferrite.FacetIndex((36, 2)), Ferrite.FacetIndex((42, 3)), Ferrite.FacetIndex((43, 3)), Ferrite.FacetIndex((47, 3))])), Dict{String, OrderedCollections.OrderedSet{Ferrite.VertexIndex}}()), 762), true))P-multigrid Configuration
1. Galerkin Coarsening Strategy
config_gal = pmultigrid_config(coarse_strategy = Galerkin())
x_gal, res_gal = solve(A, b,fe_space, config_gal;B = B, log=true, rtol = 1e-10)([-0.010066052417735043, 0.027797092959891082, -0.01240874414390502, 0.027640763158586292, -0.011854385123482435, 0.03365190209634847, -0.011198274908231555, 0.02783032460427619, -0.012028337125774731, 0.030805535642147414 … -0.002718167122881313, 0.018354231517249135, -0.006543360153482845, 0.028078201021824636, -0.0025044709085162626, 0.01695653579607672, -0.0032230431473498503, 0.01922965599077098, -0.003825764715161167, 0.020257286724516228], [3980.227721655242, 455.14707317023976, 83.10562617425046, 19.830089514541793, 5.472289795751415, 1.6740301037733865, 0.5591426284601044, 0.20049948415025493, 0.07549461363046683, 0.029266726620580375, 0.011527948122701805, 0.004579028070799461, 0.0018268659786891135, 0.0007305725473577155, 0.0002925414805076431, 0.00011723074201158901, 4.7000082133767856e-5])2. Rediscretization Coarsening Strategy
# Rediscretization Coarsening Strategy
config_red = pmultigrid_config(coarse_strategy = Rediscretization(LinearElasticityMultigrid(C)))
x_red, res_red = solve(A, b, fe_space, config_red; B = B, log=true, rtol = 1e-10)([-0.01006605241774941, 0.027797092959884813, -0.012408744143935217, 0.027640763158565375, -0.011854385123504447, 0.0336519020963316, -0.011198274908250067, 0.027830324604262383, -0.012028337125788498, 0.03080553564213324 … -0.0027181671228574743, 0.01835423151724516, -0.00654336015361581, 0.028078201022162425, -0.0025044709085287322, 0.016956535796113524, -0.0032230431473469004, 0.019229655990823005, -0.003825764715173804, 0.020257286724577155], [3980.227721655242, 457.06819374818394, 83.22542839601425, 19.868450689341543, 5.501109849713366, 1.6926739908520694, 0.5693868243804461, 0.20550627263045834, 0.07776351831307342, 0.03024997140091354, 0.01194347964883253, 0.004752219743190657, 0.0018984979075725849, 0.0007600703281681786, 0.0003046574923322057, 0.0001221995027470287, 4.9035686222845414e-5])Test the solution
using Test
@testset "Linear Elasticity Example" begin
println("Final residual with Galerkin coarsening: ", res_gal[end])
@test A * x_gal ≈ b
println("Final residual with Rediscretization coarsening: ", res_red[end])
@test A * x_red ≈ b
endTest.DefaultTestSet("Linear Elasticity Example", Any[], 2, false, false, true, 1.753141719357153e9, 1.753141719357334e9, false, "linear_elasticity.md")Plain program
Here follows a version of the program without any comments. The file is also available here: linear_elasticity.jl.
using Ferrite, FerriteGmsh, SparseArrays
using Downloads: download
Emod = 200.0e3 # Young's modulus [MPa]
ν = 0.3 # Poisson's ratio [-]
Gmod = Emod / (2(1 + ν)) # Shear modulus
Kmod = Emod / (3(1 - 2ν)) # Bulk modulus
C = gradient(ϵ -> 2 * Gmod * dev(ϵ) + 3 * Kmod * vol(ϵ), zero(SymmetricTensor{2,2}))
function assemble_external_forces!(f_ext, dh, facetset, facetvalues, prescribed_traction)
# Create a temporary array for the facet's local contributions to the external force vector
fe_ext = zeros(getnbasefunctions(facetvalues))
for facet in FacetIterator(dh, facetset)
# Update the facetvalues to the correct facet number
reinit!(facetvalues, facet)
# Reset the temporary array for the next facet
fill!(fe_ext, 0.0)
# Access the cell's coordinates
cell_coordinates = getcoordinates(facet)
for qp in 1:getnquadpoints(facetvalues)
# Calculate the global coordinate of the quadrature point.
x = spatial_coordinate(facetvalues, qp, cell_coordinates)
tₚ = prescribed_traction(x)
# Get the integration weight for the current quadrature point.
dΓ = getdetJdV(facetvalues, qp)
for i in 1:getnbasefunctions(facetvalues)
Nᵢ = shape_value(facetvalues, qp, i)
fe_ext[i] += tₚ ⋅ Nᵢ * dΓ
end
end
# Add the local contributions to the correct indices in the global external force vector
assemble!(f_ext, celldofs(facet), fe_ext)
end
return f_ext
end
function assemble_cell!(ke, cellvalues, C)
for q_point in 1:getnquadpoints(cellvalues)
# Get the integration weight for the quadrature point
dΩ = getdetJdV(cellvalues, q_point)
for i in 1:getnbasefunctions(cellvalues)
# Gradient of the test function
∇Nᵢ = shape_gradient(cellvalues, q_point, i)
for j in 1:getnbasefunctions(cellvalues)
# Symmetric gradient of the trial function
∇ˢʸᵐNⱼ = shape_symmetric_gradient(cellvalues, q_point, j)
ke[i, j] += (∇Nᵢ ⊡ C ⊡ ∇ˢʸᵐNⱼ) * dΩ
end
end
end
return ke
end
function assemble_global!(K, dh, cellvalues, C)
# Allocate the element stiffness matrix
n_basefuncs = getnbasefunctions(cellvalues)
ke = zeros(n_basefuncs, n_basefuncs)
# Create an assembler
assembler = start_assemble(K)
# Loop over all cells
for cell in CellIterator(dh)
# Update the shape function gradients based on the cell coordinates
reinit!(cellvalues, cell)
# Reset the element stiffness matrix
fill!(ke, 0.0)
# Compute element contribution
assemble_cell!(ke, cellvalues, C)
# Assemble ke into K
assemble!(assembler, celldofs(cell), ke)
end
return K
end
function linear_elasticity_2d(C)
logo_mesh = "logo.geo"
asset_url = "https://raw.githubusercontent.com/Ferrite-FEM/Ferrite.jl/gh-pages/assets/"
isfile(logo_mesh) || download(string(asset_url, logo_mesh), logo_mesh)
grid = togrid(logo_mesh)
addfacetset!(grid, "top", x -> x[2] ≈ 1.0) # facets for which x[2] ≈ 1.0 for all nodes
addfacetset!(grid, "left", x -> abs(x[1]) < 1.0e-6)
addfacetset!(grid, "bottom", x -> abs(x[2]) < 1.0e-6)
dim = 2
order = 2 # quadratic interpolation
ip = Lagrange{RefTriangle,order}()^dim # vector valued interpolation
qr = QuadratureRule{RefTriangle}(4) # 4 quadrature point
qr_face = FacetQuadratureRule{RefTriangle}(1)
cellvalues = CellValues(qr, ip)
facetvalues = FacetValues(qr_face, ip)
dh = DofHandler(grid)
add!(dh, :u, ip)
close!(dh)
ch = ConstraintHandler(dh)
add!(ch, Dirichlet(:u, getfacetset(grid, "bottom"), (x, t) -> 0.0, 2))
add!(ch, Dirichlet(:u, getfacetset(grid, "left"), (x, t) -> 0.0, 1))
close!(ch)
traction(x) = Vec(0.0, 20.0e3 * x[1])
A = allocate_matrix(dh)
assemble_global!(A, dh, cellvalues, C)
b = zeros(ndofs(dh))
assemble_external_forces!(b, dh, getfacetset(grid, "top"), facetvalues, traction)
apply!(A, b, ch)
return A, b, dh, cellvalues, ch
end
function create_nns(dh)
##Ndof = ndofs(dh)
grid = dh.grid
Ndof = 2 * (grid.nodes |> length) # nns at p = 1 for AMG
B = zeros(Float64, Ndof, 3)
B[1:2:end, 1] .= 1 # x - translation
B[2:2:end, 2] .= 1 # y - translation
# in-plane rotation (x,y) → (-y,x)
coords = reduce(hcat, grid.nodes .|> (n -> n.x |> collect))' # convert nodes to 2d array
y = coords[:, 2]
x = coords[:, 1]
B[1:2:end, 3] .= -y
B[2:2:end, 3] .= x
return B
end
using FerriteMultigrid
A, b, dh, cellvalues, ch = linear_elasticity_2d(C);
B = create_nns(dh)
fe_space = FESpace(dh, cellvalues, ch)
config_gal = pmultigrid_config(coarse_strategy = Galerkin())
x_gal, res_gal = solve(A, b,fe_space, config_gal;B = B, log=true, rtol = 1e-10)
# Rediscretization Coarsening Strategy
config_red = pmultigrid_config(coarse_strategy = Rediscretization(LinearElasticityMultigrid(C)))
x_red, res_red = solve(A, b, fe_space, config_red; B = B, log=true, rtol = 1e-10)
using Test
@testset "Linear Elasticity Example" begin
println("Final residual with Galerkin coarsening: ", res_gal[end])
@test A * x_gal ≈ b
println("Final residual with Rediscretization coarsening: ", res_red[end])
@test A * x_red ≈ b
endThis page was generated using Literate.jl.