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.000914681s, CPU 0.000909s)
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.00184661s, CPU 0.001846s)
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: 405ms / 74.8% 60.5MiB / 82.0%
Section ncalls time %tot avg alloc %tot avg
---------------------------------------------------------------------------------------
Rediscretization only 1 95.5ms 31.5% 95.5ms 13.8MiB 27.9% 13.8MiB
solve! 1 83.2ms 27.5% 83.2ms 4.73MiB 9.5% 4.73MiB
Postsmoother 100 24.9ms 8.2% 249μs 0.00B 0.0% 0.00B
Presmoother 100 24.8ms 8.2% 248μs 0.00B 0.0% 0.00B
Residual eval 100 8.05ms 2.7% 80.5μs 0.00B 0.0% 0.00B
Coarse solve 100 2.32ms 0.8% 23.2μs 230KiB 0.5% 2.30KiB
Presmoother 200 830μs 0.3% 4.15μs 0.00B 0.0% 0.00B
Postsmoother 200 793μs 0.3% 3.97μs 0.00B 0.0% 0.00B
Residual eval 200 246μs 0.1% 1.23μs 0.00B 0.0% 0.00B
Restriction 200 111μs 0.0% 553ns 0.00B 0.0% 0.00B
Prolongation 200 107μs 0.0% 534ns 0.00B 0.0% 0.00B
Coarse solve 100 15.2μs 0.0% 152ns 0.00B 0.0% 0.00B
Prolongation 100 470μs 0.2% 4.70μs 0.00B 0.0% 0.00B
Restriction 100 415μs 0.1% 4.15μs 0.00B 0.0% 0.00B
init 1 12.3ms 4.1% 12.3ms 9.10MiB 18.3% 9.10MiB
extend_hierarchy! 1 11.7ms 3.9% 11.7ms 8.48MiB 17.1% 8.48MiB
assembly 1 11.1ms 3.7% 11.1ms 4.39MiB 8.8% 4.39MiB
row normalization 1 11.0μs 0.0% 11.0μs 0.00B 0.0% 0.00B
coarse solver setup 1 443μs 0.1% 443μs 546KiB 1.1% 546KiB
extend_hierarchy! 2 414μs 0.1% 207μs 533KiB 1.0% 267KiB
improve candidates 2 106μs 0.0% 53.0μs 0.00B 0.0% 0.00B
fit candidates 2 105μs 0.0% 52.6μs 82.2KiB 0.2% 41.1KiB
RAP 2 99.0μs 0.0% 49.5μs 175KiB 0.3% 87.4KiB
restriction setup 2 55.2μs 0.0% 27.6μs 192KiB 0.4% 96.0KiB
strength 2 28.5μs 0.0% 14.2μs 51.6KiB 0.1% 25.8KiB
aggregation 2 15.4μs 0.0% 7.69μs 24.0KiB 0.0% 12.0KiB
coarse solver setup 1 21.5μs 0.0% 21.5μs 3.19KiB 0.0% 3.19KiB
ml setup 1 3.39μs 0.0% 3.39μs 224B 0.0% 224B
prologue 1 1.19μs 0.0% 1.19μs 6.93KiB 0.0% 6.93KiB
coarsen order 1 143μs 0.0% 143μs 61.4KiB 0.1% 61.4KiB
Galerkin only 1 83.9ms 27.7% 83.9ms 14.9MiB 30.0% 14.9MiB
solve! 1 70.0ms 23.1% 70.0ms 4.74MiB 9.5% 4.74MiB
Presmoother 100 25.0ms 8.2% 250μs 0.00B 0.0% 0.00B
Postsmoother 100 24.8ms 8.2% 248μs 0.00B 0.0% 0.00B
Residual eval 100 8.07ms 2.7% 80.7μs 0.00B 0.0% 0.00B
Coarse solve 100 2.30ms 0.8% 23.0μs 233KiB 0.5% 2.33KiB
Presmoother 200 806μs 0.3% 4.03μs 0.00B 0.0% 0.00B
Postsmoother 200 788μs 0.3% 3.94μs 0.00B 0.0% 0.00B
Residual eval 200 253μs 0.1% 1.27μs 0.00B 0.0% 0.00B
Restriction 200 131μs 0.0% 654ns 0.00B 0.0% 0.00B
Prolongation 200 101μs 0.0% 507ns 0.00B 0.0% 0.00B
Coarse solve 100 13.3μs 0.0% 133ns 0.00B 0.0% 0.00B
Prolongation 100 463μs 0.2% 4.63μs 0.00B 0.0% 0.00B
Restriction 100 418μs 0.1% 4.18μs 0.00B 0.0% 0.00B
init 1 13.9ms 4.6% 13.9ms 10.2MiB 20.5% 10.2MiB
extend_hierarchy! 1 13.1ms 4.3% 13.1ms 9.08MiB 18.3% 9.08MiB
build prolongator 1 12.0ms 4.0% 12.0ms 4.42MiB 8.9% 4.42MiB
assembly 1 12.0ms 4.0% 12.0ms 4.39MiB 8.8% 4.39MiB
row normalization 1 11.2μs 0.0% 11.2μs 0.00B 0.0% 0.00B
RAP 1 1.13ms 0.4% 1.13ms 4.65MiB 9.4% 4.65MiB
build restriction 1 190ns 0.0% 190ns 48.0B 0.0% 48.0B
coarse solver setup 1 563μs 0.2% 563μs 1.01MiB 2.0% 1.01MiB
extend_hierarchy! 2 513μs 0.2% 256μs 1.00MiB 2.0% 511KiB
restriction setup 2 164μs 0.1% 81.9μs 676KiB 1.3% 338KiB
fit candidates 2 114μs 0.0% 57.0μs 85.7KiB 0.2% 42.9KiB
improve candidates 2 107μs 0.0% 53.6μs 0.00B 0.0% 0.00B
RAP 2 88.7μs 0.0% 44.3μs 191KiB 0.4% 95.5KiB
strength 2 26.5μs 0.0% 13.2μs 52.1KiB 0.1% 26.1KiB
aggregation 2 6.17μs 0.0% 3.09μs 8.92KiB 0.0% 4.46KiB
coarse solver setup 1 35.8μs 0.0% 35.8μs 3.19KiB 0.0% 3.19KiB
ml setup 1 5.30μs 0.0% 5.30μs 224B 0.0% 224B
prologue 1 3.76μs 0.0% 3.76μs 6.93KiB 0.0% 6.93KiB
coarsen order 1 166μs 0.1% 166μs 61.4KiB 0.1% 61.4KiB
CG 1 59.9ms 19.8% 59.9ms 92.8KiB 0.2% 92.8KiB
Build preconditioner 2 25.5ms 8.4% 12.7ms 19.3MiB 38.8% 9.64MiB
pmultigrid hierarchy 2 25.4ms 8.4% 12.7ms 19.3MiB 38.8% 9.64MiB
extend_hierarchy! 2 24.1ms 8.0% 12.1ms 17.6MiB 35.4% 8.78MiB
assembly 1 11.2ms 3.7% 11.2ms 4.39MiB 8.8% 4.39MiB
build prolongator 1 11.2ms 3.7% 11.2ms 4.42MiB 8.9% 4.42MiB
assembly 1 11.2ms 3.7% 11.2ms 4.39MiB 8.8% 4.39MiB
row normalization 1 11.1μs 0.0% 11.1μs 0.00B 0.0% 0.00B
RAP 1 1.10ms 0.4% 1.10ms 4.65MiB 9.4% 4.65MiB
row normalization 1 11.9μs 0.0% 11.9μs 0.00B 0.0% 0.00B
build restriction 1 120ns 0.0% 120ns 48.0B 0.0% 48.0B
coarse solver setup 2 935μs 0.3% 467μs 1.54MiB 3.1% 789KiB
extend_hierarchy! 4 853μs 0.3% 213μs 1.52MiB 3.1% 388KiB
improve candidates 4 225μs 0.1% 56.2μs 0.00B 0.0% 0.00B
fit candidates 4 216μs 0.1% 54.1μs 168KiB 0.3% 42.0KiB
RAP 4 172μs 0.1% 43.1μs 366KiB 0.7% 91.4KiB
restriction setup 4 157μs 0.1% 39.2μs 868KiB 1.7% 217KiB
strength 4 55.0μs 0.0% 13.8μs 104KiB 0.2% 25.9KiB
aggregation 4 19.0μs 0.0% 4.74μs 32.9KiB 0.1% 8.22KiB
coarse solver setup 2 63.0μs 0.0% 31.5μs 6.38KiB 0.0% 3.19KiB
ml setup 2 8.92μs 0.0% 4.46μs 448B 0.0% 224B
prologue 2 3.36μs 0.0% 1.68μs 13.9KiB 0.0% 6.93KiB
coarsen order 2 303μs 0.1% 151μs 123KiB 0.2% 61.4KiB
Rediscretization CG 1 19.7ms 6.5% 19.7ms 804KiB 1.6% 804KiB
Postsmoother 28 6.99ms 2.3% 250μs 0.00B 0.0% 0.00B
Presmoother 28 6.99ms 2.3% 250μs 0.00B 0.0% 0.00B
Residual eval 28 2.25ms 0.7% 80.3μs 0.00B 0.0% 0.00B
Coarse solve 28 647μs 0.2% 23.1μs 66.0KiB 0.1% 2.36KiB
Presmoother 56 229μs 0.1% 4.10μs 0.00B 0.0% 0.00B
Postsmoother 56 229μs 0.1% 4.08μs 0.00B 0.0% 0.00B
Residual eval 56 68.8μs 0.0% 1.23μs 0.00B 0.0% 0.00B
Restriction 56 30.9μs 0.0% 552ns 0.00B 0.0% 0.00B
Prolongation 56 26.9μs 0.0% 480ns 0.00B 0.0% 0.00B
Coarse solve 28 4.25μs 0.0% 152ns 0.00B 0.0% 0.00B
Prolongation 28 129μs 0.0% 4.61μs 0.00B 0.0% 0.00B
Restriction 28 115μs 0.0% 4.11μs 0.00B 0.0% 0.00B
Galerkin CG 1 18.4ms 6.1% 18.4ms 754KiB 1.5% 754KiB
Postsmoother 26 6.55ms 2.2% 252μs 0.00B 0.0% 0.00B
Presmoother 26 6.46ms 2.1% 248μs 0.00B 0.0% 0.00B
Residual eval 26 2.13ms 0.7% 81.8μs 0.00B 0.0% 0.00B
Coarse solve 26 610μs 0.2% 23.5μs 62.3KiB 0.1% 2.39KiB
Presmoother 52 217μs 0.1% 4.18μs 0.00B 0.0% 0.00B
Postsmoother 52 208μs 0.1% 4.00μs 0.00B 0.0% 0.00B
Residual eval 52 64.4μs 0.0% 1.24μs 0.00B 0.0% 0.00B
Restriction 52 30.5μs 0.0% 586ns 0.00B 0.0% 0.00B
Prolongation 52 27.0μs 0.0% 519ns 0.00B 0.0% 0.00B
Coarse solve 26 4.46μs 0.0% 171ns 0.00B 0.0% 0.00B
Prolongation 26 121μs 0.0% 4.65μs 0.00B 0.0% 0.00B
Restriction 26 107μs 0.0% 4.12μ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.772040140780767e9, 1.772040140781106e9, false, "linear_elasticity.md", Random.Xoshiro(0x5016a5094b6f4d5a, 0x9a0b1b7f9cb75763, 0xab7a9652a60f98a1, 0x23aadadbc358c2a2, 0x7148e327b750a116))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.