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:
- Fourth-order
Lagrangeshape functions are used for field approximation:ip = Lagrange{RefTriangle,4}()^2. - High-order quadrature points are used to accommodate the fourth-order shape functions:
qr = QuadratureRule{RefTriangle}(8).
using Ferrite, FerriteGmsh, FerriteMultigrid, AlgebraicMultigrid
using Downloads: download
using IterativeSolvers
using TimerOutputs
TimerOutputs.enable_debug_timings(AlgebraicMultigrid)
TimerOutputs.enable_debug_timings(FerriteMultigrid)
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 = 4
ip = Lagrange{RefTriangle,order}()^dim # vector valued interpolation
ip_coarse = Lagrange{RefTriangle,1}()^dim
qr = QuadratureRule{RefTriangle}(8)
qr_face = FacetQuadratureRule{RefTriangle}(6)
cellvalues = CellValues(qr, ip)
facetvalues = FacetValues(qr_face, ip)
dh = DofHandler(grid)
add!(dh, :u, ip)
close!(dh)
dh_coarse = DofHandler(grid)
add!(dh_coarse, :u, ip_coarse)
close!(dh_coarse)
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, dh_coarse, 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, fieldname = first(dh.field_names))
@assert length(dh.field_names) == 1 "Only a single field is supported for now."
coords_flat = zeros(ndofs(dh))
apply_analytical!(coords_flat, dh, fieldname, x -> x)
coords = reshape(coords_flat, (length(coords_flat) ÷ 2, 2))
grid = dh.grid
B = zeros(Float64, ndofs(dh), 3)
B[1:2:end, 1] .= 1 # x - translation
B[2:2:end, 2] .= 1 # y - translation
# in-plane rotation (x,y) → (-y,x)
x = coords[:, 1]
y = coords[:, 2]
B[1:2:end, 3] .= -y
B[2:2:end, 3] .= x
return B
endcreate_nns (generic function with 2 methods)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, dh_coarse, 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.000933237s, CPU 0.000934s)
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.00191482s, CPU 0.001914s)
Info : 104 nodes 245 elementsConstruct the near null space (NNS) matrix
B = create_nns(dh_coarse)208×3 Matrix{Float64}:
1.0 0.0 -1.0
0.0 1.0 0.801535
1.0 0.0 -0.253172
0.0 1.0 0.454964
1.0 0.0 -0.612999
0.0 1.0 0.900767
1.0 0.0 -0.0858959
0.0 1.0 0.417361
1.0 0.0 -0.66085
0.0 1.0 0.820318
⋮
0.0 1.0 0.11219
1.0 0.0 -0.464069
0.0 1.0 0.592985
1.0 0.0 -0.450733
0.0 1.0 0.301241
1.0 0.0 -0.677392
0.0 1.0 1.0
1.0 0.0 -0.477795
0.0 1.0 0.126586Since 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 = 4}$
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, 4, Ferrite.Lagrange{Ferrite.RefTriangle, 4}}, 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, 4, Ferrite.Lagrange{Ferrite.RefTriangle, 4}}(Ferrite.Lagrange{Ferrite.RefTriangle, 4}())], Int64[], 30)], [:u], [1, 2, 3, 4, 5, 6, 7, 8, 9, 10 … 2519, 2520, 2517, 2518, 2909, 2910, 2911, 2912, 2913, 2914], [1, 31, 61, 91, 121, 151, 181, 211, 241, 271 … 4921, 4951, 4981, 5011, 5041, 5071, 5101, 5131, 5161, 5191], [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.3718154833954542, 0.3875412105674909]), Ferrite.Node{2, Float64}([0.5769806418976113, 0.4322841106613626]), Ferrite.Node{2, Float64}([0.6165746595005668, 0.5544955814857395]), Ferrite.Node{2, Float64}([0.4773632501523721, 0.33970530481730954]), Ferrite.Node{2, Float64}([0.4640687092382758, 0.45073264476884767]), Ferrite.Node{2, Float64}([0.5673862832010109, 0.6502733177441277]), Ferrite.Node{2, Float64}([0.36593227030528475, 0.4872171539051689]), Ferrite.Node{2, Float64}([0.6773917381820608, 0.4777949032272307]), Ferrite.Node{2, Float64}([0.5225603855864607, 0.5344788107726227]), Ferrite.Node{2, Float64}([0.4733875373049082, 0.6402608858085922])], 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}}()), 2914), Ferrite.CellValues{Ferrite.FunctionValues{1, Ferrite.VectorizedInterpolation{2, Ferrite.RefTriangle, 4, Ferrite.Lagrange{Ferrite.RefTriangle, 4}}, 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, 4, Ferrite.Lagrange{Ferrite.RefTriangle, 4}}, 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, 4, Ferrite.Lagrange{Ferrite.RefTriangle, 4}}(Ferrite.Lagrange{Ferrite.RefTriangle, 4}()), Tensors.Vec{2, Float64}[[0.02057613168724234, 0.0] [0.012133933309930203, 0.0] … [0.004244602629577666, 0.0] [-0.007887393245413425, 0.0]; [0.0, 0.02057613168724234] [0.0, 0.012133933309930203] … [0.0, 0.004244602629577666] [0.0, -0.007887393245413425]; … ; [0.39506172839504583, 0.0] [0.46009358509516385, 0.0] … [-0.04976139182375528, 0.0] [0.0985510461872007, 0.0]; [0.0, 0.39506172839504583] [0.0, 0.46009358509516385] … [0.0, -0.04976139182375528] [0.0, 0.0985510461872007]], Tensors.Vec{2, Float64}[[0.02057613168724234, 0.0] [0.012133933309930203, 0.0] … [0.004244602629577666, 0.0] [-0.007887393245413425, 0.0]; [0.0, 0.02057613168724234] [0.0, 0.012133933309930203] … [0.0, 0.004244602629577666] [0.0, -0.007887393245413425]; … ; [0.39506172839504583, 0.0] [0.46009358509516385, 0.0] … [-0.04976139182375528, 0.0] [0.0985510461872007, 0.0]; [0.0, 0.39506172839504583] [0.0, 0.46009358509516385] … [0.0, -0.04976139182375528] [0.0, 0.0985510461872007]], Tensors.Tensor{2, 2, Float64, 4}[[1.4891412422819572 -0.4169756286720052; 0.0 0.0] [-2.800822118983641 0.7842604386350973; 0.0 0.0] … [3.4343339184972983 -0.9616505836213068; 0.0 0.0] [-9.6522014634064 2.7027206412638463; 0.0 0.0]; [0.0 0.0; 1.4891412422819572 -0.4169756286720052] [0.0 0.0; -2.800822118983641 0.7842604386350973] … [0.0 0.0; 3.4343339184972983 -0.9616505836213068] [0.0 0.0; -9.6522014634064 2.7027206412638463]; … ; [-6.2386590191166515 58.569798790751925; 0.0 0.0] [-47.75823294142756 -16.006924476821816; 0.0 0.0] … [6.114845128584763 -69.47423908206534; 0.0 0.0] [124.66677855368695 -35.30736927016434; 0.0 0.0]; [0.0 0.0; -6.2386590191166515 58.569798790751925] [0.0 0.0; -47.75823294142756 -16.006924476821816] … [0.0 0.0; 6.114845128584763 -69.47423908206534] [0.0 0.0; 124.66677855368695 -35.30736927016434]], Tensors.Tensor{2, 2, Float64, 4}[[0.1358024691358129 0.0; 0.0 0.0] [-0.25542141240096017 0.0; 0.0 0.0] … [0.31319462031291806 0.0; 0.0 0.0] [-0.8802340262353081 0.0; 0.0 0.0]; [0.0 0.0; 0.1358024691358129 0.0] [0.0 0.0; -0.25542141240096017 0.0] … [0.0 0.0; 0.31319462031291806 0.0] [0.0 0.0; -0.8802340262353081 0.0]; … ; [3.57245097701588e-14 4.740740740740776; 0.0 0.0] [-4.649482244529567 -2.4511550660646066; 0.0 0.0] … [-0.12081831664643666 -5.653391729417777; 0.0 0.0] [11.365008207313393 -0.03331572693235846; 0.0 0.0]; [0.0 0.0; 3.57245097701588e-14 4.740740740740776] [0.0 0.0; -4.649482244529567 -2.4511550660646066] … [0.0 0.0; -0.12081831664643666 -5.653391729417777] [0.0 0.0; 11.365008207313393 -0.03331572693235846]], nothing, nothing), Ferrite.GeometryMapping{1, Ferrite.Lagrange{Ferrite.RefTriangle, 1}, Matrix{Float64}, Matrix{Tensors.Vec{2, Float64}}, Nothing}(Ferrite.Lagrange{Ferrite.RefTriangle, 1}(), [0.33333333333333 0.45929258829272 … 0.26311282963464 0.00839477740996; 0.33333333333333 0.45929258829272 … 0.00839477740996 0.7284923929554; 0.3333333333333401 0.08141482341455997 … 0.7284923929554 0.26311282963464], 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.072157803838895, 0.04754581713364, 0.04754581713364, 0.04754581713364, 0.05160868526736, 0.05160868526736, 0.05160868526736, 0.0162292488116, 0.0162292488116, 0.0162292488116, 0.013615157087215, 0.013615157087215, 0.013615157087215, 0.013615157087215, 0.013615157087215, 0.013615157087215], Tensors.Vec{2, Float64}[[0.33333333333333, 0.33333333333333], [0.45929258829272, 0.45929258829272], [0.45929258829272, 0.08141482341455], [0.08141482341455, 0.45929258829272], [0.17056930775176, 0.17056930775176], [0.17056930775176, 0.65886138449648], [0.65886138449648, 0.17056930775176], [0.05054722831703, 0.05054722831703], [0.05054722831703, 0.89890554336594], [0.89890554336594, 0.05054722831703], [0.26311282963464, 0.7284923929554], [0.7284923929554, 0.00839477740996], [0.00839477740996, 0.26311282963464], [0.7284923929554, 0.26311282963464], [0.26311282963464, 0.00839477740996], [0.00839477740996, 0.7284923929554]]), [0.0005490069465463841, 0.000361748591238613, 0.000361748591238613, 0.000361748591238613, 0.0003926606022706325, 0.0003926606022706325, 0.0003926606022706325, 0.0001234789566862529, 0.0001234789566862529, 0.0001234789566862529, 0.00010358984644114322, 0.00010358984644114322, 0.00010358984644114322, 0.00010358984644114322, 0.00010358984644114322, 0.00010358984644114322]), Ferrite.ConstraintHandler{Ferrite.DofHandler{2, Ferrite.Grid{2, Ferrite.Triangle, Float64}}, Float64}(Ferrite.Dirichlet[Ferrite.Dirichlet(Main.var"#linear_elasticity_2d##6#linear_elasticity_2d##7"(), 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, 10, 12, 4, 6, 14, 16, 18, 6, 2, 20, 22, 24], [1, 6, 11, 16]), Ferrite.Dirichlet(Main.var"#linear_elasticity_2d##8#linear_elasticity_2d##9"(), 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, 9, 11, 3, 5, 13, 15, 17, 5, 1, 19, 21, 23], [1, 6, 11, 16])], Ferrite.ProjectedDirichlet[], [621, 1250, 1252, 1254, 1256, 1258, 1314, 1318, 1320, 1322 … 2145, 2147, 2149, 2151, 2205, 2207, 2209, 2243, 2245, 2247], [1, 2, 3, 4, 5, 6, 7, 8, 9, 10 … 2905, 2906, 2907, 2908, 2909, 2910, 2911, 2912, 2913, 2914], [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(1900 => 41, 1252 => 3, 2207 => 62, 2013 => 50, 1845 => 32, 2015 => 51, 2061 => 55, 1793 => 28, 1941 => 47, 1852 => 35…), Ferrite.BCValues{Float64}[Ferrite.BCValues{Float64}([1.0 0.0 … 0.5 0.25; 0.0 1.0 … 0.5 0.75; 0.0 0.0 … 0.0 0.0;;; 0.0 0.0 … 0.0 0.0; 1.0 0.0 … 0.5 0.25; 0.0 1.0 … 0.5 0.75;;; 0.0 1.0 … 0.5 0.75; 0.0 0.0 … 0.0 0.0; 1.0 0.0 … 0.5 0.25], [5, 5, 5], 1), Ferrite.BCValues{Float64}([1.0 0.0 … 0.5 0.25; 0.0 1.0 … 0.5 0.75; 0.0 0.0 … 0.0 0.0;;; 0.0 0.0 … 0.0 0.0; 1.0 0.0 … 0.5 0.25; 0.0 1.0 … 0.5 0.75;;; 0.0 1.0 … 0.5 0.75; 0.0 0.0 … 0.0 0.0; 1.0 0.0 … 0.5 0.25], [5, 5, 5], 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, 4, Ferrite.Lagrange{Ferrite.RefTriangle, 4}}(Ferrite.Lagrange{Ferrite.RefTriangle, 4}())], Int64[], 30)], [:u], [1, 2, 3, 4, 5, 6, 7, 8, 9, 10 … 2519, 2520, 2517, 2518, 2909, 2910, 2911, 2912, 2913, 2914], [1, 31, 61, 91, 121, 151, 181, 211, 241, 271 … 4921, 4951, 4981, 5011, 5041, 5071, 5101, 5131, 5161, 5191], [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.3718154833954542, 0.3875412105674909]), Ferrite.Node{2, Float64}([0.5769806418976113, 0.4322841106613626]), Ferrite.Node{2, Float64}([0.6165746595005668, 0.5544955814857395]), Ferrite.Node{2, Float64}([0.4773632501523721, 0.33970530481730954]), Ferrite.Node{2, Float64}([0.4640687092382758, 0.45073264476884767]), Ferrite.Node{2, Float64}([0.5673862832010109, 0.6502733177441277]), Ferrite.Node{2, Float64}([0.36593227030528475, 0.4872171539051689]), Ferrite.Node{2, Float64}([0.6773917381820608, 0.4777949032272307]), Ferrite.Node{2, Float64}([0.5225603855864607, 0.5344788107726227]), Ferrite.Node{2, Float64}([0.4733875373049082, 0.6402608858085922])], 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}}()), 2914), true))P-multigrid Configuration
reset_timer!()
pcoarse_solver = SmoothedAggregationCoarseSolver(; B)SmoothedAggregationCoarseSolver{Tuple{}, Base.Pairs{Symbol, Matrix{Float64}, Nothing, @NamedTuple{B::Matrix{Float64}}}}((), Base.Pairs(:B => [1.0 0.0 -1.0; 0.0 1.0 0.801534880751; … ; 1.0 0.0 -0.4777949032272307; 0.0 1.0 0.1265859930873549]))0. CG as baseline
@timeit "CG" x_cg = IterativeSolvers.cg(A, b; maxiter = 1000, verbose=false)2914-element Vector{Float64}:
-0.009922773142129101
0.027908992919798835
-0.012273544906839417
0.02778766287937867
-0.01171609807787134
0.03380248712544142
-0.010480498291317988
0.027964097065637683
-0.011058650993627232
0.027963704627201345
⋮
0.02044737400079658
-0.003425201397929082
0.01966749967600462
-0.006499317726805468
0.027966501533483217
-0.007055602277993196
0.029169542096681003
-0.007090577588019493
0.0291864692341921651. Galerkin Coarsening Strategy
config_gal = pmultigrid_config(coarse_strategy = Galerkin())
@timeit "Galerkin only" x_gal, res_gal = FerriteMultigrid.solve(A, b,fe_space, config_gal; pcoarse_solver, verbose=false, log=true, rtol = 1e-10)
builder_gal = PMultigridPreconBuilder(fe_space, config_gal; pcoarse_solver)
@timeit "Build preconditioner" Pl_gal = builder_gal(A)[1]
@timeit "Galerkin CG" IterativeSolvers.cg(A, b; Pl = Pl_gal, maxiter = 1000, verbose=false)2914-element Vector{Float64}:
-0.00992277315074982
0.027908992913412128
-0.012273544888830971
0.02778766291311703
-0.011716098087832071
0.03380248711993018
-0.01048049829364379
0.027964097067214727
-0.01105865098999501
0.027963704638106202
⋮
0.02044737399268512
-0.003425201420622597
0.019667499668565236
-0.00649931770467806
0.02796650151137346
-0.007055602255972188
0.029169542075212395
-0.007090577561606642
0.0291864692119266672. Rediscretization Coarsening Strategy
# Rediscretization Coarsening Strategy
config_red = pmultigrid_config(coarse_strategy = Rediscretization(LinearElasticityMultigrid(C)))
@timeit "Rediscretization only" x_red, res_red = solve(A, b, fe_space, config_red; pcoarse_solver, log=true, rtol = 1e-10)
builder_red = PMultigridPreconBuilder(fe_space, config_red; pcoarse_solver)
@timeit "Build preconditioner" Pl_red = builder_red(A)[1]
@timeit "Rediscretization CG" IterativeSolvers.cg(A, b; Pl = Pl_red, maxiter = 1000, verbose=false)
print_timer(title = "Analysis with $(getncells(dh.grid)) elements", linechars = :ascii)---------------------------------------------------------------------------------------
Analysis with 174 elements Time Allocations
----------------------- ------------------------
Tot / % measured: 418ms / 73.1% 60.5MiB / 82.0%
Section ncalls time %tot avg alloc %tot avg
---------------------------------------------------------------------------------------
Rediscretization only 1 98.4ms 32.2% 98.4ms 13.8MiB 27.9% 13.8MiB
solve! 1 85.7ms 28.0% 85.7ms 4.73MiB 9.5% 4.73MiB
Postsmoother 100 24.7ms 8.1% 247μs 0.00B 0.0% 0.00B
Presmoother 100 24.6ms 8.1% 246μs 0.00B 0.0% 0.00B
Residual eval 100 8.07ms 2.6% 80.7μs 0.00B 0.0% 0.00B
Coarse solve 100 2.32ms 0.8% 23.2μs 230KiB 0.5% 2.30KiB
Presmoother 200 809μs 0.3% 4.05μs 0.00B 0.0% 0.00B
Postsmoother 200 800μs 0.3% 4.00μs 0.00B 0.0% 0.00B
Residual eval 200 261μs 0.1% 1.31μs 0.00B 0.0% 0.00B
Restriction 200 113μs 0.0% 565ns 0.00B 0.0% 0.00B
Prolongation 200 96.1μs 0.0% 480ns 0.00B 0.0% 0.00B
Coarse solve 100 18.2μs 0.0% 182ns 0.00B 0.0% 0.00B
Prolongation 100 463μs 0.2% 4.63μs 0.00B 0.0% 0.00B
Restriction 100 429μs 0.1% 4.29μs 0.00B 0.0% 0.00B
init 1 12.6ms 4.1% 12.6ms 9.10MiB 18.3% 9.10MiB
extend_hierarchy! 1 11.9ms 3.9% 11.9ms 8.48MiB 17.1% 8.48MiB
assembly 1 11.3ms 3.7% 11.3ms 4.39MiB 8.8% 4.39MiB
row normalization 1 11.6μs 0.0% 11.6μs 0.00B 0.0% 0.00B
coarse solver setup 1 518μs 0.2% 518μs 546KiB 1.1% 546KiB
extend_hierarchy! 2 469μs 0.2% 235μs 533KiB 1.0% 267KiB
RAP 2 127μs 0.0% 63.4μs 175KiB 0.3% 87.4KiB
fit candidates 2 117μs 0.0% 58.3μs 82.2KiB 0.2% 41.1KiB
improve candidates 2 107μs 0.0% 53.3μs 0.00B 0.0% 0.00B
restriction setup 2 65.2μs 0.0% 32.6μs 192KiB 0.4% 96.0KiB
strength 2 32.2μs 0.0% 16.1μs 51.6KiB 0.1% 25.8KiB
aggregation 2 15.5μs 0.0% 7.77μs 24.0KiB 0.0% 12.0KiB
coarse solver setup 1 38.0μs 0.0% 38.0μs 3.19KiB 0.0% 3.19KiB
ml setup 1 5.48μs 0.0% 5.48μs 224B 0.0% 224B
prologue 1 1.55μs 0.0% 1.55μs 6.93KiB 0.0% 6.93KiB
coarsen order 1 164μs 0.1% 164μs 61.4KiB 0.1% 61.4KiB
Galerkin only 1 82.9ms 27.1% 82.9ms 14.9MiB 30.0% 14.9MiB
solve! 1 69.4ms 22.7% 69.4ms 4.74MiB 9.5% 4.74MiB
Presmoother 100 24.6ms 8.1% 246μs 0.00B 0.0% 0.00B
Postsmoother 100 24.6ms 8.0% 246μs 0.00B 0.0% 0.00B
Residual eval 100 7.95ms 2.6% 79.5μs 0.00B 0.0% 0.00B
Coarse solve 100 2.36ms 0.8% 23.6μs 233KiB 0.5% 2.33KiB
Presmoother 200 817μs 0.3% 4.09μs 0.00B 0.0% 0.00B
Postsmoother 200 791μs 0.3% 3.96μs 0.00B 0.0% 0.00B
Residual eval 200 257μs 0.1% 1.28μs 0.00B 0.0% 0.00B
Restriction 200 135μs 0.0% 675ns 0.00B 0.0% 0.00B
Prolongation 200 112μs 0.0% 559ns 0.00B 0.0% 0.00B
Coarse solve 100 15.9μs 0.0% 159ns 0.00B 0.0% 0.00B
Prolongation 100 459μs 0.2% 4.59μs 0.00B 0.0% 0.00B
Restriction 100 419μs 0.1% 4.19μs 0.00B 0.0% 0.00B
init 1 13.4ms 4.4% 13.4ms 10.2MiB 20.5% 10.2MiB
extend_hierarchy! 1 12.6ms 4.1% 12.6ms 9.08MiB 18.3% 9.08MiB
build prolongator 1 11.5ms 3.8% 11.5ms 4.42MiB 8.9% 4.42MiB
assembly 1 11.5ms 3.7% 11.5ms 4.39MiB 8.8% 4.39MiB
row normalization 1 12.5μs 0.0% 12.5μs 0.00B 0.0% 0.00B
RAP 1 1.14ms 0.4% 1.14ms 4.65MiB 9.4% 4.65MiB
build restriction 1 150ns 0.0% 150ns 48.0B 0.0% 48.0B
coarse solver setup 1 582μs 0.2% 582μs 1.01MiB 2.0% 1.01MiB
extend_hierarchy! 2 518μs 0.2% 259μs 1.00MiB 2.0% 511KiB
restriction setup 2 138μs 0.0% 69.2μs 676KiB 1.3% 338KiB
fit candidates 2 117μs 0.0% 58.6μs 85.7KiB 0.2% 42.9KiB
improve candidates 2 108μs 0.0% 54.0μs 0.00B 0.0% 0.00B
RAP 2 92.1μs 0.0% 46.0μs 191KiB 0.4% 95.5KiB
strength 2 28.1μs 0.0% 14.0μs 52.1KiB 0.1% 26.1KiB
aggregation 2 7.23μs 0.0% 3.62μs 8.92KiB 0.0% 4.46KiB
coarse solver setup 1 39.0μs 0.0% 39.0μs 3.19KiB 0.0% 3.19KiB
prologue 1 13.5μs 0.0% 13.5μs 6.93KiB 0.0% 6.93KiB
ml setup 1 4.96μs 0.0% 4.96μs 224B 0.0% 224B
coarsen order 1 166μs 0.1% 166μs 61.4KiB 0.1% 61.4KiB
CG 1 60.5ms 19.8% 60.5ms 92.8KiB 0.2% 92.8KiB
Build preconditioner 2 26.2ms 8.6% 13.1ms 19.3MiB 38.8% 9.64MiB
pmultigrid hierarchy 2 26.1ms 8.5% 13.1ms 19.3MiB 38.8% 9.64MiB
extend_hierarchy! 2 24.6ms 8.1% 12.3ms 17.6MiB 35.4% 8.78MiB
build prolongator 1 11.4ms 3.7% 11.4ms 4.42MiB 8.9% 4.42MiB
assembly 1 11.4ms 3.7% 11.4ms 4.39MiB 8.8% 4.39MiB
row normalization 1 11.8μs 0.0% 11.8μs 0.00B 0.0% 0.00B
assembly 1 11.4ms 3.7% 11.4ms 4.39MiB 8.8% 4.39MiB
RAP 1 1.11ms 0.4% 1.11ms 4.65MiB 9.4% 4.65MiB
row normalization 1 11.3μs 0.0% 11.3μs 0.00B 0.0% 0.00B
build restriction 1 330ns 0.0% 330ns 48.0B 0.0% 48.0B
coarse solver setup 2 1.06ms 0.3% 530μs 1.54MiB 3.1% 789KiB
extend_hierarchy! 4 954μs 0.3% 239μs 1.52MiB 3.1% 388KiB
fit candidates 4 258μs 0.1% 64.5μs 168KiB 0.3% 42.0KiB
improve candidates 4 215μs 0.1% 53.7μs 0.00B 0.0% 0.00B
restriction setup 4 202μs 0.1% 50.6μs 868KiB 1.7% 217KiB
RAP 4 187μs 0.1% 46.9μs 366KiB 0.7% 91.4KiB
strength 4 60.5μs 0.0% 15.1μs 104KiB 0.2% 25.9KiB
aggregation 4 21.7μs 0.0% 5.43μs 32.9KiB 0.1% 8.22KiB
coarse solver setup 2 73.6μs 0.0% 36.8μs 6.38KiB 0.0% 3.19KiB
prologue 2 14.4μs 0.0% 7.20μs 13.9KiB 0.0% 6.93KiB
ml setup 2 9.77μs 0.0% 4.88μs 448B 0.0% 224B
coarsen order 2 367μs 0.1% 184μs 123KiB 0.2% 61.4KiB
Rediscretization CG 1 19.6ms 6.4% 19.6ms 804KiB 1.6% 804KiB
Presmoother 28 6.93ms 2.3% 248μs 0.00B 0.0% 0.00B
Postsmoother 28 6.89ms 2.3% 246μs 0.00B 0.0% 0.00B
Residual eval 28 2.25ms 0.7% 80.4μs 0.00B 0.0% 0.00B
Coarse solve 28 665μs 0.2% 23.8μs 66.0KiB 0.1% 2.36KiB
Presmoother 56 246μs 0.1% 4.40μs 0.00B 0.0% 0.00B
Postsmoother 56 224μs 0.1% 4.01μs 0.00B 0.0% 0.00B
Residual eval 56 69.3μs 0.0% 1.24μs 0.00B 0.0% 0.00B
Restriction 56 30.9μs 0.0% 551ns 0.00B 0.0% 0.00B
Prolongation 56 27.3μs 0.0% 488ns 0.00B 0.0% 0.00B
Coarse solve 28 5.36μs 0.0% 191ns 0.00B 0.0% 0.00B
Prolongation 28 136μs 0.0% 4.86μs 0.00B 0.0% 0.00B
Restriction 28 116μs 0.0% 4.16μs 0.00B 0.0% 0.00B
Galerkin CG 1 18.1ms 5.9% 18.1ms 754KiB 1.5% 754KiB
Presmoother 26 6.33ms 2.1% 244μs 0.00B 0.0% 0.00B
Postsmoother 26 6.33ms 2.1% 243μs 0.00B 0.0% 0.00B
Residual eval 26 2.14ms 0.7% 82.4μs 0.00B 0.0% 0.00B
Coarse solve 26 617μs 0.2% 23.7μs 62.3KiB 0.1% 2.39KiB
Postsmoother 52 216μs 0.1% 4.16μs 0.00B 0.0% 0.00B
Presmoother 52 216μs 0.1% 4.15μs 0.00B 0.0% 0.00B
Residual eval 52 64.8μs 0.0% 1.25μs 0.00B 0.0% 0.00B
Restriction 52 30.8μs 0.0% 592ns 0.00B 0.0% 0.00B
Prolongation 52 26.9μs 0.0% 517ns 0.00B 0.0% 0.00B
Coarse solve 26 4.37μs 0.0% 168ns 0.00B 0.0% 0.00B
Prolongation 26 132μs 0.0% 5.09μs 0.00B 0.0% 0.00B
Restriction 26 108μs 0.0% 4.15μs 0.00B 0.0% 0.00B
---------------------------------------------------------------------------------------Test the solution
using Test
@testset "Linear Elasticity Example" begin
println("Final residual with Galerkin coarsening: ", res_gal[end])
@test A * x_gal ≈ b atol=1e-4
println("Final residual with Rediscretization coarsening: ", res_red[end])
@test A * x_red ≈ b atol=1e-4
endTest.DefaultTestSet("Linear Elasticity Example", Any[], 2, false, false, true, 1.77506205019592e9, 1.775062050196324e9, false, "linear_elasticity.md", Random.Xoshiro(0xb9f7526810f007ab, 0x6a53ab7670a886fe, 0x9154dfdeabfb935b, 0x5f74a2a3401aedf1, 0x6bcb269b7f8fff03))Plain program
Here follows a version of the program without any comments. The file is also available here: linear_elasticity.jl.
using Ferrite, FerriteGmsh, FerriteMultigrid, AlgebraicMultigrid
using Downloads: download
using IterativeSolvers
using TimerOutputs
TimerOutputs.enable_debug_timings(AlgebraicMultigrid)
TimerOutputs.enable_debug_timings(FerriteMultigrid)
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 = 4
ip = Lagrange{RefTriangle,order}()^dim # vector valued interpolation
ip_coarse = Lagrange{RefTriangle,1}()^dim
qr = QuadratureRule{RefTriangle}(8)
qr_face = FacetQuadratureRule{RefTriangle}(6)
cellvalues = CellValues(qr, ip)
facetvalues = FacetValues(qr_face, ip)
dh = DofHandler(grid)
add!(dh, :u, ip)
close!(dh)
dh_coarse = DofHandler(grid)
add!(dh_coarse, :u, ip_coarse)
close!(dh_coarse)
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, dh_coarse, cellvalues, ch
end
function create_nns(dh, fieldname = first(dh.field_names))
@assert length(dh.field_names) == 1 "Only a single field is supported for now."
coords_flat = zeros(ndofs(dh))
apply_analytical!(coords_flat, dh, fieldname, x -> x)
coords = reshape(coords_flat, (length(coords_flat) ÷ 2, 2))
grid = dh.grid
B = zeros(Float64, ndofs(dh), 3)
B[1:2:end, 1] .= 1 # x - translation
B[2:2:end, 2] .= 1 # y - translation
# in-plane rotation (x,y) → (-y,x)
x = coords[:, 1]
y = coords[:, 2]
B[1:2:end, 3] .= -y
B[2:2:end, 3] .= x
return B
end
using FerriteMultigrid
A, b, dh, dh_coarse, cellvalues, ch = linear_elasticity_2d(C);
B = create_nns(dh_coarse)
fe_space = FESpace(dh, cellvalues, ch)
reset_timer!()
pcoarse_solver = SmoothedAggregationCoarseSolver(; B)
@timeit "CG" x_cg = IterativeSolvers.cg(A, b; maxiter = 1000, verbose=false)
config_gal = pmultigrid_config(coarse_strategy = Galerkin())
@timeit "Galerkin only" x_gal, res_gal = FerriteMultigrid.solve(A, b,fe_space, config_gal; pcoarse_solver, verbose=false, log=true, rtol = 1e-10)
builder_gal = PMultigridPreconBuilder(fe_space, config_gal; pcoarse_solver)
@timeit "Build preconditioner" Pl_gal = builder_gal(A)[1]
@timeit "Galerkin CG" IterativeSolvers.cg(A, b; Pl = Pl_gal, maxiter = 1000, verbose=false)
# Rediscretization Coarsening Strategy
config_red = pmultigrid_config(coarse_strategy = Rediscretization(LinearElasticityMultigrid(C)))
@timeit "Rediscretization only" x_red, res_red = solve(A, b, fe_space, config_red; pcoarse_solver, log=true, rtol = 1e-10)
builder_red = PMultigridPreconBuilder(fe_space, config_red; pcoarse_solver)
@timeit "Build preconditioner" Pl_red = builder_red(A)[1]
@timeit "Rediscretization CG" IterativeSolvers.cg(A, b; Pl = Pl_red, maxiter = 1000, verbose=false)
print_timer(title = "Analysis with $(getncells(dh.grid)) elements", linechars = :ascii)
using Test
@testset "Linear Elasticity Example" begin
println("Final residual with Galerkin coarsening: ", res_gal[end])
@test A * x_gal ≈ b atol=1e-4
println("Final residual with Rediscretization coarsening: ", res_red[end])
@test A * x_red ≈ b atol=1e-4
endThis page was generated using Literate.jl.