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, FerriteOperators, 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)
dhh = DofHandlerHierarchy(grid, 2)
add!(dhh, :u, [ip_coarse, ip])
close!(dhh)
chh = ConstraintHandlerHierarchy(dhh)
add!(chh, dh->Dirichlet(:u, getfacetset(dh.grid, "bottom"), (x, t) -> 0.0, 2))
add!(chh, dh->Dirichlet(:u, getfacetset(dh.grid, "left"), (x, t) -> 0.0, 1))
close!(chh)
traction(x) = Vec(0.0, 20.0e3 * x[1])
dh = dhh[end]
ch = chh[end]
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, dhh, chh
endlinear_elasticity_2d (generic function with 1 method)Rediscretization
For the rediscretization approach we can define the assembly for the levels via FerriteOperators directly.
"""
LinearElasticityIntegrator{TC, QRC} <: AbstractBilinearIntegrator
Multigrid problem for linear elasticity.
Implements the FerriteOperators `AbstractBilinearIntegrator` interface.
"""
struct LinearElasticityIntegrator{TC <: SymmetricTensor, QRC} <: FerriteOperators.AbstractBilinearIntegrator
ℂ::TC # material stiffness tensor (4th order)
qrc::QRC
end
struct LinearElasticityElementCache{CV, TC} <: FerriteOperators.AbstractVolumetricElementCache
cv::CV
ℂ::TC
end
function FerriteOperators.setup_element_cache(problem::LinearElasticityIntegrator, sdh::SubDofHandler)
qr = getquadraturerule(problem.qrc, sdh)
ip = Ferrite.getfieldinterpolation(sdh, first(Ferrite.getfieldnames(sdh)))
first_cell = getcells(Ferrite.get_grid(sdh.dh), first(sdh.cellset))
ip_geo = Ferrite.geometric_interpolation(typeof(first_cell))
cv = CellValues(qr, ip, ip_geo)
return LinearElasticityElementCache(cv, problem.ℂ)
end
function FerriteOperators.assemble_element!(Ke::AbstractMatrix, cell::CellCache, cache::LinearElasticityElementCache, p)
reinit!(cache.cv, cell)
fill!(Ke, 0.0)
ℂ = cache.ℂ
for q_point in 1:getnquadpoints(cache.cv)
dΩ = getdetJdV(cache.cv, q_point)
for i in 1:getnbasefunctions(cache.cv)
∇ˢʸᵐNᵢ = shape_symmetric_gradient(cache.cv, q_point, i)
for j in 1:getnbasefunctions(cache.cv)
∇ˢʸᵐNⱼ = shape_symmetric_gradient(cache.cv, q_point, j)
Ke[i, j] += (∇ˢʸᵐNᵢ ⊡ ℂ ⊡ ∇ˢʸᵐNⱼ) * dΩ
end
end
end
return Ke
endNear 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 4th order polynomial shape functions.
A, b, dhh, chh = 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.00102546s, CPU 0.001026s)
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.00191031s, CPU 0.001909s)
Info : 104 nodes 245 elementsConstruct the near null space (NNS) matrix
B = create_nns(dhh[1])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.
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 = solve(A, b, dhh, chh, config_gal; pcoarse_solver, log=true, maxiter = 1000, rtol = 1e-10)
builder_gal = PMultigridPreconBuilder(dhh, chh, config_gal; pcoarse_solver)
@timeit "Build preconditioner" Pl_gal = builder_gal(A)[1]
@timeit "Galerkin CG" x_gcg, res_gcg = IterativeSolvers.cg(A, b; Pl = Pl_gal, maxiter = 1000, log=true, verbose=false)([-0.00992277315765481, 0.027908992912983734, -0.012273544893880071, 0.027787662912528526, -0.011716098095084209, 0.03380248712387573, -0.010480498299943412, 0.027964097067350556, -0.011058650996367192, 0.02796370463977638 … -0.003571423563060232, 0.020447373994073623, -0.0034252014243340086, 0.01966749967042404, -0.0064993177094255315, 0.02796650150961365, -0.007055602260338829, 0.029169542071724543, -0.0070905775655036255, 0.029186469208947057], Converged after 25 iterations.)2. Rediscretization Coarsening Strategy
# Rediscretization Coarsening Strategy
config_red = pmultigrid_config(coarse_strategy = Rediscretization(LinearElasticityIntegrator(C, QuadratureRuleCollection(7))))
@timeit "Rediscretization only" x_red, res_red = solve(A, b, dhh, chh, config_red; pcoarse_solver, log=true, maxiter = 1000, rtol = 1e-10)
builder_red = PMultigridPreconBuilder(dhh, chh, config_red; pcoarse_solver)
@timeit "Build preconditioner" Pl_red = builder_red(A)[1]
@timeit "Rediscretization CG" x_rcg, res_rcg = IterativeSolvers.cg(A, b; Pl = Pl_red, maxiter = 1000, log=true, verbose=false)
print_timer(title = "Analysis with $(getncells(dhh[end].grid)) elements", linechars = :ascii)-----------------------------------------------------------------------------------------------
Analysis with 174 elements Time Allocations
----------------------- ------------------------
Tot / % measured: 1.48s / 88.0% 119MiB / 86.8%
Section ncalls time %tot avg alloc %tot avg
-----------------------------------------------------------------------------------------------
Rediscretization only 1 961ms 73.8% 961ms 61.2MiB 59.1% 61.2MiB
init 1 281ms 21.6% 281ms 29.1MiB 28.2% 29.1MiB
pmultigrid numeric 1 279ms 21.4% 279ms 20.6MiB 20.0% 20.6MiB
setup coarse operator 1 178ms 13.7% 178ms 15.5MiB 15.0% 15.5MiB
assemble coarse operator 1 100ms 7.7% 100ms 4.55MiB 4.4% 4.55MiB
coarse solver setup 1 610μs 0.0% 610μs 557KiB 0.5% 557KiB
extend_hierarchy! 2 493μs 0.0% 246μs 539KiB 0.5% 270KiB
fit candidates 2 144μs 0.0% 72.0μs 82.1KiB 0.1% 41.0KiB
improve candidates 2 116μs 0.0% 57.8μs 0.00B 0.0% 0.00B
RAP 2 96.6μs 0.0% 48.3μs 180KiB 0.2% 90.0KiB
restriction setup 2 71.7μs 0.0% 35.9μs 194KiB 0.2% 97.0KiB
strength 2 34.1μs 0.0% 17.1μs 51.9KiB 0.0% 26.0KiB
aggregation 2 23.4μs 0.0% 11.7μs 21.9KiB 0.0% 10.9KiB
smoother setup 2 551ns 0.0% 276ns 704B 0.0% 352B
coarse solver setup 1 101μs 0.0% 101μs 7.64KiB 0.0% 7.64KiB
prologue 1 5.53μs 0.0% 5.53μs 6.93KiB 0.0% 6.93KiB
smoother setup 1 601ns 0.0% 601ns 576B 0.0% 576B
pmultigrid symbolic 1 1.40ms 0.1% 1.40ms 8.48MiB 8.2% 8.48MiB
build prolongator 1 1.39ms 0.1% 1.39ms 8.48MiB 8.2% 8.48MiB
setup transfer operator 1 948μs 0.1% 948μs 8.45MiB 8.2% 8.45MiB
assemble transfer operator 1 371μs 0.0% 371μs 1.44KiB 0.0% 1.44KiB
row normalization 1 37.0μs 0.0% 37.0μs 0.00B 0.0% 0.00B
build restriction 1 90.0ns 0.0% 90.0ns 48.0B 0.0% 48.0B
solve! 1 76.0ms 5.8% 76.0ms 4.66MiB 4.5% 4.66MiB
Postsmoother 97 26.0ms 2.0% 268μs 0.00B 0.0% 0.00B
Presmoother 97 26.0ms 2.0% 268μs 0.00B 0.0% 0.00B
Residual eval 97 8.73ms 0.7% 90.0μs 0.00B 0.0% 0.00B
Coarse solve 97 2.56ms 0.2% 26.4μs 293KiB 0.3% 3.02KiB
Presmoother 194 888μs 0.1% 4.58μs 0.00B 0.0% 0.00B
Postsmoother 194 839μs 0.1% 4.32μs 0.00B 0.0% 0.00B
Residual eval 194 266μs 0.0% 1.37μs 0.00B 0.0% 0.00B
Restriction 194 129μs 0.0% 663ns 0.00B 0.0% 0.00B
Prolongation 194 110μs 0.0% 569ns 0.00B 0.0% 0.00B
Coarse solve 97 92.6μs 0.0% 954ns 69.7KiB 0.1% 736B
Restriction 97 1.58ms 0.1% 16.3μs 0.00B 0.0% 0.00B
Prolongation 97 1.39ms 0.1% 14.3μs 0.00B 0.0% 0.00B
Build preconditioner 2 139ms 10.7% 69.7ms 15.0MiB 14.5% 7.51MiB
pmultigrid hierarchy 2 134ms 10.3% 66.9ms 5.59MiB 5.4% 2.80MiB
RAP numeric 1 3.68ms 0.3% 3.68ms 752B 0.0% 752B
coarse solver setup 2 1.68ms 0.1% 842μs 1.46MiB 1.4% 748KiB
extend_hierarchy! 4 1.44ms 0.1% 361μs 1.43MiB 1.4% 366KiB
improve candidates 4 429μs 0.0% 107μs 0.00B 0.0% 0.00B
RAP 4 353μs 0.0% 88.3μs 379KiB 0.4% 94.6KiB
fit candidates 4 295μs 0.0% 73.8μs 168KiB 0.2% 41.9KiB
restriction setup 4 210μs 0.0% 52.5μs 621KiB 0.6% 155KiB
strength 4 111μs 0.0% 27.8μs 251KiB 0.2% 62.7KiB
aggregation 4 30.2μs 0.0% 7.54μs 28.6KiB 0.0% 7.16KiB
smoother setup 4 891ns 0.0% 223ns 1.38KiB 0.0% 352B
coarse solver setup 2 219μs 0.0% 110μs 15.3KiB 0.0% 7.64KiB
prologue 2 4.54μs 0.0% 2.27μs 13.9KiB 0.0% 6.93KiB
assemble coarse operator 1 334μs 0.0% 334μs 528B 0.0% 528B
setup coarse operator 1 141μs 0.0% 141μs 107KiB 0.1% 107KiB
smoother setup 2 2.88μs 0.0% 1.44μs 1.12KiB 0.0% 576B
RAP symbolic 1 4.24ms 0.3% 4.24ms 967KiB 0.9% 967KiB
build prolongator 1 1.35ms 0.1% 1.35ms 8.48MiB 8.2% 8.48MiB
setup transfer operator 1 911μs 0.1% 911μs 8.45MiB 8.2% 8.45MiB
assemble transfer operator 1 382μs 0.0% 382μs 1.44KiB 0.0% 1.44KiB
row normalization 1 36.7μs 0.0% 36.7μs 0.00B 0.0% 0.00B
build restriction 1 100ns 0.0% 100ns 48.0B 0.0% 48.0B
Galerkin only 1 91.5ms 7.0% 91.5ms 17.1MiB 16.6% 17.1MiB
solve! 1 80.7ms 6.2% 80.7ms 4.71MiB 4.6% 4.71MiB
Presmoother 98 26.2ms 2.0% 267μs 0.00B 0.0% 0.00B
Postsmoother 98 26.2ms 2.0% 267μs 0.00B 0.0% 0.00B
Residual eval 98 8.81ms 0.7% 89.9μs 0.00B 0.0% 0.00B
Coarse solve 98 6.53ms 0.5% 66.7μs 299KiB 0.3% 3.05KiB
Postsmoother 196 2.54ms 0.2% 13.0μs 0.00B 0.0% 0.00B
Presmoother 196 2.53ms 0.2% 12.9μs 0.00B 0.0% 0.00B
Residual eval 196 828μs 0.1% 4.22μs 0.00B 0.0% 0.00B
Restriction 196 142μs 0.0% 723ns 0.00B 0.0% 0.00B
Prolongation 196 123μs 0.0% 628ns 0.00B 0.0% 0.00B
Coarse solve 98 116μs 0.0% 1.19μs 70.4KiB 0.1% 736B
Restriction 98 1.66ms 0.1% 17.0μs 0.00B 0.0% 0.00B
Prolongation 98 1.40ms 0.1% 14.2μs 0.00B 0.0% 0.00B
init 1 10.8ms 0.8% 10.8ms 12.4MiB 12.0% 12.4MiB
pmultigrid symbolic 1 5.94ms 0.5% 5.94ms 11.5MiB 11.1% 11.5MiB
RAP symbolic 1 4.38ms 0.3% 4.38ms 967KiB 0.9% 967KiB
build prolongator 1 1.43ms 0.1% 1.43ms 8.48MiB 8.2% 8.48MiB
setup transfer operator 1 989μs 0.1% 989μs 8.45MiB 8.2% 8.45MiB
assemble transfer operator 1 376μs 0.0% 376μs 1.44KiB 0.0% 1.44KiB
row normalization 1 36.8μs 0.0% 36.8μs 0.00B 0.0% 0.00B
build restriction 1 100ns 0.0% 100ns 48.0B 0.0% 48.0B
pmultigrid numeric 1 4.84ms 0.4% 4.84ms 975KiB 0.9% 975KiB
RAP numeric 1 3.74ms 0.3% 3.74ms 752B 0.0% 752B
coarse solver setup 1 1.08ms 0.1% 1.08ms 944KiB 0.9% 944KiB
extend_hierarchy! 2 945μs 0.1% 472μs 926KiB 0.9% 463KiB
improve candidates 2 315μs 0.0% 158μs 0.00B 0.0% 0.00B
RAP 2 224μs 0.0% 112μs 199KiB 0.2% 99.3KiB
fit candidates 2 153μs 0.0% 76.3μs 85.4KiB 0.1% 42.7KiB
restriction setup 2 139μs 0.0% 69.6μs 427KiB 0.4% 214KiB
strength 2 95.0μs 0.0% 47.5μs 199KiB 0.2% 99.3KiB
aggregation 2 10.5μs 0.0% 5.26μs 6.78KiB 0.0% 3.39KiB
smoother setup 2 521ns 0.0% 260ns 704B 0.0% 352B
coarse solver setup 1 124μs 0.0% 124μs 7.64KiB 0.0% 7.64KiB
prologue 1 3.37μs 0.0% 3.37μs 6.93KiB 0.0% 6.93KiB
smoother setup 1 691ns 0.0% 691ns 576B 0.0% 576B
CG 1 66.5ms 5.1% 66.5ms 92.8KiB 0.1% 92.8KiB
Rediscretization CG 1 22.2ms 1.7% 22.2ms 834KiB 0.8% 834KiB
Presmoother 28 7.57ms 0.6% 270μs 0.00B 0.0% 0.00B
Postsmoother 28 7.46ms 0.6% 266μs 0.00B 0.0% 0.00B
Residual eval 28 2.59ms 0.2% 92.5μs 0.00B 0.0% 0.00B
Coarse solve 28 770μs 0.1% 27.5μs 86.1KiB 0.1% 3.08KiB
Presmoother 56 255μs 0.0% 4.55μs 0.00B 0.0% 0.00B
Postsmoother 56 250μs 0.0% 4.46μs 0.00B 0.0% 0.00B
Residual eval 56 84.6μs 0.0% 1.51μs 0.00B 0.0% 0.00B
Restriction 56 41.6μs 0.0% 744ns 0.00B 0.0% 0.00B
Prolongation 56 32.2μs 0.0% 574ns 0.00B 0.0% 0.00B
Coarse solve 28 31.3μs 0.0% 1.12μs 20.1KiB 0.0% 736B
Restriction 28 457μs 0.0% 16.3μs 0.00B 0.0% 0.00B
Prolongation 28 414μs 0.0% 14.8μs 0.00B 0.0% 0.00B
Galerkin CG 1 20.5ms 1.6% 20.5ms 757KiB 0.7% 757KiB
Presmoother 25 6.60ms 0.5% 264μs 0.00B 0.0% 0.00B
Postsmoother 25 6.58ms 0.5% 263μs 0.00B 0.0% 0.00B
Residual eval 25 2.23ms 0.2% 89.0μs 0.00B 0.0% 0.00B
Coarse solve 25 1.66ms 0.1% 66.5μs 77.9KiB 0.1% 3.12KiB
Presmoother 50 631μs 0.0% 12.6μs 0.00B 0.0% 0.00B
Postsmoother 50 627μs 0.0% 12.5μs 0.00B 0.0% 0.00B
Residual eval 50 235μs 0.0% 4.69μs 0.00B 0.0% 0.00B
Restriction 50 41.4μs 0.0% 828ns 0.00B 0.0% 0.00B
Prolongation 50 32.1μs 0.0% 641ns 0.00B 0.0% 0.00B
Coarse solve 25 25.9μs 0.0% 1.03μs 18.0KiB 0.0% 736B
Restriction 25 417μs 0.0% 16.7μs 0.00B 0.0% 0.00B
Prolongation 25 348μs 0.0% 13.9μs 0.00B 0.0% 0.00B
build prolongator 1 1.51ms 0.1% 1.51ms 8.48MiB 8.2% 8.48MiB
setup transfer operator 1 1.02ms 0.1% 1.02ms 8.45MiB 8.2% 8.45MiB
assemble transfer operator 1 427μs 0.0% 427μs 1.44KiB 0.0% 1.44KiB
row normalization 1 36.9μs 0.0% 36.9μs 0.00B 0.0% 0.00B
build restriction 1 90.0ns 0.0% 90.0ns 48.0B 0.0% 48.0B
-----------------------------------------------------------------------------------------------Plain program
Here follows a version of the program without any comments. The file is also available here: linear_elasticity.jl.
using Ferrite, FerriteGmsh, FerriteOperators, 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)
dhh = DofHandlerHierarchy(grid, 2)
add!(dhh, :u, [ip_coarse, ip])
close!(dhh)
chh = ConstraintHandlerHierarchy(dhh)
add!(chh, dh->Dirichlet(:u, getfacetset(dh.grid, "bottom"), (x, t) -> 0.0, 2))
add!(chh, dh->Dirichlet(:u, getfacetset(dh.grid, "left"), (x, t) -> 0.0, 1))
close!(chh)
traction(x) = Vec(0.0, 20.0e3 * x[1])
dh = dhh[end]
ch = chh[end]
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, dhh, chh
end
"""
LinearElasticityIntegrator{TC, QRC} <: AbstractBilinearIntegrator
Multigrid problem for linear elasticity.
Implements the FerriteOperators `AbstractBilinearIntegrator` interface.
"""
struct LinearElasticityIntegrator{TC <: SymmetricTensor, QRC} <: FerriteOperators.AbstractBilinearIntegrator
ℂ::TC # material stiffness tensor (4th order)
qrc::QRC
end
struct LinearElasticityElementCache{CV, TC} <: FerriteOperators.AbstractVolumetricElementCache
cv::CV
ℂ::TC
end
function FerriteOperators.setup_element_cache(problem::LinearElasticityIntegrator, sdh::SubDofHandler)
qr = getquadraturerule(problem.qrc, sdh)
ip = Ferrite.getfieldinterpolation(sdh, first(Ferrite.getfieldnames(sdh)))
first_cell = getcells(Ferrite.get_grid(sdh.dh), first(sdh.cellset))
ip_geo = Ferrite.geometric_interpolation(typeof(first_cell))
cv = CellValues(qr, ip, ip_geo)
return LinearElasticityElementCache(cv, problem.ℂ)
end
function FerriteOperators.assemble_element!(Ke::AbstractMatrix, cell::CellCache, cache::LinearElasticityElementCache, p)
reinit!(cache.cv, cell)
fill!(Ke, 0.0)
ℂ = cache.ℂ
for q_point in 1:getnquadpoints(cache.cv)
dΩ = getdetJdV(cache.cv, q_point)
for i in 1:getnbasefunctions(cache.cv)
∇ˢʸᵐNᵢ = shape_symmetric_gradient(cache.cv, q_point, i)
for j in 1:getnbasefunctions(cache.cv)
∇ˢʸᵐNⱼ = shape_symmetric_gradient(cache.cv, q_point, j)
Ke[i, j] += (∇ˢʸᵐNᵢ ⊡ ℂ ⊡ ∇ˢʸᵐNⱼ) * dΩ
end
end
end
return Ke
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, dhh, chh = linear_elasticity_2d(C);
B = create_nns(dhh[1])
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 = solve(A, b, dhh, chh, config_gal; pcoarse_solver, log=true, maxiter = 1000, rtol = 1e-10)
builder_gal = PMultigridPreconBuilder(dhh, chh, config_gal; pcoarse_solver)
@timeit "Build preconditioner" Pl_gal = builder_gal(A)[1]
@timeit "Galerkin CG" x_gcg, res_gcg = IterativeSolvers.cg(A, b; Pl = Pl_gal, maxiter = 1000, log=true, verbose=false)
# Rediscretization Coarsening Strategy
config_red = pmultigrid_config(coarse_strategy = Rediscretization(LinearElasticityIntegrator(C, QuadratureRuleCollection(7))))
@timeit "Rediscretization only" x_red, res_red = solve(A, b, dhh, chh, config_red; pcoarse_solver, log=true, maxiter = 1000, rtol = 1e-10)
builder_red = PMultigridPreconBuilder(dhh, chh, config_red; pcoarse_solver)
@timeit "Build preconditioner" Pl_red = builder_red(A)[1]
@timeit "Rediscretization CG" x_rcg, res_rcg = IterativeSolvers.cg(A, b; Pl = Pl_red, maxiter = 1000, log=true, verbose=false)
print_timer(title = "Analysis with $(getncells(dhh[end].grid)) elements", linechars = :ascii)This page was generated using Literate.jl.