rm Ferrite dependency

Author: Abdelrahman912

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"]

test/lin_elastic_2d.jld2

@@ -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,16 +168,18 @@ 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)
 
         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