From 101be7a0cc3ee48ac2f4a9fb83fe7f638383bfdd Mon Sep 17 00:00:00 2001 From: Abdelrahman95 Date: Sat, 3 May 2025 14:01:10 +0200 Subject: [PATCH 01/24] init elasticity test --- .gitignore | 1 + logo.geo | 45 +++++++++++++++++ test/nns_test.jl | 122 +++++++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 168 insertions(+) create mode 100644 logo.geo create mode 100644 test/nns_test.jl diff --git a/.gitignore b/.gitignore index 97e6a6f..ffd6da3 100644 --- a/.gitignore +++ b/.gitignore @@ -2,3 +2,4 @@ *.jl.*.cov *.jl.mem /Manifest.toml +.vscode diff --git a/logo.geo b/logo.geo new file mode 100644 index 0000000..22434c9 --- /dev/null +++ b/logo.geo @@ -0,0 +1,45 @@ +Point (1) = {1.000000000000,1.000000000000,0.000000000000}; +Point (2) = {1.000000000000,0.379757979263,-0.000000000000}; +Point (3) = {0.801534880751,0.454963545532,-0.000000000000}; +Point (4) = {0.656107331955,1.000000000000,-0.000000000000}; +Point (5) = {0.600672553037,0.775245709538,-0.000000000000}; +Point (6) = {0.000000000000,1.000000000000,0.000000000000}; +Point (7) = {0.392825178821,0.672136259831,-0.000000000000}; +Point (8) = {1.000000000000,0.000000000000,0.000000000000}; +Point (9) = {0.547800422194,-0.000000000000,-0.000000000000}; +Point (10) = {0.488710023938,0.224380304618,-0.000000000000}; +Point (11) = {0.000000000000,0.000000000000,0.000000000000}; +Point (12) = {-0.000000000000,0.324566579562,-0.000000000000}; +Point (13) = {0.172066723668,0.367888021869,-0.000000000000}; +Line (1) = {2,1}; +Line (2) = {1,4}; +Line (3) = {2,3}; +Line (4) = {3,5}; +Line (5) = {5,4}; +Line (6) = {4,6}; +Line (7) = {7,6}; +Line (8) = {5,7}; +Line (9) = {8,2}; +Line (10) = {9,8}; +Line (11) = {9,10}; +Line (12) = {10,3}; +Line (13) = {12,11}; +Line (14) = {11,9}; +Line (15) = {12,13}; +Line (16) = {13,10}; +Line (17) = {6,12}; +Line (18) = {7,13}; +Line Loop (1) = {-1,3,4,5,-2}; +Plane Surface (1) = {1}; Physical Surface (1) = {1}; +Line Loop (2) = {7,-6,-5,8}; +Plane Surface (2) = {2}; Physical Surface (2) = {2}; +Line Loop (3) = {-10,11,12,-3,-9}; +Plane Surface (3) = {3}; Physical Surface (3) = {3}; +Line Loop (4) = {-11,-14,-13,15,16}; +Plane Surface (4) = {4}; Physical Surface (4) = {4}; +Line Loop (5) = {-7,18,-15,-17}; +Plane Surface (5) = {5}; Physical Surface (5) = {5}; +Line Loop (6) = {-16,-18,-8,-4,-12}; +Plane Surface (6) = {6}; Physical Surface (6) = {6}; + +ReverseMesh Surface{1,2,3,4,5,6}; diff --git a/test/nns_test.jl b/test/nns_test.jl new file mode 100644 index 0000000..04190aa --- /dev/null +++ b/test/nns_test.jl @@ -0,0 +1,122 @@ +using AlgebraicMultigrid +using Test + +# Example test: https://ferrite-fem.github.io/Ferrite.jl/stable/tutorials/linear_elasticity/ +using Ferrite, FerriteGmsh, SparseArrays + +using Downloads: download +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 = 1 # linear interpolation +ip = Lagrange{RefTriangle, order}()^dim; # vector valued interpolation + +qr = QuadratureRule{RefTriangle}(1) # 1 quadrature point +qr_face = FacetQuadratureRule{RefTriangle}(1); + +cellvalues = CellValues(qr, ip) +facetvalues = FacetValues(qr_face, ip); + +dh = DofHandler(grid) +add!(dh, :u, ip) +close!(dh); + +ch = ConstraintHandler(dh) +add!(ch, Dirichlet(:u, getfacetset(grid, "bottom"), (x, t) -> 0.0, 2)) +add!(ch, Dirichlet(:u, getfacetset(grid, "left"), (x, t) -> 0.0, 1)) +close!(ch); + +traction(x) = Vec(0.0, 20.0e3 * x[1]); + +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 + +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_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 + +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) +x = A \ b; + +x_amg,residuals = solve(A, b, SmoothedAggregationAMG();log=true) + +#ml = smoothed_aggregation(A) From d550e2647b361f3823bf69eb01ce1d711d12542a Mon Sep 17 00:00:00 2001 From: Abdelrahman95 Date: Sat, 3 May 2025 17:12:45 +0200 Subject: [PATCH 02/24] write the test for nns --- .gitignore | 2 ++ test/nns_test.jl | 34 ++++++++++++++++++++++++++++++++++ 2 files changed, 36 insertions(+) diff --git a/.gitignore b/.gitignore index ffd6da3..d904ed1 100644 --- a/.gitignore +++ b/.gitignore @@ -3,3 +3,5 @@ *.jl.mem /Manifest.toml .vscode +*.mtx +*.csv diff --git a/test/nns_test.jl b/test/nns_test.jl index 04190aa..0875866 100644 --- a/test/nns_test.jl +++ b/test/nns_test.jl @@ -108,10 +108,27 @@ function assemble_global!(K, dh, cellvalues, C) return K end +function create_nns(dh) + Ndof = ndofs(dh) + grid = dh.grid + B = zeros(Float64, Ndof, 3) + B[1:2:end,1] .= 1 # x - translation + B[2:2:end,2] .= 1 # y - translation + + # in-plane rotation (x,y) → (-y,x) + coords = reduce(hcat,grid.nodes .|> (n -> n.x |> collect))' # convert nodes to 2d array + y = coords[:,2] + x = coords[:,1] + B[1:2:end,3] .= -y + B[2:2:end,3] .= x + return B +end + A = allocate_matrix(dh) assemble_global!(A, dh, cellvalues, C); b = zeros(ndofs(dh)) +B = create_nns(dh) assemble_external_forces!(b, dh, getfacetset(grid, "top"), facetvalues, traction); apply!(A, b, ch) @@ -120,3 +137,20 @@ x = A \ b; x_amg,residuals = solve(A, b, SmoothedAggregationAMG();log=true) #ml = smoothed_aggregation(A) +@test A * x_amg ≈ b atol=1e-6 + +## show some results +ml = smoothed_aggregation(A) +println("Number of iterations: $(length(residuals))") +using Printf + +# assuming `residuals` is your Vector of floats +for (i, r) in enumerate(residuals) + @printf("residual at iteration %2d: %6.2e\n", i-1, r) +end + +## Write for python +using DelimitedFiles, MatrixMarket +mmwrite("A.mtx", A) +writedlm("b.csv", b, ',') +writedlm("B.csv", B, ',') From 847c9faba9b1d0b6b3123ce148a50867db237eab Mon Sep 17 00:00:00 2001 From: Abdelrahman95 Date: Sun, 4 May 2025 03:21:32 +0200 Subject: [PATCH 03/24] add qr user defined nns --- src/aggregation.jl | 108 +++++++++++++++++++++++++-------------------- src/multilevel.jl | 13 ++++-- test/nns_test.jl | 48 +++++++++++++++++--- 3 files changed, 110 insertions(+), 59 deletions(-) diff --git a/src/aggregation.jl b/src/aggregation.jl index 4258cdc..efee3df 100644 --- a/src/aggregation.jl +++ b/src/aggregation.jl @@ -1,4 +1,4 @@ -function smoothed_aggregation(A::TA, +function smoothed_aggregation(A::TA, B = nothing, ::Type{Val{bs}}=Val{1}; symmetry = HermitianSymmetry(), strength = SymmetricStrength(), @@ -14,8 +14,8 @@ function smoothed_aggregation(A::TA, coarse_solver = Pinv, kwargs...) where {T,V,bs,TA<:SparseMatrixCSC{T,V}} n = size(A, 1) - # B = kron(ones(n, 1), eye(1)) - B = ones(T,n) + B = isnothing(B) ? ones(T,n) : B + @assert size(A, 1) == size(B, 1) #=max_levels, max_coarse, strength = levelize_strength_or_aggregation(max_levels, max_coarse, strength) @@ -91,56 +91,66 @@ construct_R(::HermitianSymmetry, P) = P' function fit_candidates(AggOp, B, tol = 1e-10) A = adjoint(AggOp) - n_fine, n_coarse = size(A) - n_col = n_coarse - - R = zeros(eltype(B), n_coarse) - Qx = zeros(eltype(B), nnz(A)) - # copy!(Qx, B) - for i = 1:size(Qx, 1) - Qx[i] = B[i] - end - # copy!(A.nzval, B) - for i = 1:n_col - for j in nzrange(A,i) - row = A.rowval[j] - A.nzval[j] = B[row] - end - end - k = 1 - for i = 1:n_col - norm_i = norm_col(A, Qx, i) - threshold_i = tol * norm_i - if norm_i > threshold_i - scale = 1 / norm_i - R[i] = norm_i - else - scale = 0 - R[i] = 0 - end - for j in nzrange(A, i) - row = A.rowval[j] - # Qx[row] *= scale - #@show k - # Qx[k] *= scale - # k += 1 - A.nzval[j] *= scale + n_fine, m = size(B) + n_fine2, n_agg = size(A) + @assert n_fine2 == n_fine + + n_coarse = m * n_agg + T = eltype(B) + Qs = spzeros(T, n_fine, n_coarse) + R = zeros(T, n_coarse, m) + + for agg in 1:n_agg + # fine‐node indices in this aggregate + rows = A.rowval[A.colptr[agg] : A.colptr[agg+1]-1] + + # local near‐nullspace block (length(rows)×m) + M = @view B[rows, :] + + # thin QR ⇒ Qj(length(rows)×m), Rj(m×m) + Qj, Rj = qr(M) + + # offset in global Qs/R + offset = (agg - 1) * m + + # scatter dense Qj into sparse Qs + for local_i in 1:length(rows), local_j in 1:m + val = Qj[local_i, local_j] + if abs(val) >= tol + Qs[rows[local_i], offset + local_j] = val + end end + dropzeros!(Qs) + # stack Rj into R + R[offset+1 : offset+m, :] .= Rj end - # SparseMatrixCSC(size(A)..., A.colptr, A.rowval, Qx), R - A, R + return Qs, R end -function norm_col(A, Qx, i) - s = zero(eltype(A)) - for j in nzrange(A, i) - if A.rowval[j] > length(Qx) - val = 1 - else - val = Qx[A.rowval[j]] + +function qr(A, tol = 1e-10) + T = eltype(A) + m, n = size(A) + Q = similar(A) # m×n, will hold the orthonormal vectors + R = zeros(T, n, n) # n×n upper triangular + + for j in 1:n + # start with the j-th column of A + v = @view A[:, j] # creates a copy of the column + + # subtract off components along previously computed q_i + for i in 1:j-1 + q = @view Q[:,i] + r_ij = q ⋅ v + R[i, j] = r_ij > tol ? r_ij : zero(T) + v = v - R[i, j] * q end - # val = A.nzval[A.rowval[j]] - s += val*val + + # normalize to get q_j + R[j, j] = norm(v) + @assert R[j, j] > tol "Matrix is rank-deficient at column $j" + Q[:, j] = v / R[j, j] end - sqrt(s) + + return Q, R end diff --git a/src/multilevel.jl b/src/multilevel.jl index 40b7b15..6bf9f62 100644 --- a/src/multilevel.jl +++ b/src/multilevel.jl @@ -287,7 +287,12 @@ end abstract type AMGAlg end struct RugeStubenAMG <: AMGAlg end -struct SmoothedAggregationAMG <: AMGAlg end +struct SmoothedAggregationAMG <: AMGAlg + B::Union{<:AbstractMatrix,Nothing} + function SmoothedAggregationAMG(B::Union{AbstractMatrix,Nothing}=nothing) + new(B) + end +end function solve(A::AbstractMatrix, b::Vector, s::AMGAlg, args...; kwargs...) solt = init(s, A, b, args...; kwargs...) @@ -296,9 +301,9 @@ end function init(::RugeStubenAMG, A, b, args...; kwargs...) AMGSolver(ruge_stuben(A; kwargs...), b) end -function init(::SmoothedAggregationAMG, A, b; kwargs...) - AMGSolver(smoothed_aggregation(A; kwargs...), b) +function init(sa::SmoothedAggregationAMG, A, b; kwargs...) + AMGSolver(smoothed_aggregation(A,sa.B; kwargs...), b) end function solve!(solt::AMGSolver, args...; kwargs...) _solve(solt.ml, solt.b, args...; kwargs...) -end \ No newline at end of file +end diff --git a/test/nns_test.jl b/test/nns_test.jl index 0875866..4fcaf9f 100644 --- a/test/nns_test.jl +++ b/test/nns_test.jl @@ -1,6 +1,42 @@ using AlgebraicMultigrid using Test +import AlgebraicMultigrid as AMG +using SparseArrays, LinearAlgebra + +## Test QR factorization +# tolerance and common dimensions +n_fine = 9 +n_agg = 3 +rows = collect(1:n_fine) +agg_id = repeat(1:n_agg, inner=n_fine ÷ n_agg) + +# Agg_input is n_agg×n_fine; we then transpose it +Agg = sparse(agg_id, rows, ones(n_fine), n_agg, n_fine) + +for m in (1, 2, 3) + # build B (n_fine×m) for m = 1, 2, 3 + B = m == 1 ? ones(n_fine,1) : + m == 2 ? hcat(ones(n_fine), collect(1:n_fine)) : + hcat(ones(n_fine), collect(1:n_fine), collect(1:n_fine).^2) + + Qs, R = AMG.fit_candidates(Agg, B) + + @test size(Qs) == (n_fine, m * n_agg) + @test size(R) == (m * n_agg, m) + + for agg in 1:n_agg + idx = findall(x->x==agg, agg_id) + M = B[idx, :] + Qj = Array(Qs[idx, (agg-1)*m+1 : agg*m]) + Rj = R[(agg-1)*m+1 : agg*m, :] + + @test isapprox(M, Qj * Rj; atol=1e-8) + @test isapprox(Qj' * Qj, I(m); atol=1e-8) + @test istriu(Rj) + end +end +## Test Convergance of AMG for linear elasticity # Example test: https://ferrite-fem.github.io/Ferrite.jl/stable/tutorials/linear_elasticity/ using Ferrite, FerriteGmsh, SparseArrays @@ -134,10 +170,10 @@ assemble_external_forces!(b, dh, getfacetset(grid, "top"), facetvalues, traction apply!(A, b, ch) x = A \ b; -x_amg,residuals = solve(A, b, SmoothedAggregationAMG();log=true) +x_amg,residuals = solve(A, b, SmoothedAggregationAMG(B);log=true) #ml = smoothed_aggregation(A) -@test A * x_amg ≈ b atol=1e-6 +@test A * x_amg ≈ b atol=1e-4 ## show some results ml = smoothed_aggregation(A) @@ -150,7 +186,7 @@ for (i, r) in enumerate(residuals) end ## Write for python -using DelimitedFiles, MatrixMarket -mmwrite("A.mtx", A) -writedlm("b.csv", b, ',') -writedlm("B.csv", B, ',') +# using DelimitedFiles, MatrixMarket +# mmwrite("A.mtx", A) +# writedlm("b.csv", b, ',') +# writedlm("B.csv", B, ',') From 99e55506903a787559d30f433cd470297d3b61c9 Mon Sep 17 00:00:00 2001 From: Abdelrahman95 Date: Sat, 17 May 2025 02:15:35 +0200 Subject: [PATCH 04/24] some changes --- Project.toml | 2 ++ src/aggregation.jl | 5 +++-- test/nns_test.jl | 13 ++++++++++--- 3 files changed, 15 insertions(+), 5 deletions(-) diff --git a/Project.toml b/Project.toml index 6891f2b..be57c74 100644 --- a/Project.toml +++ b/Project.toml @@ -6,6 +6,7 @@ version = "1.0.0" CommonSolve = "38540f10-b2f7-11e9-35d8-d573e4eb0ff2" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae" +MatrixMarket = "4d4711f2-db25-561a-b6b3-d35e7d4047d3" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" Reexport = "189a3867-3050-52da-a836-e630ba90ab69" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" @@ -13,6 +14,7 @@ SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" [compat] CommonSolve = "0.2" LinearSolve = "2, 3" +MatrixMarket = "0.5.2" Reexport = "1" julia = "1.6" diff --git a/src/aggregation.jl b/src/aggregation.jl index efee3df..6ca92ab 100644 --- a/src/aggregation.jl +++ b/src/aggregation.jl @@ -14,7 +14,7 @@ function smoothed_aggregation(A::TA, B = nothing, coarse_solver = Pinv, kwargs...) where {T,V,bs,TA<:SparseMatrixCSC{T,V}} n = size(A, 1) - B = isnothing(B) ? ones(T,n) : B + B = isnothing(B) ? ones(T,n,1) : B @assert size(A, 1) == size(B, 1) #=max_levels, max_coarse, strength = @@ -91,10 +91,11 @@ construct_R(::HermitianSymmetry, P) = P' function fit_candidates(AggOp, B, tol = 1e-10) A = adjoint(AggOp) + @show AggOp |> size n_fine, m = size(B) n_fine2, n_agg = size(A) @assert n_fine2 == n_fine - + n_coarse = m * n_agg T = eltype(B) Qs = spzeros(T, n_fine, n_coarse) diff --git a/test/nns_test.jl b/test/nns_test.jl index 4fcaf9f..0790461 100644 --- a/test/nns_test.jl +++ b/test/nns_test.jl @@ -170,13 +170,13 @@ assemble_external_forces!(b, dh, getfacetset(grid, "top"), facetvalues, traction apply!(A, b, ch) x = A \ b; -x_amg,residuals = solve(A, b, SmoothedAggregationAMG(B);log=true) +x_amg,residuals = solve(A, b, SmoothedAggregationAMG(B);log=true); #ml = smoothed_aggregation(A) @test A * x_amg ≈ b atol=1e-4 ## show some results -ml = smoothed_aggregation(A) +ml = smoothed_aggregation(A,B) println("Number of iterations: $(length(residuals))") using Printf @@ -189,4 +189,11 @@ end # using DelimitedFiles, MatrixMarket # mmwrite("A.mtx", A) # writedlm("b.csv", b, ',') -# writedlm("B.csv", B, ',') +# writedlm("B_nns.csv", B, ',') +grid +reduce(hcat,grid.nodes .|> (n -> n.x |> collect))' +E2V = reduce(hcat, grid.cells .|> (c -> c.nodes |> collect) )' +V = reduce(hcat, grid.nodes .|> (n -> n.x |> collect))' + +writedlm("E2V.csv", E2V, ',') +writedlm("V.csv", V, ',') From 59d2c3f3984117cc81cdf637f1225e4391053dec Mon Sep 17 00:00:00 2001 From: Abdelrahman95 Date: Tue, 20 May 2025 00:35:16 +0200 Subject: [PATCH 05/24] rm dead code --- src/aggregation.jl | 2 +- test/nns_test.jl | 26 ++------------------------ 2 files changed, 3 insertions(+), 25 deletions(-) diff --git a/src/aggregation.jl b/src/aggregation.jl index 6ca92ab..30bcdf1 100644 --- a/src/aggregation.jl +++ b/src/aggregation.jl @@ -64,7 +64,7 @@ function extend_hierarchy!(levels, strength, aggregate, smooth, else S, _T = strength(adjoint(A), bsr_flag) end - + # Aggregation operator AggOp = aggregate(S) # b = zeros(eltype(A), size(A, 1)) diff --git a/test/nns_test.jl b/test/nns_test.jl index 0790461..e980a4b 100644 --- a/test/nns_test.jl +++ b/test/nns_test.jl @@ -170,30 +170,8 @@ assemble_external_forces!(b, dh, getfacetset(grid, "top"), facetvalues, traction apply!(A, b, ch) x = A \ b; -x_amg,residuals = solve(A, b, SmoothedAggregationAMG(B);log=true); +x_amg,residuals = solve(A, b, SmoothedAggregationAMG(B);log=true,reltol = 1e-10); #ml = smoothed_aggregation(A) -@test A * x_amg ≈ b atol=1e-4 +@test A * x_amg ≈ b -## show some results -ml = smoothed_aggregation(A,B) -println("Number of iterations: $(length(residuals))") -using Printf - -# assuming `residuals` is your Vector of floats -for (i, r) in enumerate(residuals) - @printf("residual at iteration %2d: %6.2e\n", i-1, r) -end - -## Write for python -# using DelimitedFiles, MatrixMarket -# mmwrite("A.mtx", A) -# writedlm("b.csv", b, ',') -# writedlm("B_nns.csv", B, ',') -grid -reduce(hcat,grid.nodes .|> (n -> n.x |> collect))' -E2V = reduce(hcat, grid.cells .|> (c -> c.nodes |> collect) )' -V = reduce(hcat, grid.nodes .|> (n -> n.x |> collect))' - -writedlm("E2V.csv", E2V, ',') -writedlm("V.csv", V, ',') From ed798f6576dca62489d26260572362072428035c Mon Sep 17 00:00:00 2001 From: Abdelrahman95 Date: Tue, 20 May 2025 00:40:39 +0200 Subject: [PATCH 06/24] minor fixes --- Project.toml | 2 -- src/aggregation.jl | 1 - 2 files changed, 3 deletions(-) diff --git a/Project.toml b/Project.toml index be57c74..6891f2b 100644 --- a/Project.toml +++ b/Project.toml @@ -6,7 +6,6 @@ version = "1.0.0" CommonSolve = "38540f10-b2f7-11e9-35d8-d573e4eb0ff2" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae" -MatrixMarket = "4d4711f2-db25-561a-b6b3-d35e7d4047d3" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" Reexport = "189a3867-3050-52da-a836-e630ba90ab69" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" @@ -14,7 +13,6 @@ SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" [compat] CommonSolve = "0.2" LinearSolve = "2, 3" -MatrixMarket = "0.5.2" Reexport = "1" julia = "1.6" diff --git a/src/aggregation.jl b/src/aggregation.jl index 30bcdf1..d7a15bf 100644 --- a/src/aggregation.jl +++ b/src/aggregation.jl @@ -91,7 +91,6 @@ construct_R(::HermitianSymmetry, P) = P' function fit_candidates(AggOp, B, tol = 1e-10) A = adjoint(AggOp) - @show AggOp |> size n_fine, m = size(B) n_fine2, n_agg = size(A) @assert n_fine2 == n_fine From 674a76bfce7c9e81113ae615679e14d40e96b182 Mon Sep 17 00:00:00 2001 From: Abdelrahman95 Date: Wed, 21 May 2025 11:27:27 +0200 Subject: [PATCH 07/24] some changes --- src/aggregation.jl | 4 +-- test/nns_test.jl | 71 ++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 73 insertions(+), 2 deletions(-) diff --git a/src/aggregation.jl b/src/aggregation.jl index d7a15bf..6882408 100644 --- a/src/aggregation.jl +++ b/src/aggregation.jl @@ -1,4 +1,4 @@ -function smoothed_aggregation(A::TA, B = nothing, +function smoothed_aggregation(A::TA, _B = nothing, ::Type{Val{bs}}=Val{1}; symmetry = HermitianSymmetry(), strength = SymmetricStrength(), @@ -14,7 +14,7 @@ function smoothed_aggregation(A::TA, B = nothing, coarse_solver = Pinv, kwargs...) where {T,V,bs,TA<:SparseMatrixCSC{T,V}} n = size(A, 1) - B = isnothing(B) ? ones(T,n,1) : B + B = isnothing(B) ? ones(T,n,1) : copy(_B) @assert size(A, 1) == size(B, 1) #=max_levels, max_coarse, strength = diff --git a/test/nns_test.jl b/test/nns_test.jl index e980a4b..19a6b77 100644 --- a/test/nns_test.jl +++ b/test/nns_test.jl @@ -36,6 +36,51 @@ for m in (1, 2, 3) end end +@testset "QR examples tests" begin + # Example 1: four nodes, two aggregates, constant vector + AggOp = sparse([1,2,3,4], [1,1,2,2], [1,1,1,1], 4, 2) |> adjoint + B1 = ones(4,1) + Q1, R1 = fit_candidates(AggOp, B1) + Q1a = [0.70710678 0.0; + 0.70710678 0.0; + 0.0 0.70710678; + 0.0 0.70710678] + R1a = [1.41421356; + 1.41421356] + @test Q1 ≈ Q1a + @test R1 ≈ R1a + + # Example 2: constant vector + linear function + B2 = Float64.([1 0; + 1 1; + 1 2; + 1 3]) + Q2, R2 = fit_candidates(AggOp, B2) + Q2a = [ 0.70710678 -0.70710678 0.0 0.0; + 0.70710678 0.70710678 0.0 0.0; + 0.0 0.0 0.70710678 -0.70710678; + 0.0 0.0 0.70710678 0.70710678 ] + R2a = [1.41421356 0.70710678; + 0.0 0.70710678; + 1.41421356 3.53553391; + 0.0 0.70710678] + @test Q2 ≈ Q2a + @test R2 ≈ R2a + + # Example 3: aggregation excludes third node + AggOp3 = sparse([1,2,4], [1,1,2], [1,1,1], 4, 2) |> adjoint + B3 = ones(4,1) + Q3, R3 = fit_candidates(AggOp3, B3) + Q3a = [0.70710678 0.0; + 0.70710678 0.0; + 0.0 0.0; + 0.0 1.0] + R3a = [1.41421356; + 1.0] + @test Q3 ≈ Q3a + @test R3 ≈ R3a +end + ## Test Convergance of AMG for linear elasticity # Example test: https://ferrite-fem.github.io/Ferrite.jl/stable/tutorials/linear_elasticity/ using Ferrite, FerriteGmsh, SparseArrays @@ -175,3 +220,29 @@ x_amg,residuals = solve(A, b, SmoothedAggregationAMG(B);log=true,reltol = 1e-10) #ml = smoothed_aggregation(A) @test A * x_amg ≈ b + + +## DEBUG + +ml = smoothed_aggregation(A,B) +aggregate = StandardAggregation() +improve_candidates = GaussSeidel(iter=4) + + + +strength = SymmetricStrength() +S , _= strength(A, false) +AggOp = aggregate(S) +b = zeros(size(A,1)) +improve_candidates(A, B, b) +T, B = fit_candidates(AggOp, B) + + + +using DelimitedFiles + + + +# write with comma delimiter +#writedlm("B_nns.csv", B, ',') +B_py = readdlm("B_nns_1.csv", ',', Float64) From a3e680ca4e174c9afd631fa3bd4a5431e97a2f36 Mon Sep 17 00:00:00 2001 From: Abdelrahman95 Date: Thu, 22 May 2025 00:47:44 +0200 Subject: [PATCH 08/24] init fix --- src/AlgebraicMultigrid.jl | 2 +- src/aggregation.jl | 53 +++---------- test/nns_test.jl | 163 ++++++++++++++------------------------ 3 files changed, 73 insertions(+), 145 deletions(-) diff --git a/src/AlgebraicMultigrid.jl b/src/AlgebraicMultigrid.jl index 1f91621..c47fc4f 100644 --- a/src/AlgebraicMultigrid.jl +++ b/src/AlgebraicMultigrid.jl @@ -9,7 +9,7 @@ using Printf @reexport import CommonSolve: solve, solve!, init using Reexport -using LinearAlgebra: rmul! +using LinearAlgebra: rmul!,qr include("utils.jl") export approximate_spectral_radius diff --git a/src/aggregation.jl b/src/aggregation.jl index 6882408..3fc1e0f 100644 --- a/src/aggregation.jl +++ b/src/aggregation.jl @@ -14,7 +14,7 @@ function smoothed_aggregation(A::TA, _B = nothing, coarse_solver = Pinv, kwargs...) where {T,V,bs,TA<:SparseMatrixCSC{T,V}} n = size(A, 1) - B = isnothing(B) ? ones(T,n,1) : copy(_B) + B = isnothing(_B) ? ones(T,n,1) : copy(_B) @assert size(A, 1) == size(B, 1) #=max_levels, max_coarse, strength = @@ -88,32 +88,28 @@ function extend_hierarchy!(levels, strength, aggregate, smooth, end construct_R(::HermitianSymmetry, P) = P' -function fit_candidates(AggOp, B, tol = 1e-10) - +function fit_candidates(AggOp, B; tol = 1e-10) A = adjoint(AggOp) - n_fine, m = size(B) - n_fine2, n_agg = size(A) + n_fine, m = size(B) + n_fine2, n_agg = size(A) @assert n_fine2 == n_fine - + n_coarse = m * n_agg T = eltype(B) Qs = spzeros(T, n_fine, n_coarse) R = zeros(T, n_coarse, m) for agg in 1:n_agg - # fine‐node indices in this aggregate rows = A.rowval[A.colptr[agg] : A.colptr[agg+1]-1] + M = @view B[rows, :] # size(rows) × m - # local near‐nullspace block (length(rows)×m) - M = @view B[rows, :] + F = qr(M) + Qfull = Matrix(F.Q) + Qj = Qfull[:, 1:m] + Rj = F.R - # thin QR ⇒ Qj(length(rows)×m), Rj(m×m) - Qj, Rj = qr(M) - - # offset in global Qs/R offset = (agg - 1) * m - # scatter dense Qj into sparse Qs for local_i in 1:length(rows), local_j in 1:m val = Qj[local_i, local_j] if abs(val) >= tol @@ -121,36 +117,9 @@ function fit_candidates(AggOp, B, tol = 1e-10) end end dropzeros!(Qs) - # stack Rj into R + R[offset+1 : offset+m, :] .= Rj end return Qs, R end - -function qr(A, tol = 1e-10) - T = eltype(A) - m, n = size(A) - Q = similar(A) # m×n, will hold the orthonormal vectors - R = zeros(T, n, n) # n×n upper triangular - - for j in 1:n - # start with the j-th column of A - v = @view A[:, j] # creates a copy of the column - - # subtract off components along previously computed q_i - for i in 1:j-1 - q = @view Q[:,i] - r_ij = q ⋅ v - R[i, j] = r_ij > tol ? r_ij : zero(T) - v = v - R[i, j] * q - end - - # normalize to get q_j - R[j, j] = norm(v) - @assert R[j, j] > tol "Matrix is rank-deficient at column $j" - Q[:, j] = v / R[j, j] - end - - return Q, R -end diff --git a/test/nns_test.jl b/test/nns_test.jl index 19a6b77..2cf4395 100644 --- a/test/nns_test.jl +++ b/test/nns_test.jl @@ -4,83 +4,69 @@ import AlgebraicMultigrid as AMG using SparseArrays, LinearAlgebra ## Test QR factorization -# tolerance and common dimensions -n_fine = 9 -n_agg = 3 -rows = collect(1:n_fine) -agg_id = repeat(1:n_agg, inner=n_fine ÷ n_agg) - -# Agg_input is n_agg×n_fine; we then transpose it -Agg = sparse(agg_id, rows, ones(n_fine), n_agg, n_fine) - -for m in (1, 2, 3) - # build B (n_fine×m) for m = 1, 2, 3 - B = m == 1 ? ones(n_fine,1) : - m == 2 ? hcat(ones(n_fine), collect(1:n_fine)) : - hcat(ones(n_fine), collect(1:n_fine), collect(1:n_fine).^2) - - Qs, R = AMG.fit_candidates(Agg, B) - - @test size(Qs) == (n_fine, m * n_agg) - @test size(R) == (m * n_agg, m) - - for agg in 1:n_agg - idx = findall(x->x==agg, agg_id) - M = B[idx, :] - Qj = Array(Qs[idx, (agg-1)*m+1 : agg*m]) - Rj = R[(agg-1)*m+1 : agg*m, :] - - @test isapprox(M, Qj * Rj; atol=1e-8) - @test isapprox(Qj' * Qj, I(m); atol=1e-8) - @test istriu(Rj) +@testset "fit_candidates Unit tests" begin + cases = Vector{Tuple{SparseMatrixCSC{Float64,Int}, Matrix{Float64}}}() + + # 1. Aggregates include all dofs, one candidate + push!(cases, ( + sparse([1,2,3,4,5], [1,1,1,2,2], ones(5), 5, 2), + ones(5, 1) + )) + push!(cases, ( + sparse([1,2,3,4,5], [2,2,1,1,1], ones(5), 5, 2), + ones(5, 1) + )) + push!(cases, ( + sparse([1,2,3,4,5,6,7,8,9], + repeat([1,2,3], inner=3), + ones(9), 9, 3), + ones(9, 1) + )) + push!(cases, ( + sparse([1,2,3,4,5,6,7,8,9], [3,2,1,1,2,3,2,1,3], + ones(9), 9, 3), + reshape(Float64.(0:8), 9, 1) + )) + + # 2. Aggregates include all dofs, two candidates + push!(cases, ( + sparse([1,2,3,4], [1,1,2,2], ones(4), 4, 2), + hcat(ones(4), collect(0:3)) + )) + push!(cases, ( + sparse([1,2,3,4,5,6,7,8,9], repeat([1,2,3], inner=3), ones(9), 9, 3), + hcat(ones(9), collect(0:8)) + )) + push!(cases, ( + sparse([1,2,3,4,5,6,7,8,9], [1,1,2,2,3,3,4,4,4], ones(9), 9, 4), + hcat(ones(9), collect(0:8)) + )) + + # 3. Small norms + push!(cases, ( + sparse([1,2,3,4], [1,1,2,2], ones(4), 4, 2), + hcat(ones(4), 1e-20 .* collect(0:3)) + )) + push!(cases, ( + sparse([1,2,3,4], [1,1,2,2], ones(4), 4, 2), + 1e-20 .* hcat(ones(4), collect(0:3)) + )) + + # Run tests + for (AggOp, fine) in cases + # mask dofs not in aggregation + for i in 1:size(AggOp,1) + if nnz(AggOp[i, :]) == 0 + fine[i, :] .= 0.0 + end + end + Q, R = fit_candidates(AggOp |> adjoint, fine) + # fit exactly and via projection + @test fine ≈ Q * R + @test fine ≈ Q * (Q' * fine) end end -@testset "QR examples tests" begin - # Example 1: four nodes, two aggregates, constant vector - AggOp = sparse([1,2,3,4], [1,1,2,2], [1,1,1,1], 4, 2) |> adjoint - B1 = ones(4,1) - Q1, R1 = fit_candidates(AggOp, B1) - Q1a = [0.70710678 0.0; - 0.70710678 0.0; - 0.0 0.70710678; - 0.0 0.70710678] - R1a = [1.41421356; - 1.41421356] - @test Q1 ≈ Q1a - @test R1 ≈ R1a - - # Example 2: constant vector + linear function - B2 = Float64.([1 0; - 1 1; - 1 2; - 1 3]) - Q2, R2 = fit_candidates(AggOp, B2) - Q2a = [ 0.70710678 -0.70710678 0.0 0.0; - 0.70710678 0.70710678 0.0 0.0; - 0.0 0.0 0.70710678 -0.70710678; - 0.0 0.0 0.70710678 0.70710678 ] - R2a = [1.41421356 0.70710678; - 0.0 0.70710678; - 1.41421356 3.53553391; - 0.0 0.70710678] - @test Q2 ≈ Q2a - @test R2 ≈ R2a - - # Example 3: aggregation excludes third node - AggOp3 = sparse([1,2,4], [1,1,2], [1,1,1], 4, 2) |> adjoint - B3 = ones(4,1) - Q3, R3 = fit_candidates(AggOp3, B3) - Q3a = [0.70710678 0.0; - 0.70710678 0.0; - 0.0 0.0; - 0.0 1.0] - R3a = [1.41421356; - 1.0] - @test Q3 ≈ Q3a - @test R3 ≈ R3a -end - ## Test Convergance of AMG for linear elasticity # Example test: https://ferrite-fem.github.io/Ferrite.jl/stable/tutorials/linear_elasticity/ using Ferrite, FerriteGmsh, SparseArrays @@ -217,32 +203,5 @@ x = A \ b; x_amg,residuals = solve(A, b, SmoothedAggregationAMG(B);log=true,reltol = 1e-10); -#ml = smoothed_aggregation(A) +ml = smoothed_aggregation(A) @test A * x_amg ≈ b - - - -## DEBUG - -ml = smoothed_aggregation(A,B) -aggregate = StandardAggregation() -improve_candidates = GaussSeidel(iter=4) - - - -strength = SymmetricStrength() -S , _= strength(A, false) -AggOp = aggregate(S) -b = zeros(size(A,1)) -improve_candidates(A, B, b) -T, B = fit_candidates(AggOp, B) - - - -using DelimitedFiles - - - -# write with comma delimiter -#writedlm("B_nns.csv", B, ',') -B_py = readdlm("B_nns_1.csv", ',', Float64) From f718a6cba1d1e6431885ef970c5c06c335950feb Mon Sep 17 00:00:00 2001 From: Abdelrahman95 Date: Thu, 22 May 2025 02:03:29 +0200 Subject: [PATCH 09/24] fix test --- src/aggregation.jl | 1 - test/nns_test.jl | 84 +++++++++++++++++++++++++++------------------- 2 files changed, 49 insertions(+), 36 deletions(-) diff --git a/src/aggregation.jl b/src/aggregation.jl index 3fc1e0f..ab7915a 100644 --- a/src/aggregation.jl +++ b/src/aggregation.jl @@ -93,7 +93,6 @@ function fit_candidates(AggOp, B; tol = 1e-10) n_fine, m = size(B) n_fine2, n_agg = size(A) @assert n_fine2 == n_fine - n_coarse = m * n_agg T = eltype(B) Qs = spzeros(T, n_fine, n_coarse) diff --git a/test/nns_test.jl b/test/nns_test.jl index 2cf4395..5b144b1 100644 --- a/test/nns_test.jl +++ b/test/nns_test.jl @@ -2,68 +2,70 @@ using AlgebraicMultigrid using Test import AlgebraicMultigrid as AMG using SparseArrays, LinearAlgebra +using Printf + ## Test QR factorization @testset "fit_candidates Unit tests" begin - cases = Vector{Tuple{SparseMatrixCSC{Float64,Int}, Matrix{Float64}}}() + cases = Vector{Tuple{SparseMatrixCSC{Float64,Int},Matrix{Float64}}}() # 1. Aggregates include all dofs, one candidate push!(cases, ( - sparse([1,2,3,4,5], [1,1,1,2,2], ones(5), 5, 2), + sparse([1, 2, 3, 4, 5], [1, 1, 1, 2, 2], ones(5), 5, 2), ones(5, 1) )) push!(cases, ( - sparse([1,2,3,4,5], [2,2,1,1,1], ones(5), 5, 2), + sparse([1, 2, 3, 4, 5], [2, 2, 1, 1, 1], ones(5), 5, 2), ones(5, 1) )) push!(cases, ( - sparse([1,2,3,4,5,6,7,8,9], - repeat([1,2,3], inner=3), - ones(9), 9, 3), + sparse([1, 2, 3, 4, 5, 6, 7, 8, 9], + repeat([1, 2, 3], inner=3), + ones(9), 9, 3), ones(9, 1) )) push!(cases, ( - sparse([1,2,3,4,5,6,7,8,9], [3,2,1,1,2,3,2,1,3], - ones(9), 9, 3), + sparse([1, 2, 3, 4, 5, 6, 7, 8, 9], [3, 2, 1, 1, 2, 3, 2, 1, 3], + ones(9), 9, 3), reshape(Float64.(0:8), 9, 1) )) # 2. Aggregates include all dofs, two candidates push!(cases, ( - sparse([1,2,3,4], [1,1,2,2], ones(4), 4, 2), + sparse([1, 2, 3, 4], [1, 1, 2, 2], ones(4), 4, 2), hcat(ones(4), collect(0:3)) )) push!(cases, ( - sparse([1,2,3,4,5,6,7,8,9], repeat([1,2,3], inner=3), ones(9), 9, 3), + sparse([1, 2, 3, 4, 5, 6, 7, 8, 9], repeat([1, 2, 3], inner=3), ones(9), 9, 3), hcat(ones(9), collect(0:8)) )) push!(cases, ( - sparse([1,2,3,4,5,6,7,8,9], [1,1,2,2,3,3,4,4,4], ones(9), 9, 4), + sparse([1, 2, 3, 4, 5, 6, 7, 8, 9], [1, 1, 2, 2, 3, 3, 4, 4, 4], ones(9), 9, 4), hcat(ones(9), collect(0:8)) )) # 3. Small norms push!(cases, ( - sparse([1,2,3,4], [1,1,2,2], ones(4), 4, 2), + sparse([1, 2, 3, 4], [1, 1, 2, 2], ones(4), 4, 2), hcat(ones(4), 1e-20 .* collect(0:3)) )) push!(cases, ( - sparse([1,2,3,4], [1,1,2,2], ones(4), 4, 2), + sparse([1, 2, 3, 4], [1, 1, 2, 2], ones(4), 4, 2), 1e-20 .* hcat(ones(4), collect(0:3)) )) # Run tests for (AggOp, fine) in cases # mask dofs not in aggregation - for i in 1:size(AggOp,1) + for i in 1:size(AggOp, 1) if nnz(AggOp[i, :]) == 0 fine[i, :] .= 0.0 end end Q, R = fit_candidates(AggOp |> adjoint, fine) # fit exactly and via projection - @test fine ≈ Q * R - @test fine ≈ Q * (Q' * fine) + @test fine ≈ Q * R + @test fine ≈ Q * (Q' * fine) end end @@ -84,7 +86,7 @@ addfacetset!(grid, "bottom", x -> abs(x[2]) < 1.0e-6); dim = 2 order = 1 # linear interpolation -ip = Lagrange{RefTriangle, order}()^dim; # vector valued interpolation +ip = Lagrange{RefTriangle,order}()^dim; # vector valued interpolation qr = QuadratureRule{RefTriangle}(1) # 1 quadrature point qr_face = FacetQuadratureRule{RefTriangle}(1); @@ -136,7 +138,7 @@ Emod = 200.0e3 # Young's modulus [MPa] 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})); +C = gradient(ϵ -> 2 * Gmod * dev(ϵ) + 3 * Kmod * vol(ϵ), zero(SymmetricTensor{2,2})); function assemble_cell!(ke, cellvalues, C) for q_point in 1:getnquadpoints(cellvalues) @@ -179,29 +181,41 @@ function create_nns(dh) Ndof = ndofs(dh) grid = dh.grid B = zeros(Float64, Ndof, 3) - B[1:2:end,1] .= 1 # x - translation - B[2:2:end,2] .= 1 # y - translation + B[1:2:end, 1] .= 1 # x - translation + B[2:2:end, 2] .= 1 # y - translation # in-plane rotation (x,y) → (-y,x) - coords = reduce(hcat,grid.nodes .|> (n -> n.x |> collect))' # convert nodes to 2d array - y = coords[:,2] - x = coords[:,1] - B[1:2:end,3] .= -y - B[2:2:end,3] .= x + coords = reduce(hcat, grid.nodes .|> (n -> n.x |> collect))' # convert nodes to 2d array + y = coords[:, 2] + x = coords[:, 1] + B[1:2:end, 3] .= -y + B[2:2:end, 3] .= x return B end -A = allocate_matrix(dh) -assemble_global!(A, dh, cellvalues, C); -b = zeros(ndofs(dh)) -B = create_nns(dh) -assemble_external_forces!(b, dh, getfacetset(grid, "top"), facetvalues, traction); +@testset "Linear elaticity test" begin + A = allocate_matrix(dh) + assemble_global!(A, dh, cellvalues, C) + + b = zeros(ndofs(dh)) + B = create_nns(dh) + assemble_external_forces!(b, dh, getfacetset(grid, "top"), facetvalues, traction) + + apply!(A, b, ch) + #x = A \ b -apply!(A, b, ch) -x = A \ b; + x_nns, residuals_nns = solve(A, b, SmoothedAggregationAMG(B); log=true, reltol=1e-10) + x_wonns, residuals_wonns = solve(A, b, SmoothedAggregationAMG(); log=true, reltol=1e-10) -x_amg,residuals = solve(A, b, SmoothedAggregationAMG(B);log=true,reltol = 1e-10); + ml = smoothed_aggregation(A, B) + @show ml -ml = smoothed_aggregation(A) -@test A * x_amg ≈ b + # assuming `residuals` is your Vector of floats + for (i, r) in enumerate(residuals_nns) + @printf("residual at iteration %2d: %6.2e\n", i - 1, r) + end + + @test !(A * x_wonns ≈ b) + @test A * x_nns ≈ b +end From ca23f843f2167b89674667146486cd4ea2017de4 Mon Sep 17 00:00:00 2001 From: Abdelrahman95 Date: Tue, 27 May 2025 11:45:53 +0200 Subject: [PATCH 10/24] init beam test --- test/nns_test.jl | 251 ++++++++++++++++++++++++++++++++++++----------- 1 file changed, 193 insertions(+), 58 deletions(-) diff --git a/test/nns_test.jl b/test/nns_test.jl index 5b144b1..1e8d3bb 100644 --- a/test/nns_test.jl +++ b/test/nns_test.jl @@ -1,12 +1,13 @@ using AlgebraicMultigrid -using Test import AlgebraicMultigrid as AMG +using Test using SparseArrays, LinearAlgebra +using Ferrite, FerriteGmsh, SparseArrays using Printf - +using Downloads: download ## Test QR factorization -@testset "fit_candidates Unit tests" begin +@testset "fit_candidates unit test cases" begin cases = Vector{Tuple{SparseMatrixCSC{Float64,Int},Matrix{Float64}}}() # 1. Aggregates include all dofs, one candidate @@ -69,42 +70,7 @@ using Printf end end -## Test Convergance of AMG for linear elasticity -# Example test: https://ferrite-fem.github.io/Ferrite.jl/stable/tutorials/linear_elasticity/ -using Ferrite, FerriteGmsh, SparseArrays - -using Downloads: download -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 = 1 # linear interpolation -ip = Lagrange{RefTriangle,order}()^dim; # vector valued interpolation - -qr = QuadratureRule{RefTriangle}(1) # 1 quadrature point -qr_face = FacetQuadratureRule{RefTriangle}(1); - -cellvalues = CellValues(qr, ip) -facetvalues = FacetValues(qr_face, ip); - -dh = DofHandler(grid) -add!(dh, :u, ip) -close!(dh); - -ch = ConstraintHandler(dh) -add!(ch, Dirichlet(:u, getfacetset(grid, "bottom"), (x, t) -> 0.0, 2)) -add!(ch, Dirichlet(:u, getfacetset(grid, "left"), (x, t) -> 0.0, 1)) -close!(ch); - -traction(x) = Vec(0.0, 20.0e3 * x[1]); - +## Test Convergance of AMG for linear elasticity & bending beam 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)) @@ -132,14 +98,6 @@ function assemble_external_forces!(f_ext, dh, facetset, facetvalues, prescribed_ return f_ext end -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_cell!(ke, cellvalues, C) for q_point in 1:getnquadpoints(cellvalues) # Get the integration weight for the quadrature point @@ -193,8 +151,44 @@ function create_nns(dh) return B end +function linear_elasticity_2d() + # Example test: https://ferrite-fem.github.io/Ferrite.jl/stable/tutorials/linear_elasticity/ + 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 = 1 # linear interpolation + ip = Lagrange{RefTriangle,order}()^dim # vector valued interpolation + + qr = QuadratureRule{RefTriangle}(1) # 1 quadrature point + qr_face = FacetQuadratureRule{RefTriangle}(1) + + cellvalues = CellValues(qr, ip) + facetvalues = FacetValues(qr_face, ip) + + dh = DofHandler(grid) + add!(dh, :u, ip) + close!(dh) + + ch = ConstraintHandler(dh) + add!(ch, Dirichlet(:u, getfacetset(grid, "bottom"), (x, t) -> 0.0, 2)) + add!(ch, Dirichlet(:u, getfacetset(grid, "left"), (x, t) -> 0.0, 1)) + close!(ch) + + traction(x) = Vec(0.0, 20.0e3 * x[1]) + Emod = 200.0e3 # Young's modulus [MPa] + ν = 0.3 # Poisson's ratio [-] -@testset "Linear elaticity test" begin + 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})) A = allocate_matrix(dh) assemble_global!(A, dh, cellvalues, C) @@ -203,19 +197,160 @@ end assemble_external_forces!(b, dh, getfacetset(grid, "top"), facetvalues, traction) apply!(A, b, ch) - #x = A \ b - x_nns, residuals_nns = solve(A, b, SmoothedAggregationAMG(B); log=true, reltol=1e-10) - x_wonns, residuals_wonns = solve(A, b, SmoothedAggregationAMG(); log=true, reltol=1e-10) + return A, b, B # return the assembled matrix, force vector, and NNS matrix +end + +@testset "Mechanics test cases" begin + @testset "Linear elasticity 2D" begin + A, b, B = linear_elasticity_2d() + + x_nns, residuals_nns = solve(A, b, SmoothedAggregationAMG(B); log=true, reltol=1e-10) + x_wonns, residuals_wonns = solve(A, b, SmoothedAggregationAMG(); log=true, reltol=1e-10) + + ml = smoothed_aggregation(A, B) + @show ml + + # assuming `residuals` is your Vector of floats + for (i, r) in enumerate(residuals_nns) + @printf("residual at iteration %2d: %6.2e\n", i - 1, r) + end + + #test QR factorization on linear elasticity + aggregate = StandardAggregation() + AggOp = aggregate(A) + Q, R = fit_candidates(AggOp, B) + # fit exactly and via projection + @test B ≈ Q * R + @test B ≈ Q * (Q' * B) - ml = smoothed_aggregation(A, B) - @show ml + # Check convergence + @test !(A * x_wonns ≈ b) + @test A * x_nns ≈ b - # assuming `residuals` is your Vector of floats - for (i, r) in enumerate(residuals_nns) - @printf("residual at iteration %2d: %6.2e\n", i - 1, r) end +end + + + +@testset "fit_candidates on cantilever frame beam" begin + # Beam parameters + n_elem = 10 + E = 210e9 # Young's modulus + A = 1e-4 # Cross-section area (for axial) + I = 1e-6 # Moment of inertia (for bending) + L = 1.0 # Total length + le = L / n_elem + n_nodes = n_elem + 1 + dofs_per_node = 3 # u, w, theta + n_dofs = n_nodes * dofs_per_node + + # Element stiffness matrix + function frame_element_stiffness(EA, EI, le) + l2 = le^2 + l3 = le^3 + Ke = zeros(6, 6) + + # Axial (u-u) + Ke[1, 1] = EA / le + Ke[1, 4] = -EA / le + Ke[4, 1] = -EA / le + Ke[4, 4] = EA / le + + # Bending (w, theta) + Kb = EI / l3 * [ + 12.0 6*le -12.0 6*le; + 6*le 4*l2 -6*le 2*l2; + -12.0 -6*le 12.0 -6*le; + 6*le 2*l2 -6*le 4*l2 + ] + idx = [2, 3, 5, 6] + for i in 1:4, j in 1:4 + Ke[idx[i], idx[j]] = Kb[i, j] + end + return Ke + end + + Ke = frame_element_stiffness(E * A, E * I, le) + + # Assemble global stiffness matrix + row = Int[] + col = Int[] + val = Float64[] + for e in 1:n_elem + n1 = e + n2 = e + 1 + dofmap = [ + 3 * n1 - 2, # u₁ + 3 * n1 - 1, # w₁ + 3 * n1, # θ₁ + 3 * n2 - 2, # u₂ + 3 * n2 - 1, # w₂ + 3 * n2 # θ₂ + ] + for i in 1:6, j in 1:6 + push!(row, dofmap[i]) + push!(col, dofmap[j]) + push!(val, Ke[i, j]) + end + end + A = sparse(row, col, val, n_dofs, n_dofs) + # rhs + b = zeros(n_dofs) + force_dof = 3 * (n_nodes - 1) + 2 # w at last node + b[force_dof] = -1000 # Apply downward force at the last node + # Boundary conditions: clamp left end + fixed_dofs = [1, 2, 3] # u₁, w₁, θ₁ + free_dofs = setdiff(1:n_dofs, fixed_dofs) + A_free = A[free_dofs, free_dofs] + b_free = b[free_dofs] + + # NNS matrix (3 columns) + # i = 1 (u) : 1 0 -y + # i = 2 (v) : 0 1 x + # i = 3 (θ) : 0 0 1 + function create_nns_frame(x_coords::Vector{Float64}, dofmap::Vector{Int}) + N = length(dofmap) + B = zeros(N, 3) # 3 rigid body modes: x-translation, y-translation, rotation + + for (i, dof) in enumerate(dofmap) + node = div(dof - 1, 3) + 1 + offset = mod(dof - 1, 3) + x = x_coords[node] + y = 0.0 # 1D beam along x + + if offset == 0 # u (x-translation) + B[i, 1] = 1.0 # translation in x + B[i, 3] = -y # rotation effect on u (−y) + elseif offset == 1 # v (y-translation) + B[i, 2] = 1.0 # translation in y + B[i, 3] = x # rotation effect on v (+x) + elseif offset == 2 # θ (rotation DOF) + B[i, 3] = 1.0 # direct contribution to rigid rotation + end + end + + return B + end + + + # x-coordinates of nodes + x_coords = [le * (i - 1) for i in 1:n_nodes] + B = create_nns_frame(x_coords, free_dofs) + + # Aggregation + aggregate = StandardAggregation() + AggOp = aggregate(A_free) + + # Fit candidates + Q, R = fit_candidates(AggOp, B) + + # Test projection and reconstruction + @test B ≈ Q * R + @test B ≈ Q * (Q' * B) - @test !(A * x_wonns ≈ b) - @test A * x_nns ≈ b + # test solution + # Analaytical solution for cantilever beam + u = A_free \ b_free + @test u[end - 1] ≈ (-1000 * L^3)/ (3*E * I) # vertical disp. at the end of the beam end From b559e12ae414d9b6be144155c9c6dfa045998930 Mon Sep 17 00:00:00 2001 From: Abdelrahman95 Date: Wed, 28 May 2025 08:55:48 +0200 Subject: [PATCH 11/24] add cant. beam test --- src/aggregation.jl | 70 +++++++-------- test/nns_test.jl | 211 ++++++++++++++++++++++++--------------------- 2 files changed, 148 insertions(+), 133 deletions(-) diff --git a/src/aggregation.jl b/src/aggregation.jl index ab7915a..4fe91de 100644 --- a/src/aggregation.jl +++ b/src/aggregation.jl @@ -1,20 +1,20 @@ -function smoothed_aggregation(A::TA, _B = nothing, - ::Type{Val{bs}}=Val{1}; - symmetry = HermitianSymmetry(), - strength = SymmetricStrength(), - aggregate = StandardAggregation(), - smooth = JacobiProlongation(4.0/3.0), - presmoother = GaussSeidel(), - postsmoother = GaussSeidel(), - improve_candidates = GaussSeidel(iter=4), - max_levels = 10, - max_coarse = 10, - diagonal_dominance = false, - keep = false, - coarse_solver = Pinv, kwargs...) where {T,V,bs,TA<:SparseMatrixCSC{T,V}} +function smoothed_aggregation(A::TA, _B=nothing, + ::Type{Val{bs}}=Val{1}; + symmetry=HermitianSymmetry(), + strength=SymmetricStrength(), + aggregate=StandardAggregation(), + smooth=JacobiProlongation(4.0 / 3.0), + presmoother=GaussSeidel(), + postsmoother=GaussSeidel(), + improve_candidates=GaussSeidel(iter=4), + max_levels=10, + max_coarse=10, + diagonal_dominance=false, + keep=false, + coarse_solver=Pinv, kwargs...) where {T,V,bs,TA<:SparseMatrixCSC{T,V}} n = size(A, 1) - B = isnothing(_B) ? ones(T,n,1) : copy(_B) + B = isnothing(_B) ? ones(T, n, 1) : copy(_B) @assert size(A, 1) == size(B, 1) #=max_levels, max_coarse, strength = @@ -28,15 +28,15 @@ function smoothed_aggregation(A::TA, _B = nothing, # agg = [aggregate for _ in 1:max_levels - 1] # sm = [smooth for _ in 1:max_levels] - levels = Vector{Level{TA, TA, Adjoint{T, TA}}}() + levels = Vector{Level{TA,TA,Adjoint{T,TA}}}() bsr_flag = false w = MultiLevelWorkspace(Val{bs}, eltype(A)) residual!(w, size(A, 1)) while length(levels) + 1 < max_levels && size(A, 1) > max_coarse A, B, bsr_flag = extend_hierarchy!(levels, strength, aggregate, smooth, - improve_candidates, diagonal_dominance, - keep, A, B, symmetry, bsr_flag) + improve_candidates, diagonal_dominance, + keep, A, B, symmetry, bsr_flag) coarse_x!(w, size(A, 1)) coarse_b!(w, size(A, 1)) #=if size(A, 1) <= max_coarse @@ -54,9 +54,9 @@ struct HermitianSymmetry end function extend_hierarchy!(levels, strength, aggregate, smooth, - improve_candidates, diagonal_dominance, keep, - A, B, - symmetry, bsr_flag) + improve_candidates, diagonal_dominance, keep, + A, B, + symmetry, bsr_flag) # Calculate strength of connection matrix if symmetry isa HermitianSymmetry @@ -64,13 +64,13 @@ function extend_hierarchy!(levels, strength, aggregate, smooth, else S, _T = strength(adjoint(A), bsr_flag) end - + # Aggregation operator AggOp = aggregate(S) # b = zeros(eltype(A), size(A, 1)) # Improve candidates - b = zeros(size(A,1)) + b = zeros(size(A, 1)) improve_candidates(A, B, b) T, B = fit_candidates(AggOp, B) @@ -88,36 +88,38 @@ function extend_hierarchy!(levels, strength, aggregate, smooth, end construct_R(::HermitianSymmetry, P) = P' -function fit_candidates(AggOp, B; tol = 1e-10) +function fit_candidates(AggOp, B; tol=1e-10) A = adjoint(AggOp) - n_fine, m = size(B) + n_fine, m = size(B) n_fine2, n_agg = size(A) @assert n_fine2 == n_fine n_coarse = m * n_agg T = eltype(B) Qs = spzeros(T, n_fine, n_coarse) - R = zeros(T, n_coarse, m) + R = zeros(T, n_coarse, m) for agg in 1:n_agg - rows = A.rowval[A.colptr[agg] : A.colptr[agg+1]-1] + rows = A.rowval[A.colptr[agg]:A.colptr[agg+1]-1] M = @view B[rows, :] # size(rows) × m - F = qr(M) - Qfull = Matrix(F.Q) - Qj = Qfull[:, 1:m] - Rj = F.R + + F = qr(M) + r = min(length(rows), m) + Qfull = Matrix(F.Q) + Qj = Qfull[:, 1:r] + Rj = F.R offset = (agg - 1) * m - for local_i in 1:length(rows), local_j in 1:m + for local_i in 1:length(rows), local_j in 1:r val = Qj[local_i, local_j] if abs(val) >= tol - Qs[rows[local_i], offset + local_j] = val + Qs[rows[local_i], offset+local_j] = val end end dropzeros!(Qs) - R[offset+1 : offset+m, :] .= Rj + R[offset+1:offset+r, :] .= Rj[1:r, :] end return Qs, R diff --git a/test/nns_test.jl b/test/nns_test.jl index 1e8d3bb..9435d1f 100644 --- a/test/nns_test.jl +++ b/test/nns_test.jl @@ -201,76 +201,70 @@ function linear_elasticity_2d() return A, b, B # return the assembled matrix, force vector, and NNS matrix end -@testset "Mechanics test cases" begin - @testset "Linear elasticity 2D" begin - A, b, B = linear_elasticity_2d() - x_nns, residuals_nns = solve(A, b, SmoothedAggregationAMG(B); log=true, reltol=1e-10) - x_wonns, residuals_wonns = solve(A, b, SmoothedAggregationAMG(); log=true, reltol=1e-10) +## Helper functions for cantilever frame beam ## + +# Element stiffness matrix +function frame_element_stiffness(EA, EI, le) + l2 = le^2 + l3 = le^3 + Ke = zeros(6, 6) + + # Axial (u-u) + Ke[1, 1] = EA / le + Ke[1, 4] = -EA / le + Ke[4, 1] = -EA / le + Ke[4, 4] = EA / le + + # Bending (w, theta) + Kb = EI / l3 * [ + 12.0 6*le -12.0 6*le; + 6*le 4*l2 -6*le 2*l2; + -12.0 -6*le 12.0 -6*le; + 6*le 2*l2 -6*le 4*l2 + ] + idx = [2, 3, 5, 6] + for i in 1:4, j in 1:4 + Ke[idx[i], idx[j]] = Kb[i, j] + end + return Ke +end - ml = smoothed_aggregation(A, B) - @show ml - # assuming `residuals` is your Vector of floats - for (i, r) in enumerate(residuals_nns) - @printf("residual at iteration %2d: %6.2e\n", i - 1, r) +function create_nns_frame(x_coords::Vector{Float64}, dofmap::Vector{Int}) + # NNS matrix (3 columns) + # i = 1 (u) : 1 0 -y + # i = 2 (v) : 0 1 x + # i = 3 (θ) : 0 0 1 + N = length(dofmap) + B = zeros(N, 3) # 3 rigid body modes: x-translation, y-translation, rotation + + for (i, dof) in enumerate(dofmap) + node = div(dof - 1, 3) + 1 + offset = mod(dof - 1, 3) + x = x_coords[node] + y = 0.0 # 1D beam along x + + if offset == 0 # u (x-translation) + B[i, 1] = 1.0 # translation in x + B[i, 3] = -y # rotation effect on u (−y) + elseif offset == 1 # v (y-translation) + B[i, 2] = 1.0 # translation in y + B[i, 3] = x # rotation effect on v (+x) + elseif offset == 2 # θ (rotation DOF) + B[i, 3] = 1.0 # direct contribution to rigid rotation end - - #test QR factorization on linear elasticity - aggregate = StandardAggregation() - AggOp = aggregate(A) - Q, R = fit_candidates(AggOp, B) - # fit exactly and via projection - @test B ≈ Q * R - @test B ≈ Q * (Q' * B) - - # Check convergence - @test !(A * x_wonns ≈ b) - @test A * x_nns ≈ b - end -end - + return B +end -@testset "fit_candidates on cantilever frame beam" begin - # Beam parameters - n_elem = 10 - E = 210e9 # Young's modulus - A = 1e-4 # Cross-section area (for axial) - I = 1e-6 # Moment of inertia (for bending) - L = 1.0 # Total length +function cantilever_beam(P, E,A, I, L, n_elem) le = L / n_elem n_nodes = n_elem + 1 dofs_per_node = 3 # u, w, theta n_dofs = n_nodes * dofs_per_node - # Element stiffness matrix - function frame_element_stiffness(EA, EI, le) - l2 = le^2 - l3 = le^3 - Ke = zeros(6, 6) - - # Axial (u-u) - Ke[1, 1] = EA / le - Ke[1, 4] = -EA / le - Ke[4, 1] = -EA / le - Ke[4, 4] = EA / le - - # Bending (w, theta) - Kb = EI / l3 * [ - 12.0 6*le -12.0 6*le; - 6*le 4*l2 -6*le 2*l2; - -12.0 -6*le 12.0 -6*le; - 6*le 2*l2 -6*le 4*l2 - ] - idx = [2, 3, 5, 6] - for i in 1:4, j in 1:4 - Ke[idx[i], idx[j]] = Kb[i, j] - end - return Ke - end - Ke = frame_element_stiffness(E * A, E * I, le) # Assemble global stiffness matrix @@ -298,59 +292,78 @@ end # rhs b = zeros(n_dofs) force_dof = 3 * (n_nodes - 1) + 2 # w at last node - b[force_dof] = -1000 # Apply downward force at the last node + b[force_dof] = P # Apply downward force at the last node # Boundary conditions: clamp left end fixed_dofs = [1, 2, 3] # u₁, w₁, θ₁ free_dofs = setdiff(1:n_dofs, fixed_dofs) A_free = A[free_dofs, free_dofs] b_free = b[free_dofs] - # NNS matrix (3 columns) - # i = 1 (u) : 1 0 -y - # i = 2 (v) : 0 1 x - # i = 3 (θ) : 0 0 1 - function create_nns_frame(x_coords::Vector{Float64}, dofmap::Vector{Int}) - N = length(dofmap) - B = zeros(N, 3) # 3 rigid body modes: x-translation, y-translation, rotation - - for (i, dof) in enumerate(dofmap) - node = div(dof - 1, 3) + 1 - offset = mod(dof - 1, 3) - x = x_coords[node] - y = 0.0 # 1D beam along x - - if offset == 0 # u (x-translation) - B[i, 1] = 1.0 # translation in x - B[i, 3] = -y # rotation effect on u (−y) - elseif offset == 1 # v (y-translation) - B[i, 2] = 1.0 # translation in y - B[i, 3] = x # rotation effect on v (+x) - elseif offset == 2 # θ (rotation DOF) - B[i, 3] = 1.0 # direct contribution to rigid rotation - end - end + # x-coordinates of nodes + x_coords = [le * (i - 1) for i in 1:n_nodes] + B = create_nns_frame(x_coords, free_dofs) + + return A_free, b_free, B +end + +@testset "Mechanics test cases" begin + @testset "Linear elasticity 2D" begin + A, b, B = linear_elasticity_2d() + + x_nns, residuals_nns = solve(A, b, SmoothedAggregationAMG(B); log=true, reltol=1e-10) + x_wonns, residuals_wonns = solve(A, b, SmoothedAggregationAMG(); log=true, reltol=1e-10) + + ml = smoothed_aggregation(A, B) + @show ml + + @printf("No NNS: final residual at iteration %2d: %6.2e\n", length(residuals_wonns), residuals_nns[end]) + @printf("With NNS: final residual at iteration %2d: %6.2e\n", length(residuals_nns), residuals_wonns[end]) + + #test QR factorization on linear elasticity + aggregate = StandardAggregation() + AggOp = aggregate(A) + Q, R = fit_candidates(AggOp, B) + # fit exactly and via projection + @test B ≈ Q * R + @test B ≈ Q * (Q' * B) + + # Check convergence + @test !(A * x_wonns ≈ b) + @test A * x_nns ≈ b - return B end + @testset "fit_candidates on cantilever frame beam" begin + # Beam parameters + P = -1000.0 # Applied force at the end of the beam + n_elem = 10 + E = 210e9 # Young's modulus + A = 1e-4 # Cross-section area (for axial) + I = 1e-6 # Moment of inertia (for bending) + L = 1.0 # Total length + A, b, B = cantilever_beam(P, E,A, I, L, n_elem) + # test solution + # Analaytical solution for cantilever beam + u = A \ b + @test u[end-1] ≈ (P * L^3) / (3 * E * I) # vertical disp. at the end of the beam - # x-coordinates of nodes - x_coords = [le * (i - 1) for i in 1:n_nodes] - B = create_nns_frame(x_coords, free_dofs) - # Aggregation - aggregate = StandardAggregation() - AggOp = aggregate(A_free) + x_nns, residuals_nns = solve(A, b, SmoothedAggregationAMG(B); log=true, reltol=1e-10) + x_wonns, residuals_wonns = solve(A, b, SmoothedAggregationAMG(); log=true, reltol=1e-10) - # Fit candidates - Q, R = fit_candidates(AggOp, B) + @printf("No NNS: final residual at iteration %2d: %6.2e\n", length(residuals_wonns), residuals_nns[end]) + @printf("With NNS: final residual at iteration %2d: %6.2e\n", length(residuals_nns), residuals_wonns[end]) - # Test projection and reconstruction - @test B ≈ Q * R - @test B ≈ Q * (Q' * B) + # test QR factorization on bending beam + # Aggregation + aggregate = StandardAggregation() + AggOp = aggregate(A) + Q, R = fit_candidates(AggOp, B) + @test B ≈ Q * R + @test B ≈ Q * (Q' * B) - # test solution - # Analaytical solution for cantilever beam - u = A_free \ b_free - @test u[end - 1] ≈ (-1000 * L^3)/ (3*E * I) # vertical disp. at the end of the beam + # Check convergence + @test !(A * x_wonns ≈ b) + @test A * x_nns ≈ b + end end From 517bb95fcd40811172dd1705f2688855c7b040cc Mon Sep 17 00:00:00 2001 From: Abdelrahman95 Date: Wed, 28 May 2025 09:33:24 +0200 Subject: [PATCH 12/24] minor fix --- .gitignore | 1 + logo.geo | 45 --------------------------------------------- 2 files changed, 1 insertion(+), 45 deletions(-) delete mode 100644 logo.geo diff --git a/.gitignore b/.gitignore index d904ed1..93c87aa 100644 --- a/.gitignore +++ b/.gitignore @@ -5,3 +5,4 @@ .vscode *.mtx *.csv +*.geo diff --git a/logo.geo b/logo.geo deleted file mode 100644 index 22434c9..0000000 --- a/logo.geo +++ /dev/null @@ -1,45 +0,0 @@ -Point (1) = {1.000000000000,1.000000000000,0.000000000000}; -Point (2) = {1.000000000000,0.379757979263,-0.000000000000}; -Point (3) = {0.801534880751,0.454963545532,-0.000000000000}; -Point (4) = {0.656107331955,1.000000000000,-0.000000000000}; -Point (5) = {0.600672553037,0.775245709538,-0.000000000000}; -Point (6) = {0.000000000000,1.000000000000,0.000000000000}; -Point (7) = {0.392825178821,0.672136259831,-0.000000000000}; -Point (8) = {1.000000000000,0.000000000000,0.000000000000}; -Point (9) = {0.547800422194,-0.000000000000,-0.000000000000}; -Point (10) = {0.488710023938,0.224380304618,-0.000000000000}; -Point (11) = {0.000000000000,0.000000000000,0.000000000000}; -Point (12) = {-0.000000000000,0.324566579562,-0.000000000000}; -Point (13) = {0.172066723668,0.367888021869,-0.000000000000}; -Line (1) = {2,1}; -Line (2) = {1,4}; -Line (3) = {2,3}; -Line (4) = {3,5}; -Line (5) = {5,4}; -Line (6) = {4,6}; -Line (7) = {7,6}; -Line (8) = {5,7}; -Line (9) = {8,2}; -Line (10) = {9,8}; -Line (11) = {9,10}; -Line (12) = {10,3}; -Line (13) = {12,11}; -Line (14) = {11,9}; -Line (15) = {12,13}; -Line (16) = {13,10}; -Line (17) = {6,12}; -Line (18) = {7,13}; -Line Loop (1) = {-1,3,4,5,-2}; -Plane Surface (1) = {1}; Physical Surface (1) = {1}; -Line Loop (2) = {7,-6,-5,8}; -Plane Surface (2) = {2}; Physical Surface (2) = {2}; -Line Loop (3) = {-10,11,12,-3,-9}; -Plane Surface (3) = {3}; Physical Surface (3) = {3}; -Line Loop (4) = {-11,-14,-13,15,16}; -Plane Surface (4) = {4}; Physical Surface (4) = {4}; -Line Loop (5) = {-7,18,-15,-17}; -Plane Surface (5) = {5}; Physical Surface (5) = {5}; -Line Loop (6) = {-16,-18,-8,-4,-12}; -Plane Surface (6) = {6}; Physical Surface (6) = {6}; - -ReverseMesh Surface{1,2,3,4,5,6}; From d61951194e21ab02c3e88766c773c65ed9340763 Mon Sep 17 00:00:00 2001 From: Abdelrahman95 Date: Wed, 28 May 2025 10:36:20 +0200 Subject: [PATCH 13/24] fix tests + add nns tests to ci --- Project.toml | 5 ++++- src/aggregation.jl | 2 +- test/nns_test.jl | 15 ++++++++------- test/runtests.jl | 1 + 4 files changed, 14 insertions(+), 9 deletions(-) diff --git a/Project.toml b/Project.toml index 6891f2b..f690ceb 100644 --- a/Project.toml +++ b/Project.toml @@ -24,6 +24,9 @@ JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819" LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" +Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6" +Ferrite = "c061ca5d-56c9-439f-9c0e-210fe06d3992" +FerriteGmsh = "4f95f4f8-b27c-4ae5-9a39-ea55e634e36b" [targets] -test = ["DelimitedFiles", "FileIO", "IterativeSolvers", "JLD2", "LinearSolve", "Random", "Test"] +test = ["DelimitedFiles", "FileIO", "IterativeSolvers", "JLD2", "LinearSolve", "Random", "Test","Ferrite","FerriteGmsh","Downloads"] diff --git a/src/aggregation.jl b/src/aggregation.jl index 4fe91de..640a3a0 100644 --- a/src/aggregation.jl +++ b/src/aggregation.jl @@ -90,7 +90,7 @@ construct_R(::HermitianSymmetry, P) = P' function fit_candidates(AggOp, B; tol=1e-10) A = adjoint(AggOp) - n_fine, m = size(B) + n_fine, m = ndims(B) == 1 ? (length(B), 1) : size(B) n_fine2, n_agg = size(A) @assert n_fine2 == n_fine n_coarse = m * n_agg diff --git a/test/nns_test.jl b/test/nns_test.jl index 9435d1f..5466c38 100644 --- a/test/nns_test.jl +++ b/test/nns_test.jl @@ -3,7 +3,6 @@ import AlgebraicMultigrid as AMG using Test using SparseArrays, LinearAlgebra using Ferrite, FerriteGmsh, SparseArrays -using Printf using Downloads: download ## Test QR factorization @@ -259,7 +258,7 @@ function create_nns_frame(x_coords::Vector{Float64}, dofmap::Vector{Int}) return B end -function cantilever_beam(P, E,A, I, L, n_elem) +function cantilever_beam(P, E, A, I, L, n_elem) le = L / n_elem n_nodes = n_elem + 1 dofs_per_node = 3 # u, w, theta @@ -316,8 +315,9 @@ end ml = smoothed_aggregation(A, B) @show ml - @printf("No NNS: final residual at iteration %2d: %6.2e\n", length(residuals_wonns), residuals_nns[end]) - @printf("With NNS: final residual at iteration %2d: %6.2e\n", length(residuals_nns), residuals_wonns[end]) + println("No NNS: final residual at iteration ", length(residuals_wonns), ": ", residuals_nns[end]) + println("With NNS: final residual at iteration ", length(residuals_nns), ": ", residuals_wonns[end]) + #test QR factorization on linear elasticity aggregate = StandardAggregation() @@ -341,7 +341,7 @@ end A = 1e-4 # Cross-section area (for axial) I = 1e-6 # Moment of inertia (for bending) L = 1.0 # Total length - A, b, B = cantilever_beam(P, E,A, I, L, n_elem) + A, b, B = cantilever_beam(P, E, A, I, L, n_elem) # test solution # Analaytical solution for cantilever beam u = A \ b @@ -351,8 +351,9 @@ end x_nns, residuals_nns = solve(A, b, SmoothedAggregationAMG(B); log=true, reltol=1e-10) x_wonns, residuals_wonns = solve(A, b, SmoothedAggregationAMG(); log=true, reltol=1e-10) - @printf("No NNS: final residual at iteration %2d: %6.2e\n", length(residuals_wonns), residuals_nns[end]) - @printf("With NNS: final residual at iteration %2d: %6.2e\n", length(residuals_nns), residuals_wonns[end]) + println("No NNS: final residual at iteration ", length(residuals_wonns), ": ", residuals_nns[end]) + println("With NNS: final residual at iteration ", length(residuals_nns), ": ", residuals_wonns[end]) + # test QR factorization on bending beam # Aggregation diff --git a/test/runtests.jl b/test/runtests.jl index 0adca40..fa2cd0b 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -9,6 +9,7 @@ using Random: seed! include("sa_tests.jl") include("cycle_tests.jl") +include("nns_test.jl") @testset "AlgebraicMultigrid Tests" begin From 959f8847f4935cfa66eadd2785f97e562d21906f Mon Sep 17 00:00:00 2001 From: Abdelrahman95 Date: Wed, 28 May 2025 10:51:10 +0200 Subject: [PATCH 14/24] undo formatting --- src/aggregation.jl | 44 ++++++++++++++++++++++---------------------- 1 file changed, 22 insertions(+), 22 deletions(-) diff --git a/src/aggregation.jl b/src/aggregation.jl index 640a3a0..b5c37cb 100644 --- a/src/aggregation.jl +++ b/src/aggregation.jl @@ -1,20 +1,20 @@ -function smoothed_aggregation(A::TA, _B=nothing, - ::Type{Val{bs}}=Val{1}; - symmetry=HermitianSymmetry(), - strength=SymmetricStrength(), - aggregate=StandardAggregation(), - smooth=JacobiProlongation(4.0 / 3.0), - presmoother=GaussSeidel(), - postsmoother=GaussSeidel(), - improve_candidates=GaussSeidel(iter=4), - max_levels=10, - max_coarse=10, - diagonal_dominance=false, - keep=false, - coarse_solver=Pinv, kwargs...) where {T,V,bs,TA<:SparseMatrixCSC{T,V}} +function smoothed_aggregation(A::TA, _B = nothing, + ::Type{Val{bs}}=Val{1}; + symmetry = HermitianSymmetry(), + strength = SymmetricStrength(), + aggregate = StandardAggregation(), + smooth = JacobiProlongation(4.0/3.0), + presmoother = GaussSeidel(), + postsmoother = GaussSeidel(), + improve_candidates = GaussSeidel(iter=4), + max_levels = 10, + max_coarse = 10, + diagonal_dominance = false, + keep = false, + coarse_solver = Pinv, kwargs...) where {T,V,bs,TA<:SparseMatrixCSC{T,V}} n = size(A, 1) - B = isnothing(_B) ? ones(T, n, 1) : copy(_B) + B = isnothing(_B) ? ones(T,n,1) : copy(_B) @assert size(A, 1) == size(B, 1) #=max_levels, max_coarse, strength = @@ -28,15 +28,15 @@ function smoothed_aggregation(A::TA, _B=nothing, # agg = [aggregate for _ in 1:max_levels - 1] # sm = [smooth for _ in 1:max_levels] - levels = Vector{Level{TA,TA,Adjoint{T,TA}}}() + levels = Vector{Level{TA, TA, Adjoint{T, TA}}}() bsr_flag = false w = MultiLevelWorkspace(Val{bs}, eltype(A)) residual!(w, size(A, 1)) while length(levels) + 1 < max_levels && size(A, 1) > max_coarse A, B, bsr_flag = extend_hierarchy!(levels, strength, aggregate, smooth, - improve_candidates, diagonal_dominance, - keep, A, B, symmetry, bsr_flag) + improve_candidates, diagonal_dominance, + keep, A, B, symmetry, bsr_flag) coarse_x!(w, size(A, 1)) coarse_b!(w, size(A, 1)) #=if size(A, 1) <= max_coarse @@ -54,9 +54,9 @@ struct HermitianSymmetry end function extend_hierarchy!(levels, strength, aggregate, smooth, - improve_candidates, diagonal_dominance, keep, - A, B, - symmetry, bsr_flag) + improve_candidates, diagonal_dominance, keep, + A, B, + symmetry, bsr_flag) # Calculate strength of connection matrix if symmetry isa HermitianSymmetry @@ -70,7 +70,7 @@ function extend_hierarchy!(levels, strength, aggregate, smooth, # b = zeros(eltype(A), size(A, 1)) # Improve candidates - b = zeros(size(A, 1)) + b = zeros(size(A,1)) improve_candidates(A, B, b) T, B = fit_candidates(AggOp, B) From 0a822151e9d133c4e5007341fe8c510b47b2389b Mon Sep 17 00:00:00 2001 From: Abdelrahman95 Date: Wed, 28 May 2025 13:26:29 +0200 Subject: [PATCH 15/24] rm Ferrite dependency --- Project.toml | 5 +- test/lin_elastic_2d.jld2 | Bin 0 -> 56535 bytes test/nns_test.jl | 151 ++------------------------------------- 3 files changed, 8 insertions(+), 148 deletions(-) create mode 100644 test/lin_elastic_2d.jld2 diff --git a/Project.toml b/Project.toml index f690ceb..6891f2b 100644 --- a/Project.toml +++ b/Project.toml @@ -24,9 +24,6 @@ JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819" LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" -Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6" -Ferrite = "c061ca5d-56c9-439f-9c0e-210fe06d3992" -FerriteGmsh = "4f95f4f8-b27c-4ae5-9a39-ea55e634e36b" [targets] -test = ["DelimitedFiles", "FileIO", "IterativeSolvers", "JLD2", "LinearSolve", "Random", "Test","Ferrite","FerriteGmsh","Downloads"] +test = ["DelimitedFiles", "FileIO", "IterativeSolvers", "JLD2", "LinearSolve", "Random", "Test"] diff --git a/test/lin_elastic_2d.jld2 b/test/lin_elastic_2d.jld2 new file mode 100644 index 0000000000000000000000000000000000000000..02e07b99bc5ed8d46c36406da3bc80cdd91b0ada GIT binary patch literal 56535 zcmeHw2Ur!!()J*-rIZqSg))5-{%f9rPnR=_LtE+2zn$cpYd9SW>Y?s>GX^mXA z(7{$~sI8l=mbr_ov#nbft>yNvOC4MmX>~X5VceakW!QL%iLr@^ad)lWrd{VaxM^7q z>%{y0`Cs*br^QfnC3#JGaama&PrT8JFYrl7ZETNEZp5hU7VE}X3>_!V+a@uH$MhDS zt|dIO$ky4O$Lpsi##0rxx(k2Agg-ohyg@+2r$HGd;ZI~kL-=VF+W0&EzxOdY%f{b_4;f?m#~BSA z!jBl7-yoskGsMN!-nb#QCp54$NN<$a_K%}H-gu#|OYK=#SKF0KjT;+Qwr;Kt?nA~8 zX*fj!&R2x<|Bz?C(8boRw`s#yl5kx9UyqMi1V?#1DL4-5@iwS6o|p*;m+=@C{OvdX z_P}os{Pw_a5B&DPZx8(Tz;6%y_P}os{Pw_a5B&DPZx1x(fyQ-;6B0JUbqXmS&ro>2 zvv6@jxX2;=PxxVYi+HmC+MMgMaIu>!PZifETwPWyw_V6q>~MUc>r=7O!FL0WX&j{{E*A zE-w7b0kO1$Yrcb`Vfg|ngu49s^B*X>nhqChq?rAmtVo1QWV1RMcq1K6;p)bLr>2)s>z`v6Y^UI!EaJ_k$yd1F#cdKfvLDHh{K( z&VVZbHvsMgJOp?a@H(IX@Ht=_U;$t~pc0J2c7R<0`vDFEv<933=nS|Pa3|nVzzcx4 z03QR!0KNe%0IUU+Ysgy-{R^lM*bT5B;1ED7z$t+CfG&V50XG5e2J{9z4R{qW1TYNn zIbafC7GNP@HJ}^}eGR~lfV}{R0onlC0y+b(2J{5<20RCN9WVqi3h*`Hd%$AAT0kk8 zMt?Q~)CTMZI1q3+pf%t$KzqO?fNKG_0PY6#20RTI0C*ShG2nB+B)}}dQa~};Mtv0k zTL5+d>;%{ka0K8aKzqO?fNKGF10Dyw0C*ix0Qej*1uz$|98gTIQ9l*Hc7UA$`vDFE zv<9>VTmrZra5tb2;3dFdz%alVz%;--z-mAlcrefcP#4eya3J6a!2dYQS-`>nJj+$V z1j>wSbip5wTi{m&qZ2EQ|M(r!^&cfk!xxZlgfn#G2Zjv~!jz{E-x~mq2b=&n3vf2z zOF*3E|9|IG0)3SMasJ051boZWgyXdRli>TwfZuEHKh3WIdN*b=Z6U~9mRfQEqoYy7l5bbp0`UEzRF z0iXSs_WmfJ6!2nx`Cfn1@k#f88<3CBcWHn8)BH{4TROgJd+Gc{+l%)t?ZFTKG+$Ht z)A_Y2{eM&*9dES1>3E{^3GHuP@UtG^ua@7GzI6Vh_3Z?D8v%9(>;h;4*aMKZ@1N$U z`-ir-Bj}C!owlzl9B-_Scs$y^pDmxB-{|~G>)#9X?G0!JI2w@7_q07fn~%;nbiSg; z`y|jC>n8_jI|BaQ^62qS=My@f>G4kIrw*W}79c&oel;I$5AFZ&wde1a_r3hu;8&e* z`1kLYNB3t_{M|sW?tpZ@{JZ7R@k!?=+8^J`PwW4y`RMx#x?ZF67d^h{dhaRp3)Z_Y zq5ZA<)1S>p=VN-lYAXMC1^scp>Iv;XHX{mur{{-X%}2*$C(!R(^+i+T=V!~K{ojupeN5zyW|}fWrXI z0fz%x0RGeR&Ve51|4V&-^!g#dI~0)i7cKt^oag^v>ieU7-y6Ry=s6gW9*=Z<|7`hu z(1QV_^`rem>o*3j8w(f!7zlU^@HXJD=A-l9_wv*6M~~O0`tP6Sr}NK`+E0%+y1#!l zADu60{-)x&ss8=d@@ajW(vOZmTHha)N9P|po@x0_jgP-uUQ_i#Q}IscGrFFj^V_fH zYXx}|pFh$3->ZKZ$iwG55zvkV{MCGP|9r1K^!#!c?6?O=+ZP4jKLPxw`K>_@d=@_* z+W5TO2HFz=CjsL6=X_`{0CWJv_0Xlzb^}}i=nnW#^NWEVu)XR}`wMG=KG^6Bx_lz(S{y)yxw02cx} z117g-&DT&UVqTz%O32W2l%V`zE^)+&~pypTtGWO+8(<9e>ESSzi9u` z`Jt))5rCZp@Fn24*1vx?A3ff_w|{nm{<{FZ0QUlNfR6zGf9Ip)Z6)Zv3UC|Xc0k(x z@3rqo`RMqg`RRD1`+a=@nQu^*ME0X&+3w7&HG+Z5m5 zEsxHZbbQeHmySnzzNN>*ua@7Gesn$1l)itrJlbEhK1(1@TmhE?;{5~d&mWcdz5Gp$ zpTAokJzkos*S?p(sqz1_`RMUZ_kUCV{?Y4c|28F$j>o3-``Pk;G(T;F@rCpC-z{%7 z$X^4v7H}QldO#1r4S*W~Hvw)2+yb~2&=YV6;6A|pfd6a!v^@vFj)Q=Q01pEm2Rs4z zz4rcUKDz(BLGL4gK7hV}w0$Sx`%{45YwyqIZz_K2_&f`C_yJx8y!H+I{%QGiyzB-& z_W+&-JOfDENB37#_Wa#^biSeEVJGPQt@|~)|LFLo?fcn$biStRIXb`7_R;#&AF@OR6j{XzSi zj(6I>^mx7x{(k_-0pj}1ujZrUcLV5y^$~4vQ})sG3vKVO=A*~^_l|G69{Snx==?(S zABFxo21wU0w0=J-uc`T^seD7{7drp{Z23+3&j!efru_Az^5}S> z_5D%%f3|#DUpgM?_@n*Zlz)Fz9z7p@Z#{b*;`0VzFyI})pDmxB4`{yApvSlBae6+u z3Fij^()RpnzJ8!bf4~8NW`O4aalhyVXkP@p0*L!F|899r@m~V|%iqw4j-S6<9_^2& z^a%hxasMQ(9~~dJ;C$R)Nyp1S%}>Y65YTrh;E(#7_SdiGn)abpE65Z7QE}AP?6ooX|HtrsDBx4TXMnUlFX8)Wz<-+G7W9|{NZU`x6WyOp`G>auSM$;O(ft($ zc83E-07e4R@!3@W{%Srs&{H0;6<}*XE5Ol!-)lc@@2}>g^&bOv;Qal)_R#(Lv*pwJ z()p2|uV{Okibr~U{c1is9`vDq3;@4(Jkas>669dM{MCGPzN7n#wy&xAfsQ}A|9&-} z82DElknX=(@cnE++TN!6>u1aV_jgJwHok>f^nS|bJTVh3zUWokrV2dMFA~C6Z>u-F zjrr8VMw|Yx_f&=$tsd~lS&lyI{rfTctV_(PTNjR1Ft#HmSbrI;#vPu2cj5Sj&zXp$ zS=|@Rm*#rkN?P%9yEr%XSddud(G;fU&gQVN@uf`mfs;oLJxyXVtrq&aOnt*>uGnVQ zr=phZ%&eHI4d*G4TbhPvFrIH--81_Dsvq-y`AkX5YMcs;0qnO!X&q3hf zYMytC7^Tmab+fzIVeC8N>#M0GF(H8jMO>L@bGe+^TlC2D3qOw;^X|NsMCvONUf8}@ zyI3)9U9|l5o|7txzLRv9!&~II18R@o@2f9mmb+%9hZz--l}lHK_ncEj_Soys>#vwZ z6qJ+gq$P@&1C)KpD=t&8KsrZ@MP3^t>keW!-;!~_+uGaNFf05sSr5G7@#`nr>@N*`tsv}j(kJAknG~*_Jwz!hJ&YSe2 z`U~xu2=<`dQ#O?g;>2`WlxtEbyX5&tG0q(1TAtNZ_D=mm#)4ecr>CM<9myeeBWHi8 z%hhF72TPjzF1yLRc$M~`&4DOVWKy^&V?Y^6=0_17^=M|m^o`17bUE{4%xR7(@v1Q9Rp$OlgaB_B)AI4XN62$Qzjj@=p zj-N$7J}kI(CaRJ-uO9X4c8Dss{`Taq4~^nUcqbp$(k+Zc4fYZ^94iy~RT=C>Y;Pe|017T0*(I^Y@5S*desvm_BuY0aRK ziidgJQ;4_I5O0alpXk>aU|$9B)|cMdn46%>cBu~xogk=Sp1kO zmVS`<@MS{@8GY#1o$LE6NqlI(iTgUraq=&7NX)Ea5?GLYZq0*I;`?NKjk8xZq2m$b z=q~gp`W5X(-c3^DHaTQzv&g&m<**euymyIsHLE98WUZ4G@s_`NtWaN&OniXXXYKqC zo27Xi<{KRMINoqv(DqdXu5?a{H(+;2UR@EdUCOvE95uS@UKuWZXrG9oew7Sw*W zFK%S3KH0RcLpm9guAW?UT$;12REo4Ptt6P&Fux!Vj=%Td$FX3KC-CU(9GH6gx;C2z zJXIlOr7=!od7E>@WVZVYV#5pFvO%Vj(D{OnAK#vqzIE;OS=p(9eWJxaGDGc^ z$E79;m=jBCJ>yLxn5!rD?@Y?CCzm%%GwW7WGxq~erG0#=!m;-Aj@TAeky}(gz_`Tt zkOsLceGaHA&(dL0Zs(q#eM9PGIh4DM>tj4QqlBPbpJ^pZuU3~4I=`7v`qeBOUT*BH z%ev%LIP9PBp1iBe&_1wNftwO0=fA@26Ql50#oE5ChRpUK(5KDvO6EvepxtrB^Gt<` z$?($QrWh#T)$G#`!z-y&-J1SYOQp6jy)mf`;u9%4TYR4`|BreAg zXUONLAg^YCeN%zwTDA0wNJTyNfv5G-a-DES>)!Qy#&4CmIp9BJ_u9;oJEA{nR zBo^eVZg7{rGgOLu4RVbis^_f^=5aXx;5>l&8}q76%TALjx9GBxr&l_sO^PBO(Hk?g zRcx8;tJ@x5*_cHRd-i_0aHtp;K4n*U(Arw&VFg#M5GBKf-<~nuM7fAqc3Yx5sOA;( z5#vXXI|^W2M+DhW;21O^-BjqMv<(Dr`y520^Ma51bGs= z`@Xmy{ES2bkKd5$^{afvIJ#cK@r~n>bx+&q{-`-CEErc>RCzETl_cxSCP;F{TJzk; zEUG5k)8_375=&=dtas$hIe3p?-G%iJ<_R3%yj_iQM)r+#dsU^%mVlhtF$T68Yzf)a zFGI{^SR}LmtbNE{x1vAI4<3*|F|P%}`1XhKZPv@8_}OS(*3`Sb<`(m(%oG0utsLe{ zaK4jG3`kcxG1IeiZOnn>_|R)uEeu)%5#7J^;_q*cmq!hv$-oz{4B*#pz+Z1UHt-|I(08 z(;*L{{wSB%Um_`eMKgARu4jDpkUPxI2_Lq4ggjxKFAZK<8vKC4d1NZ&_W_WXF`uHG z8qg2rr1fucIeeNLdq-*O-32`($nsT7o87UkAn31buwxJKJ%ID>?&zy~@mMp~u#JJ! z_)XUsO^vat%SSwA%1=MD6nA>SoR^%lJbtt!_q?q4*qF_g49Y=&ECsz+0G}$zQMw~H zpgUiUbyefWyqgk1Y|84U&Kg!ho`IY_swx)db{|M9;KTDaW*Nw=o8N-n(og(_^2szJ zUop63UA-9BN}(iz`&>>YjOn1^;Tp@Fp6-5Q8{KNsqgNQafkJ`Pc)4c%|E^{ z|KPlV`5SpK{*Y%k__HPW3GM3xJiR@;$SPE8vtxm$@6v*nrxFyn$#B0}4fmUsaR0%; zeF9zIVtt2oR=%C0yIF!B>oA=Zj_jSySQeU_s;-gYMow*8uKhTiD7NqO+~GEC5# zGvZoJ@P4x!+;3V?_anJbM`L|^jH>4xf#-SuOBUwMbl5S#bKh6}(wR2#4Du*GTf8GH zT$U4q`wgu7ael&igcsYG?=cUrk-2EJ{(~-ipvcr#YF|EM5mj?>XTocee0Yb~kxwNg zQM%WQK)+Aq#PN-cl6Do-QNNTsH&UKUcXSB~i7qC3)cpm{gM(nc%7T2a0r~zA@SIyQ z`?Py|9Ts`&2d#7L4kI2wrYAps?7CG#Kv#Ij_^)b$=zK^jvbbyF)lLgkaQ;U-=YW0@Ag6Gj&6L19>g@devFE$HMUu##dP~b9V;H&~q5B*89pSoxt(i^B zdZ@FcMR&7PC0Ce}s%_rywXPv9?PgWJz8fQ|Yq0*n{`?I6iFP1A${AKo1x0*yu7lMcpgvm zJifA$V4lSMCJX(J{e{=P=ShX$N#wEH;kvYelFRaj6%(vSO`#s8^DX8t%qtjAg%Ag* z&km5gxYO2xv{YSIA~9mimzomN3*=shdXvt-D7PQvE0l}zMeB=lrJUEc+UKXsCR5`T z`wRViq_^wx%)!lAoguR~8?TLIE;#W9#AswP$zSBW;^w^|PbTzM>(H`_3@nW1d)QQv znH|&GES*_S$Vv~*W_oWKdVDU1IKcixxqVz^uba+NWl^qEv*Du#sfu$b_rWc@76~st zk%OC)OLKWo8H226V@LH2Vd(nPfXX+hZ#?Mhacfy*a*Q7P=t@BCmW%Jn*3)_YUMy1J zPF`5DAu7I>@e{ALc@-_e^|ISi?=`=KB&$5#bWS>sNm(*f!!rLh!8#Y~Pt0qWZ&HB= z^+%q=&+D0aqqW&nz_Z8C;lipnRRnobXUw?SK9(cNz+*aM;_KTIGF;tmPfPwy74}-% zWw8$SPnbUzO@;rPJL`t$mGq~?{()E_b%l`W}!oNpSM$0-WKg2G8K$O()Y73g5xnniidrs_O+W)#gBftcEfZovONepk{{{ z3*(gOdsEF!$GU{zaj&?Enp_6ZXQ`;@D*=Tf`T8qbC2oCPNy>ucr!J0)XY#p0u`Qd$iTsTI zK|a(EdCWGMp0``4%lZLN+okzoQ^zL|Gewba=n$e8(xGFPwu+A^uR1)26AzH;W)tfJPUHCwGN2L`l!o30=X5@Q~EA< z2`B48uK0$^69FzN+;))bH>8uq`s?zXt*O)NE!w(l#(>sF3wU{qV$k41mAn^>*zCAf zcKPqfmv2Zl0ia%M{;KP{nhFvWIFJ~?w6}6cYG`2k-WRVC;x&X zhyF&t$bx<1VBh|pCBbbU=&&1KwUX<8E|OeqJ2gL(C(do0^;u?V$8zF2px2CrvJ%`- zuLo6u=0!|;`rtOXS6`6xMx!3uTqtJdS#uv8(|V~GTL<-4RtCyswg&H*o}pLVA* zK#D`&y*byKNn2Nv+rT?QAz*N)4TS{9(^-gfx<6@qPeng#e|xJo>n85*GEXL-;g6Xc zbz$EJqA$11f8LU}B)9i^x6BY}Zp+sAxg)Muk{xor4UaBQCKFD(J<8I^5RKb*FrF}O zjzE8+eI8&R^2~S|?9|U(pVa{#<97Z_OVrf4&A^lWaDLX5i64m#@PyTN-g3rLg_DMS z%0V8a^99|1Xpb)NL~je$6a;Cr7Dm>|;S;i%_!oVvE!aA8EUomFK|(E|^J^{SAirn37rbF!r}GsZ|ELejDPC}Dj$&7J_H+96x|nSrMRf<(2blLHAU|O|dqF&-U1ggZ z^+7qZmd&fj7Bpw+dY}UG=MKnobiTsziR0wZte~_RCv?~w>mHB!nEjF&ZZ*sBtdTgE zy88U2RO=c>v-5I|3oFB!Q8u09N)M#|sUE;QjrkPwBIc_*)OfrGJooL&&q&VDVbg$T ztodM1CxWI1}h4(pvHs6)&teQUgToT-_q!`fwRT(ErPD{@}H^Pa5G zFJy4Q&GsJaOG&@Fop)B%)sp@mD<@lr)iOahed5-ieoMY=>DS47#aj{q>w+n;UY7_w zJ)qt}zo5RzlNH-lf9hmi7I}=F7wm}(2x0~RkK0bgxznqv32$p7Pm0Fl=)&PD9L^t@ z*U{hT2h^Vh{U!aTIo$6jeBKxM@V-G_Jh6OoWz);ulH8}UWtw}>Cz7kOnnBC2yd$lr zZ<%%KQ3cUmyXE1r;Y!@mkCiaby=zye2fSc3PWrp-TvyECyoK|L73F`lr#IN6 z2)wh6j*_Rb+AQ*(?|Xc|{D;>J@*cdx4D*PrCr5zSa&`y#9M3X>ywhgx%Z}#rxOt!Y zpPOSVtV0e&Fqg(s-?{h0m3MS0mrvuK(U}(i(Ufc@rbq(+MKp;z)=3jjrWH9(Zr(xn}LGljP|A z0rtkp7@Y{iHRWVxL7Hy2S_#fGT1z>5QMqVef-LM8z~}V%yc~Ja zAIKZLD)Nexp>RD4c#ZSMI&PaN!_@$9^J}-485qAI$ZK}w{;j-iiDU=xK51*%W8%1U zg8i!we%1ziP;OMvgJP|{dTei!ba!>17_zb?dHrkaFC@N1RW-ds4Y`n&+jV;@;q#g* zUSR#LD#mb|_9O%QY*P08Av3H~IU`(u{3C9#|IkmUuPMmQT+C@H_Sa!i?!7F#q{${# zWHQK&cPdGcuPGJD-E~epdGt_O&I|T4$isdHT%X5v^=S~VVh}gjKV{&jBo8Ua#fS9R zy;EF5rIKTaM^_Us;lDP?t#Nui;SE>AD1uzvuYhvnV7;D>N8awnenPj+JjGCD%P4H z%RRo;d2Pnh*9`6N{=h5wdfvd;Yg+6o;N2U4m1ia{&UpfFtlBfFp3cc+Kk%m7gwHq@ zm`(;a>pSD7aNh$DKHIMrAGW(t&R{zfc>AtmtXfQs_|y0}PsNiv*pGU+N7UX?*VSW> zsYG^PlM_z%u}d<>ZIj~WxT&`7@uY@yd)n;8g+7lNPQI+%KOln9$Y^(U{thWly(Tne zSg$e$$HM>^4^tt&UPFALpS{4|^T2c8Dy7Gp?bRlJ zN^$e7B>BU2%HnTS6Rj&axgl>QIk1MiKU;0C(E$vP8N7B7;T+?K0ulK z(|#(*Zsk2nCcG?;F*5M<3nDRE_)Mr6g`{| z&*2vkWWghKNbVerHHE4D(`gITfn1`JZ9l6NqO!p@bI)U z#0QuZig@;~2+|m@T}hBf>Vv4=I(mq@#6m@`ReQ%U zxAbC$j;}kEJtu)@{*kFItM}=!w}B^T_QfvqJ>C-o@47Jsb4wZId(>~Pk8?ofS3clz zyU{B6gF*&rO^tsZHGb&$!~R3RTDd>ll;WhrvP^3Ut6swSp{loulf^?OFz3sG=*MrF z=+wm7lB%+t-OG@yfOVw|=9w8V{^!E@IRfJc<6Ify590*;4|#?UF4R6TU5E7tp1${5 zSQy!V<%v~X<^NuoZ|VALT&;?i_`Mcvim^9)a6=9=)_qMg$NCV4eO~A&dE_kP!_69k%AEaPF&J2gvx>it4yFP*0_ToQF0h4DXg2>lJ(};M#;};#Q{6#nH2h z$n*O=dOQ7&$p0uOrPYYHT@-3W&&M~wbL|qUKE(NUEzGx*Kwsopr(dNt7wdOy`MtVmqI@qz{YB>rCY~tH8JqmoUkhV1 z@7}j8BgX;~Z{!AlTfe`m59*Id-#^jy^$@7%IH;@8UwAzO`eXaZT(QMg4|Unp*seWm z$Apu#%(`d;Sv4;6?b^$Ck3VJbJ`L}eN}-O&dYi7_=zc*zAzvKu%u7DtmRqIGegK|j zyG~77b~lY+K9YhwLi;yk%v|+1m&8T=G4bUA{xO4cd_Kn@PgD3U|M!LaQ5m=|)PQ`r z9P%OdFUn5^e`SLms1F|Zm&#clx>cWD5Mt5BadZhI2XceqeiQR!GUP`(Zm=KR!9VbP zjzRhJLGHmiHsphy4yz3FFV2gYe=twb<2zUO)jM;sukl?sx=by-wu+$q{}#XH-#Vp- z<45xJ*hTfZ$GmDY7+qz7y@Bv~^vY|V*K;MSiFflakNg&VVEpF3OWig+nR!1aH}TSn zyQJdPL;bw47!oF7wqINC9jS7fYnzZ(!KhnVa7p4aoJZdlR()ERGACYYTCJ*3;Kq)V zT2V3f0kKr-uv34f1hvs^{8i=(EUE;d#g` z$sm)20gr@Y#$&U=Djf2i*?&r^w`3y0_D0tu9oOwkB$Bp`JSr72&PVpEaYv~BqwTR; z+gC7gx~}kf|B8~xlXsZ>t53!xEUhDkwVNdGoS~y+nRZO;LjAIK~b3AKF(A_T2&A2g6pxDYVmOI|8o` zyUuXt#}Z-zyn7?skJfbCMuwQaF`B9J_5Kxkn~$@9!)l&m&XzeWQ2qS%ewL2c0lmz7 zD|FUpFLsM)R?_Dhb8wxJqouAo=kOwSo!Iq!a?xm5?qLIYZuz_pOQg0IG0(a@D7%zW zPSTgFRK7{BVTO9&@Ov#^PTHI9=Q$j!7L9`q5Ra!ou8Z_?tLC;kEXv)hwSTv#dkI0g zs?kOR(%PnokpFYcvOB;wiOSr${D z&g1BOkw?W}#Okd#oBQgpcO1suU7z)k8C!UFca5V2XIL<2zvZrCX5;4%y9?ySx&Eb> z+YRhnOV)JOUBApdjtR8#a_hP%hB>9QJnE! z{N#r$==L?==Z!Jx;~ARBoX{B^zHweOGa`I?tN0t0#A2G2S*Uy^W4LlrOq|PCzoLJe z!f)mOrmn;MZHf6C@*n0Uj5mw}rrkS}LDf3!MC+@)c8+<+p!}Lcm;Do;3je+UK6|h0 z-oCc>+B)R6pvfL`n+tG^y&-h;Av&;(eFJBImzh;#XMv z^l_8t4)AK8J*lW6$bo|_lW|sob{uw68eQ!&O zc=q(4pX+n*>wPs|PviIYJLU--zu0fIfA62uUR_nFD)N6{y`xi29(}F5hi@o5XHX*C zZ%VDF(Bl>TiFPLWrObP*qst~rSa!6s6EL4{jh~??kmfQ;!sWB`OPNnfd-@h`<(P$s z90vGnzh?HzJ|A~xk8nMzi|*RPm#hD@ovTQ+nWArjx;b z*Pfx7jK!=2x6i63GvlP+_I&bEitAjP^KskMx(KwluPqHkNi}V+4gJR#`C5I@7UStSQQg`ZLI=7xkyw$4Tm}g=cD;B zPmF~8j@M&8ISBbg{h+Y(0h!v7>koV#XCpwa1=L3lAh!j~KR9oo ze3ZKz@(0Rgr#hQ-*VJR@9rDtdnEY+Oey=}~m&U&;2kXtjRNaU2Bu|{EQ0w+JE;c+j ztxdgFOUysIgirlc%Zz=IYNtHzYu$s_HHCjwKR(p=JgB>|9$pRgF!D@?xWG0d>Wgi^ z{gbo;>a^LYsm21SwM9gXs^?dMTve!d@j8q*Y$FZ>KVTc>@&_7aEK@OHv0lfznXWJT zL%d;J(Eij>80oU0`B(oV&yV6?_5WWz?_2s{9gOuUo`>@nw)yW0x2FtLWAjzb9=uEN zW=@>GG^Kw^41>qtyt!(;)3}uMZ*`~EO_(=ZtGp?lz?2gD{y~qbxAFQiko!(y-Qv*M zI&9vcfjgwEK9Q(Z>y*FP3-`-CG>U5ZK86`~EHCiroETESSZ33l#BwtBaYRMY`7#o^ z)7xc8eKmvkDR@7Dbv@SG$afs{&jy}~&>HihzPjve;OTtHK;r)7AO?9l86<4BHTz7$ zfJZPor_v4 zZ|Aa5bNNhy=Og~N`r!45$cOFfb4%H`FVt8or^CzI#*{P9I#??jj!z*e(`-9A3ZI)5 zKD!;6WA^oatc!#~1|VTR%`Qpy2EEFIPJ}9@MRzJh<8Vw%dpgWWv?qRuh*A{czR(MfCB{v{)x zdP|7cD#D4dBXiKxDOcj z=OPd8A4K`{srbSELmu=8@|=0Jy!+h|+H5=Iv0Jr%;n4R4_q*agRNN1V`xNQ^r1kx5 zJ1Jp&pbmTH*1j%Z?xhKz#~WC<%`GGktyW*~38-Xd!2Z@4*x!o#EOGxM_SY%sFVrs& z^g97^?K2OCN=?vX1t6F4@nY3Iq&RxL!w~9ItQ#XCKVjZT^blC>FVtaYn`KQl%el{7 z^mJaF)U}AT>_QH%sj6aj4lCRhJ^mT#KlGqhSnG14+;QTT7S$C*+2_Qx4E-`DvRdP6 zpDWT_#rTB|=2qFFb(V5i53!``)#;G0Odwy41zxu+_e?Lo)Mk;lRPbt-!ZaBUdCwm& z5`UWUMZ}xG-L>u0s7p*0@Jbr1tx0M9nqWNqRs5^@1?LNt+ZFObG>l^$Pk6i)o9ezy zNu6~~`n=40P%SfV!Lz|Hb*qUv%m;NaAJF6f66OE0Ks~QeD;-vH+)+ah*DBD>1nVMS~iITNP6R@SG6@rHfqxIY~C z?_GrbdzeQtpW!&<(t$_obm6EE;|*BkDf?18r!Yxgls~aO)_Ytzo2tqv#Yvdd3P^?h4VL@<7^9>F>R=Woo%s*sQIe2g!Yi~Zvcat-d7 zwu|r*?z5xzQ{wm@0OK3u2m6)w_ksa_A)dwBY`4k}g9nc+CKEI>KVIokOe|*Sojcf~ zikalN<+}fJSFVX`=m;2Vnms^7z4c?n1>6`nw$bjXWXO&hN6z*JY7sSF1T! zD{|$z5y128gZXiPv&cWKH{kua5!{c{`75>g#)LlNI&3S4DSI;NJ}}N=U)qa%z9Z6e z%-0yyydam_yy`9?eDBbWiE;a%6THQ7U*z6j^@p#w8VPN^oiu7T?fc;Aor!XMR^Kodm(yX}_c2MA4vS;C;9rjrNGk0vxMUWs*$uYWpN=g3e^YQY+_iJs-oJKKFjqiACqxc?W+cZkj=Ll^in5q#E6$u!c`D8IY{PEuP+KC}pEMTs zCtaoXgPZ`~s4Dk!%B!_mu*n5NAN3suixrb_fuFU6g$eZbZaohBYdcp@@+2l7` zSDO)nyqBu%ZPIOX$wk0DA`q)K?Ka30XM;`d&`EplBRZSgM=hDL*zmjCq zqjaeM+VxM!n*6h849&6_+y{mGlW<=R?svg;I9wmYILG*-?MEKNVQ-!UWazMKfyZos zrbgNh1KIlq$h*B9LpYxJKpT!ct zS^8-=gUT7laEoLuU1e^{k#WtQ+s84V7q51oeOHZ3ICEGcK`NF+Pi?#W;MQ(i(Wr(P((zMijgU1NYxwQ*;JGKe`&g%*F_W%=NKEiooIOM-4kf$)-Fb-&aea%Kp zaB8K;Iz<no9hi@2Rv9G#0BK(onRjV)6PF>!%Oj^R@N+ zw;%Svse|iRGr=br6y(t&d5WyB&)2rfX}NP5M?`IF?M~*4GWLqvcpaj-yZnwf&X0(G|qFg zrG)R3mXH%`E_@%eI8RLY%H=;ke{H_&8{T-fM>DwWKg@S|xiX93Xaf}S|E2M}{o5bE z1#P9L(Z|?X#K7Nmne+!5o9Umy{N9Fv!9zEHJzy0ysAOGrs-|_IT-uPYj{*DL8 zldwcs#1Cj~<4|UmC~8|?9N^TUeVQQbj`ON6=W}7+FA{k640s;*Ig`K7OS4tib@F9OEpEd!DVkF-xG@SE{d~UJ@MVc>piAc@?QFIk@c||89_m)4Ds)1YJ&V zTp&IC4IYQ{dAyg^`X6qD#R)*Z$Cmb|d#>7TwaFC86L*}nB<<;2QCp+GgwZn-^iPps@Xm(g(W2uQy&BqtA3s8~;Fxxf;7hH9 zkz7hF-z8>VsT%iD(0*fsf1~)84mYc}&n^(`TGPb-y};e_PByFs?_Uz_4$SQrVbn{v)@v^dUE%0w0fO zcdVg(U(`l^%fNNvr%F4p^8LSy*5#tzz1?-V2Go|0R4Z%ZQILLvCt%z!)*-wBOHwjonLh}sw1>Z`Yv zeJ^Tzw7j#W+%Zkm-V^eoa^jE-fmSnio5P|wK5z0#wdQi^g7dek{dX^U$#*>X!qr*t zqhRXgNlBB7^ZA`EecHIk=L$4ezqDB&kfi!8yQt`9Rl+2r!~y&L4PV`_#V z-=xNHaaumVEJW-)&m~`Aczi+irH&u@Wtw(*3p!^Dh|<-E=AZKU!&4S{h1(Pe3QjD} zxmT3UKcr_g`3H85aq=EnAb8>9)T-FFfUjrNL8^_xJ3+`;BS}T( zBY!MwwK-fhPoUr|IsK0FTfW_Zd+zzO@&$``t7Vx-zv45iH|C9-mn#SextCe;A&Gyg zRp_b>Hwy%rN6)_Z+Zn|l;j}JoYFWNObNhugXC|lckGfuH8(o|(L8LGeig$9rqC`DJU*NISHP6S$q%Xm6$Rk^kl5{C1lU zMhLE%4UL=p%uu0l%uPcrxEC+0uvY5hjR^KCEv0CYWy$zAp8(c$C2C;5LaN zZ({kEss$$MO7{iIzTV^J7Dw?fui07i%061q$Jxq7{7Dd>c`>N|vv{!Jp~VyHPP5|p z8)G5novx_-WNqs-Ml>7C#af4zQ+YYP9f&6x8mpC<~8 z+u3w}^)!L+@^D4{>6tKtjWoLZj2_vCNi?ELVZAg=A>h?;xv`RB`yTi*|PE!do1 zaJc>Mbbh&=_Jik%se(@@tOr@$OyghKtY0tpGE)$-$$GeqVFuss^5&Gi#%}~OB4>AV zaZ2Fd^?k2yaxh&mbll*@CIMM|yIqGyFV%b}*i#U>C;HA?{8na K*y+1L>;C~WO_VwS literal 0 HcmV?d00001 diff --git a/test/nns_test.jl b/test/nns_test.jl index 5466c38..678146e 100644 --- a/test/nns_test.jl +++ b/test/nns_test.jl @@ -1,10 +1,3 @@ -using AlgebraicMultigrid -import AlgebraicMultigrid as AMG -using Test -using SparseArrays, LinearAlgebra -using Ferrite, FerriteGmsh, SparseArrays -using Downloads: download - ## Test QR factorization @testset "fit_candidates unit test cases" begin cases = Vector{Tuple{SparseMatrixCSC{Float64,Int},Matrix{Float64}}}() @@ -69,138 +62,6 @@ using Downloads: download end end -## Test Convergance of AMG for linear elasticity & bending beam -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 create_nns(dh) - Ndof = ndofs(dh) - grid = dh.grid - B = zeros(Float64, Ndof, 3) - B[1:2:end, 1] .= 1 # x - translation - B[2:2:end, 2] .= 1 # y - translation - - # in-plane rotation (x,y) → (-y,x) - coords = reduce(hcat, grid.nodes .|> (n -> n.x |> collect))' # convert nodes to 2d array - y = coords[:, 2] - x = coords[:, 1] - B[1:2:end, 3] .= -y - B[2:2:end, 3] .= x - return B -end - -function linear_elasticity_2d() - # Example test: https://ferrite-fem.github.io/Ferrite.jl/stable/tutorials/linear_elasticity/ - 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 = 1 # linear interpolation - ip = Lagrange{RefTriangle,order}()^dim # vector valued interpolation - - qr = QuadratureRule{RefTriangle}(1) # 1 quadrature point - qr_face = FacetQuadratureRule{RefTriangle}(1) - - cellvalues = CellValues(qr, ip) - facetvalues = FacetValues(qr_face, ip) - - dh = DofHandler(grid) - add!(dh, :u, ip) - close!(dh) - - ch = ConstraintHandler(dh) - add!(ch, Dirichlet(:u, getfacetset(grid, "bottom"), (x, t) -> 0.0, 2)) - add!(ch, Dirichlet(:u, getfacetset(grid, "left"), (x, t) -> 0.0, 1)) - close!(ch) - - traction(x) = Vec(0.0, 20.0e3 * x[1]) - 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})) - A = allocate_matrix(dh) - assemble_global!(A, dh, cellvalues, C) - - b = zeros(ndofs(dh)) - B = create_nns(dh) - assemble_external_forces!(b, dh, getfacetset(grid, "top"), facetvalues, traction) - - apply!(A, b, ch) - - return A, b, B # return the assembled matrix, force vector, and NNS matrix -end - - ## Helper functions for cantilever frame beam ## # Element stiffness matrix @@ -307,7 +168,9 @@ end @testset "Mechanics test cases" begin @testset "Linear elasticity 2D" begin - A, b, B = linear_elasticity_2d() + # load linear elasticity test data + @load "lin_elastic_2d.jld2" A b B + A = SparseMatrixCSC(A.m, A.n, A.colptr, A.rowval, A.nzval) x_nns, residuals_nns = solve(A, b, SmoothedAggregationAMG(B); log=true, reltol=1e-10) x_wonns, residuals_wonns = solve(A, b, SmoothedAggregationAMG(); log=true, reltol=1e-10) @@ -315,8 +178,8 @@ end ml = smoothed_aggregation(A, B) @show ml - println("No NNS: final residual at iteration ", length(residuals_wonns), ": ", residuals_nns[end]) - println("With NNS: final residual at iteration ", length(residuals_nns), ": ", residuals_wonns[end]) + println("No NNS: final residual at iteration ", length(residuals_wonns), ": ", residuals_wonns[end]) + println("With NNS: final residual at iteration ", length(residuals_nns), ": ", residuals_nns[end]) #test QR factorization on linear elasticity @@ -351,8 +214,8 @@ end x_nns, residuals_nns = solve(A, b, SmoothedAggregationAMG(B); log=true, reltol=1e-10) x_wonns, residuals_wonns = solve(A, b, SmoothedAggregationAMG(); log=true, reltol=1e-10) - println("No NNS: final residual at iteration ", length(residuals_wonns), ": ", residuals_nns[end]) - println("With NNS: final residual at iteration ", length(residuals_nns), ": ", residuals_wonns[end]) + println("No NNS: final residual at iteration ", length(residuals_wonns), ": ", residuals_wonns[end]) + println("With NNS: final residual at iteration ", length(residuals_nns), ": ", residuals_nns[end]) # test QR factorization on bending beam From 3c215204c5d6e97617304acc76cdcb3c8d7cd72e Mon Sep 17 00:00:00 2001 From: Abdelrahman95 Date: Wed, 28 May 2025 13:36:38 +0200 Subject: [PATCH 16/24] undo .gitignore added extensions --- .gitignore | 3 --- 1 file changed, 3 deletions(-) diff --git a/.gitignore b/.gitignore index 93c87aa..ffd6da3 100644 --- a/.gitignore +++ b/.gitignore @@ -3,6 +3,3 @@ *.jl.mem /Manifest.toml .vscode -*.mtx -*.csv -*.geo From 3ac65287a6b636912e1074a95a2ac2c1a73cb463 Mon Sep 17 00:00:00 2001 From: Abdelrahman95 Date: Wed, 28 May 2025 21:37:30 +0200 Subject: [PATCH 17/24] fix ci --- src/aggregation.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/aggregation.jl b/src/aggregation.jl index b5c37cb..517a573 100644 --- a/src/aggregation.jl +++ b/src/aggregation.jl @@ -70,7 +70,7 @@ function extend_hierarchy!(levels, strength, aggregate, smooth, # b = zeros(eltype(A), size(A, 1)) # Improve candidates - b = zeros(size(A,1)) + b = zeros(size(A,1),size(B,2)) improve_candidates(A, B, b) T, B = fit_candidates(AggOp, B) From c030852e5160edc264ab6bcd969c74c72719619c Mon Sep 17 00:00:00 2001 From: Abdelrahman95 Date: Wed, 28 May 2025 22:55:13 +0200 Subject: [PATCH 18/24] add nns to prec --- src/precs.jl | 9 +++++---- src/smoother.jl | 1 + 2 files changed, 6 insertions(+), 4 deletions(-) diff --git a/src/precs.jl b/src/precs.jl index 59b683f..8d1dc23 100644 --- a/src/precs.jl +++ b/src/precs.jl @@ -4,17 +4,18 @@ Return callable object constructing a left smoothed aggregation algebraic multigrid preconditioner to be used with the `precs` API of LinearSolve. """ -struct SmoothedAggregationPreconBuilder{Tk} +struct SmoothedAggregationPreconBuilder{Tk,TB<:Union{Nothing,<:AbstractArray}} blocksize::Int + B::TB # near null space basis kwargs::Tk end -function SmoothedAggregationPreconBuilder(; blocksize = 1, kwargs...) - return SmoothedAggregationPreconBuilder(blocksize, kwargs) +function SmoothedAggregationPreconBuilder(; blocksize = 1,B = nothing, kwargs...) + return SmoothedAggregationPreconBuilder(blocksize,B, kwargs) end function (b::SmoothedAggregationPreconBuilder)(A::AbstractSparseMatrixCSC, p) - return (aspreconditioner(smoothed_aggregation(SparseMatrixCSC(A), Val{b.blocksize}; b.kwargs...)), I) + return (aspreconditioner(smoothed_aggregation(SparseMatrixCSC(A),b.B, Val{b.blocksize}; b.kwargs...)), I) end diff --git a/src/smoother.jl b/src/smoother.jl index fdd6657..829c6ec 100644 --- a/src/smoother.jl +++ b/src/smoother.jl @@ -29,6 +29,7 @@ end function gs!(A, b, x, start, step, stop) n = size(A, 1) z = zero(eltype(A)) + @assert size(x,2) == size(b, 2) "x and b must have the same number of columns" @inbounds for col in 1:size(x, 2) for i in start:step:stop rsum = z From 53ddeb71166c708db17c994859e2c13ddbae156b Mon Sep 17 00:00:00 2001 From: Abdelrahman95 Date: Thu, 29 May 2025 19:38:53 +0200 Subject: [PATCH 19/24] add B to kwargs --- src/precs.jl | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/precs.jl b/src/precs.jl index 8d1dc23..3800ae8 100644 --- a/src/precs.jl +++ b/src/precs.jl @@ -4,18 +4,18 @@ Return callable object constructing a left smoothed aggregation algebraic multigrid preconditioner to be used with the `precs` API of LinearSolve. """ -struct SmoothedAggregationPreconBuilder{Tk,TB<:Union{Nothing,<:AbstractArray}} +struct SmoothedAggregationPreconBuilder{Tk} blocksize::Int - B::TB # near null space basis kwargs::Tk end -function SmoothedAggregationPreconBuilder(; blocksize = 1,B = nothing, kwargs...) - return SmoothedAggregationPreconBuilder(blocksize,B, kwargs) +function SmoothedAggregationPreconBuilder(; blocksize = 1, kwargs...) + return SmoothedAggregationPreconBuilder(blocksize, kwargs) end function (b::SmoothedAggregationPreconBuilder)(A::AbstractSparseMatrixCSC, p) - return (aspreconditioner(smoothed_aggregation(SparseMatrixCSC(A),b.B, Val{b.blocksize}; b.kwargs...)), I) + B = get(b.kwargs, :B, nothing) # extract nns from kwargs, default to `nothing` + return (aspreconditioner(smoothed_aggregation(SparseMatrixCSC(A), B, Val{b.blocksize}; b.kwargs...)),I) end From 04005d5843976c28d6b63d59e3b653c9d1e9c56e Mon Sep 17 00:00:00 2001 From: Abdelrahman95 Date: Tue, 3 Jun 2025 10:11:16 +0200 Subject: [PATCH 20/24] minor fix --- src/multilevel.jl | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/src/multilevel.jl b/src/multilevel.jl index 6bf9f62..8e09c81 100644 --- a/src/multilevel.jl +++ b/src/multilevel.jl @@ -287,9 +287,12 @@ end abstract type AMGAlg end struct RugeStubenAMG <: AMGAlg end + +""" +""" struct SmoothedAggregationAMG <: AMGAlg - B::Union{<:AbstractMatrix,Nothing} - function SmoothedAggregationAMG(B::Union{AbstractMatrix,Nothing}=nothing) + B::Union{<:AbstractArray,Nothing} # `B` can be `Vector`, `Matrix`, or `nothing` + function SmoothedAggregationAMG(B::Union{AbstractArray,Nothing}=nothing) new(B) end end From 9ea61903e5140be87308ade64844df77476f3c54 Mon Sep 17 00:00:00 2001 From: Abdelrahman95 Date: Tue, 3 Jun 2025 11:09:29 +0200 Subject: [PATCH 21/24] add B to docstring --- src/aggregation.jl | 2 +- src/multilevel.jl | 32 ++++++++++++++++++++++++++++++++ test/nns_test.jl | 22 ++++++++++++++++++++++ 3 files changed, 55 insertions(+), 1 deletion(-) diff --git a/src/aggregation.jl b/src/aggregation.jl index 517a573..abc8bdb 100644 --- a/src/aggregation.jl +++ b/src/aggregation.jl @@ -14,7 +14,7 @@ function smoothed_aggregation(A::TA, _B = nothing, coarse_solver = Pinv, kwargs...) where {T,V,bs,TA<:SparseMatrixCSC{T,V}} n = size(A, 1) - B = isnothing(_B) ? ones(T,n,1) : copy(_B) + B = isnothing(_B) ? ones(T,n) : copy(_B) @assert size(A, 1) == size(B, 1) #=max_levels, max_coarse, strength = diff --git a/src/multilevel.jl b/src/multilevel.jl index 8e09c81..d9123bb 100644 --- a/src/multilevel.jl +++ b/src/multilevel.jl @@ -289,6 +289,38 @@ abstract type AMGAlg end struct RugeStubenAMG <: AMGAlg end """ + SmoothedAggregationAMG(B=nothing) + +Smoothed Aggregation AMG (Algebraic Multigrid) algorithm configuration. + +# Arguments +- `B::Union{AbstractArray, Nothing}`: The **near null space** in SA-AMG represents the low energy that cannot be attenuated by relaxtion, and thus needs to be perserved across the coarse grid. + + - `B` can be: + - a `Vector` (e.g., for scalar PDEs), + - a `Matrix` (e.g., for vector PDEs or systems with multiple equations), + - or `nothing`. + +# Notes +If `B` is set to `nothing`, it will be internally defaulted to `B = ones(T, n)`, where `T = eltype(A)` and `n = size(A, 1)`. + +# Examples + +**Poisson equation (scalar PDE):** +```julia +n = size(A, 1) +B = ones(Float64, n) +amg = SmoothedAggregationAMG(B) +``` + +**Linear elasticity equation in 2d (vector PDE):** +```julia +n = size(A, 1) # Ndof = 2 * number of nodes +B = zeros(Float64, n, 3) +B[1:2:end, :] = [1.0, 0.0, -y] +B[2:2:end, :] = [0.0, 1.0, x] +amg = SmoothedAggregationAMG(B) +``` """ struct SmoothedAggregationAMG <: AMGAlg B::Union{<:AbstractArray,Nothing} # `B` can be `Vector`, `Matrix`, or `nothing` diff --git a/test/nns_test.jl b/test/nns_test.jl index 678146e..abf9dc9 100644 --- a/test/nns_test.jl +++ b/test/nns_test.jl @@ -1,3 +1,25 @@ +## Test `B` as an argument for `SmoothedAggregationAMG` +@testset "Different `B` as an argument for `SmoothedAggregationAMG` " begin + A = poisson(100); + b = rand(100); + n = size(A,1) + T = eltype(A) + + #1. pass `B`` as nothing (Default case) + x_nothing = solve(A, b, SmoothedAggregationAMG(), maxiter = 1, abstol = 1e-6) + + #2. pass `B` as vector + B = ones(T,n) + x_vec = solve(A, b, SmoothedAggregationAMG(B), maxiter = 1, abstol = 1e-6) + @test x_vec ≈ x_nothing + + #3. pass `B` as matrix + B = ones(T,n,1) + x_mat = solve(A, b, SmoothedAggregationAMG(B), maxiter = 1, abstol = 1e-6) + @test x_mat ≈ x_nothing +end + + ## Test QR factorization @testset "fit_candidates unit test cases" begin cases = Vector{Tuple{SparseMatrixCSC{Float64,Int},Matrix{Float64}}}() From 700038195c680daac534dfb44eba6b52f686c4ca Mon Sep 17 00:00:00 2001 From: Abdelrahman95 Date: Tue, 10 Jun 2025 01:28:14 +0200 Subject: [PATCH 22/24] add B as kwargs --- src/aggregation.jl | 6 +++--- src/multilevel.jl | 44 ++------------------------------------------ src/precs.jl | 3 +-- test/nns_test.jl | 15 ++++++--------- test/runtests.jl | 1 - 5 files changed, 12 insertions(+), 57 deletions(-) diff --git a/src/aggregation.jl b/src/aggregation.jl index abc8bdb..c5e4e65 100644 --- a/src/aggregation.jl +++ b/src/aggregation.jl @@ -1,5 +1,6 @@ -function smoothed_aggregation(A::TA, _B = nothing, +function smoothed_aggregation(A::TA, ::Type{Val{bs}}=Val{1}; + B = nothing, symmetry = HermitianSymmetry(), strength = SymmetricStrength(), aggregate = StandardAggregation(), @@ -12,9 +13,8 @@ function smoothed_aggregation(A::TA, _B = nothing, diagonal_dominance = false, keep = false, coarse_solver = Pinv, kwargs...) where {T,V,bs,TA<:SparseMatrixCSC{T,V}} - n = size(A, 1) - B = isnothing(_B) ? ones(T,n) : copy(_B) + B = isnothing(B) ? ones(T,n) : copy(B) @assert size(A, 1) == size(B, 1) #=max_levels, max_coarse, strength = diff --git a/src/multilevel.jl b/src/multilevel.jl index d9123bb..242e0c0 100644 --- a/src/multilevel.jl +++ b/src/multilevel.jl @@ -287,47 +287,7 @@ end abstract type AMGAlg end struct RugeStubenAMG <: AMGAlg end - -""" - SmoothedAggregationAMG(B=nothing) - -Smoothed Aggregation AMG (Algebraic Multigrid) algorithm configuration. - -# Arguments -- `B::Union{AbstractArray, Nothing}`: The **near null space** in SA-AMG represents the low energy that cannot be attenuated by relaxtion, and thus needs to be perserved across the coarse grid. - - - `B` can be: - - a `Vector` (e.g., for scalar PDEs), - - a `Matrix` (e.g., for vector PDEs or systems with multiple equations), - - or `nothing`. - -# Notes -If `B` is set to `nothing`, it will be internally defaulted to `B = ones(T, n)`, where `T = eltype(A)` and `n = size(A, 1)`. - -# Examples - -**Poisson equation (scalar PDE):** -```julia -n = size(A, 1) -B = ones(Float64, n) -amg = SmoothedAggregationAMG(B) -``` - -**Linear elasticity equation in 2d (vector PDE):** -```julia -n = size(A, 1) # Ndof = 2 * number of nodes -B = zeros(Float64, n, 3) -B[1:2:end, :] = [1.0, 0.0, -y] -B[2:2:end, :] = [0.0, 1.0, x] -amg = SmoothedAggregationAMG(B) -``` -""" -struct SmoothedAggregationAMG <: AMGAlg - B::Union{<:AbstractArray,Nothing} # `B` can be `Vector`, `Matrix`, or `nothing` - function SmoothedAggregationAMG(B::Union{AbstractArray,Nothing}=nothing) - new(B) - end -end +struct SmoothedAggregationAMG <: AMGAlg end function solve(A::AbstractMatrix, b::Vector, s::AMGAlg, args...; kwargs...) solt = init(s, A, b, args...; kwargs...) @@ -337,7 +297,7 @@ function init(::RugeStubenAMG, A, b, args...; kwargs...) AMGSolver(ruge_stuben(A; kwargs...), b) end function init(sa::SmoothedAggregationAMG, A, b; kwargs...) - AMGSolver(smoothed_aggregation(A,sa.B; kwargs...), b) + AMGSolver(smoothed_aggregation(A; kwargs...), b) end function solve!(solt::AMGSolver, args...; kwargs...) _solve(solt.ml, solt.b, args...; kwargs...) diff --git a/src/precs.jl b/src/precs.jl index 3800ae8..59b683f 100644 --- a/src/precs.jl +++ b/src/precs.jl @@ -14,8 +14,7 @@ function SmoothedAggregationPreconBuilder(; blocksize = 1, kwargs...) end function (b::SmoothedAggregationPreconBuilder)(A::AbstractSparseMatrixCSC, p) - B = get(b.kwargs, :B, nothing) # extract nns from kwargs, default to `nothing` - return (aspreconditioner(smoothed_aggregation(SparseMatrixCSC(A), B, Val{b.blocksize}; b.kwargs...)),I) + return (aspreconditioner(smoothed_aggregation(SparseMatrixCSC(A), Val{b.blocksize}; b.kwargs...)), I) end diff --git a/test/nns_test.jl b/test/nns_test.jl index abf9dc9..3e86d2d 100644 --- a/test/nns_test.jl +++ b/test/nns_test.jl @@ -10,12 +10,12 @@ #2. pass `B` as vector B = ones(T,n) - x_vec = solve(A, b, SmoothedAggregationAMG(B), maxiter = 1, abstol = 1e-6) + x_vec = solve(A, b, SmoothedAggregationAMG(), maxiter = 1, abstol = 1e-6;B=B) @test x_vec ≈ x_nothing #3. pass `B` as matrix B = ones(T,n,1) - x_mat = solve(A, b, SmoothedAggregationAMG(B), maxiter = 1, abstol = 1e-6) + x_mat = solve(A, b, SmoothedAggregationAMG(), maxiter = 1, abstol = 1e-6;B=B) @test x_mat ≈ x_nothing end @@ -194,11 +194,8 @@ end @load "lin_elastic_2d.jld2" A b B A = SparseMatrixCSC(A.m, A.n, A.colptr, A.rowval, A.nzval) - x_nns, residuals_nns = solve(A, b, SmoothedAggregationAMG(B); log=true, reltol=1e-10) - x_wonns, residuals_wonns = solve(A, b, SmoothedAggregationAMG(); log=true, reltol=1e-10) - - ml = smoothed_aggregation(A, B) - @show ml + x_nns, residuals_nns = solve(A, b, SmoothedAggregationAMG(), log=true, reltol=1e-10;B=B) + x_wonns, residuals_wonns = solve(A, b, SmoothedAggregationAMG(), log=true, reltol=1e-10) println("No NNS: final residual at iteration ", length(residuals_wonns), ": ", residuals_wonns[end]) println("With NNS: final residual at iteration ", length(residuals_nns), ": ", residuals_nns[end]) @@ -233,8 +230,8 @@ end @test u[end-1] ≈ (P * L^3) / (3 * E * I) # vertical disp. at the end of the beam - x_nns, residuals_nns = solve(A, b, SmoothedAggregationAMG(B); log=true, reltol=1e-10) - x_wonns, residuals_wonns = solve(A, b, SmoothedAggregationAMG(); log=true, reltol=1e-10) + x_nns, residuals_nns = solve(A, b, SmoothedAggregationAMG(), log=true, reltol=1e-10,B=B) + x_wonns, residuals_wonns = solve(A, b, SmoothedAggregationAMG(), log=true, reltol=1e-10) println("No NNS: final residual at iteration ", length(residuals_wonns), ": ", residuals_wonns[end]) println("With NNS: final residual at iteration ", length(residuals_nns), ": ", residuals_nns[end]) diff --git a/test/runtests.jl b/test/runtests.jl index fa2cd0b..543084c 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -348,7 +348,6 @@ for sz in [ (10,10), (20,20), (50,50)] strategy = KrylovJL_CG(precs = RugeStubenPreconBuilder()) @test solve(prob, strategy, atol=1.0e-14) ≈ u0 rtol = 1.0e-8 - strategy = KrylovJL_CG(precs = SmoothedAggregationPreconBuilder()) @test solve(prob, strategy, atol=1.0e-14) ≈ u0 rtol = 1.0e-8 From 97976ce04472d920bf6ad21e21ab79101b00a778 Mon Sep 17 00:00:00 2001 From: Abdelrahman95 Date: Tue, 10 Jun 2025 01:29:26 +0200 Subject: [PATCH 23/24] fix format --- src/AlgebraicMultigrid.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/AlgebraicMultigrid.jl b/src/AlgebraicMultigrid.jl index c47fc4f..ea62549 100644 --- a/src/AlgebraicMultigrid.jl +++ b/src/AlgebraicMultigrid.jl @@ -9,7 +9,7 @@ using Printf @reexport import CommonSolve: solve, solve!, init using Reexport -using LinearAlgebra: rmul!,qr +using LinearAlgebra: rmul!, qr include("utils.jl") export approximate_spectral_radius From 25c959d1551bc168edd5d9de36dbd66d468d9d3d Mon Sep 17 00:00:00 2001 From: Abdelrahman95 Date: Tue, 10 Jun 2025 21:09:38 +0200 Subject: [PATCH 24/24] add docstring and make ruge stuben fail for B --- src/classical.jl | 4 ++++ src/multilevel.jl | 6 ++++++ 2 files changed, 10 insertions(+) diff --git a/src/classical.jl b/src/classical.jl index 8925471..3b4abc2 100644 --- a/src/classical.jl +++ b/src/classical.jl @@ -9,6 +9,10 @@ function ruge_stuben(_A::Union{TA, Symmetric{Ti, TA}, Hermitian{Ti, TA}}, max_coarse = 10, coarse_solver = Pinv, kwargs...) where {Ti,Tv,bs,TA<:SparseMatrixCSC{Ti,Tv}} + + # fails if near null space `B` is provided + haskey(kwargs, :B) && kwargs[:B] !== nothing && error("near null space `B` is only supported for smoothed aggregation AMG, not Ruge-Stüben AMG.") + if _A isa Symmetric && Ti <: Real || _A isa Hermitian A = _A.data symmetric = true diff --git a/src/multilevel.jl b/src/multilevel.jl index 242e0c0..18146e0 100644 --- a/src/multilevel.jl +++ b/src/multilevel.jl @@ -187,7 +187,13 @@ Keyword Arguments * maxiter::Int64 - maximum number of iterations to execute * verbose::Bool - display residual at each iteration * log::Bool - return vector of residuals along with solution +* B::AbstractArray - the **near null space** in SA-AMG, which represents the low energy that cannot be attenuated by relaxtion, and thus needs to be perserved across the coarse grid. +!!! note + `B` can be: + - a `Vector` (e.g., for scalar PDEs), + - a `Matrix` (e.g., for vector PDEs or systems with multiple equations), + If `B` is not provided, it defaults to a vector of ones. """ function _solve(ml::MultiLevel, b::AbstractArray, args...; kwargs...) n = length(ml) == 1 ? size(ml.final_A, 1) : size(ml.levels[1].A, 1)