diff --git a/src/ModelReduction.jl b/src/ModelReduction.jl index 9c4c014..74ec066 100644 --- a/src/ModelReduction.jl +++ b/src/ModelReduction.jl @@ -8,4 +8,6 @@ include("global_mass.jl") include("guyan_reduction.jl") include("craig_bampton.jl") include("sort_nodes.jl") +include("model_reduction_craig_bampton.jl") +export CraigBampton end diff --git a/src/model_reduction_craig_bampton.jl b/src/model_reduction_craig_bampton.jl new file mode 100644 index 0000000..ba49413 --- /dev/null +++ b/src/model_reduction_craig_bampton.jl @@ -0,0 +1,143 @@ +# This file is a part of JuliaFEM. +# License is MIT: see https://github.com/JuliaFEM/ModelReduction.jl/blob/master/LICENSE + +using FEMBase + +# TODO: abstract elimination of boundary conditions must be introduced +using JuliaFEM: eliminate_boundary_conditions!, + get_field_assembly, + get_boundary_problems + +type CraigBampton <: AbstractAnalysis + r_nodes :: Vector{Int64} + l_nodes :: Vector{Int64} + K :: Matrix{Float64} + M :: Matrix{Float64} + B :: Matrix{Float64} + dim :: Int64 + nev :: Int64 + time ::Float64 +end + +function CraigBampton() + return CraigBampton( + Vector{Int64}(0), + Vector{Int64}(0), + Matrix{Float64}(0,0), + Matrix{Float64}(0,0), + Matrix{Float64}(0,0), + 3, 10, 0.0) +end + +function FEMBase.run!(cb::Analysis{CraigBampton}) + time = cb.properties.time + dim = cb.properties.dim + r_nodes = cb.properties.r_nodes + l_nodes = cb.properties.l_nodes + nev = cb.properties.nev + r_dofs = [dim*(i-1)+j for i in r_nodes for j=1:dim] + info("number of retained nodes = ", length(r_nodes)) + info("number of retained dofs = ", length(r_dofs)) + l_dofs = [dim*(i-1)+j for i in l_nodes for j=1:dim] + info("number of other nodes = ", length(l_nodes)) + info("number of other dofs = ", length(l_dofs)) + + # Construct global matrices M and K (manually eliminate tie constraints) + assemble!(cb, time; with_mass_matrix=true) + M_red, K_red, Kg, f = get_field_assembly(cb) + dim1 = size(K_red, 1) + for boundary_problem in get_boundary_problems(cb) + eliminate_boundary_conditions!(K_red, M_red, boundary_problem, dim1) + end + SparseArrays.dropzeros!(K_red) + SparseArrays.dropzeros!(M_red) + + K_LL = K_red[l_dofs, l_dofs] + K_LL = 1/2*(K_LL + K_LL') + K_RR = K_red[r_dofs, r_dofs] + K_RL = K_red[r_dofs, l_dofs] + K_LR = K_red[l_dofs, r_dofs] + + M_LL = M_red[l_dofs, l_dofs] + M_LL = 1/2*(M_LL + M_LL') + M_RR = M_red[r_dofs, r_dofs] + M_RL = M_red[r_dofs, l_dofs] + M_LR = M_red[l_dofs, r_dofs] + + nz = get_nonzero_rows(K_LL) + nz2 = get_nonzero_columns(K_LR) + dim = size(K_LL, 1) + info("Dimension of K_LL: ", dim) + info("Number of nonzeros: ", length(nz)) + info("Zero rows should match to the 3*slave_nodes: ", dim-length(nz)) + + info("1. Static condensation") + + t0 = Base.time() + info("Construct recovery matrix B = K_LL^-1 * K_LR") + info("Factorize K_LL using Cholesky factorization") + CF = cholfact(K_LL[nz,nz]) + info("Calculate B = K_LL^-1 * K_LR blockwise") + dim1, dim2 = size(K_LR) + cb.properties.B = B = zeros(dim1, dim2) + for i=1:length(nz2) + if mod(i,10) == 0 + info("Condensate dof $i / $(length(nz2))...") + end + B[nz,nz2[i]] = CF \ full(K_LR[nz,nz2[i]]) + end + dt = round(Base.time() - t0, 2) + info("Calculation of recovery matrix B done in $dt seconds") + + t0 = Base.time() + info("Calculate boundary stiffness matrix K_BB") + K_BB = K_RR - K_RL*B + dt = round(Base.time() - t0, 2) + info("Calculating K_BB done in $dt seconds") + + t0 = Base.time() + info("Calculate boundary mass matrix M_BB") + info("Calculate Y = M_RL*B") + Y = -M_RL*B + Z = M_LL*B + info("Calculate M_BB = M_RR + Y + Y' + B'*M_LL*B") + M_BB = M_RR + Y + Y' + B'*Z + dt = round(Base.time() - t0, 2) + info("Calculating M_BB done in $dt seconds") + + info("2. Dynamic condensation") + + info("Calculate eigenvalue problem K_LL + w^2*M_LL = 0") + t0 = Base.time() + X_LL = zeros(size(K_LL, 1), nev) + w_LL, X_LL[nz,:] = eigs(K_LL[nz,nz], M_LL[nz,nz]; which=:SM, nev=nev) + w_LL = real(w_LL) + dt = round(Base.time() - t0, 2) + info("Eigenvalues solved in $dt seconds") + + info("Construct the final system M_red*u'' + K_red*u = f_red") + t0 = Base.time() + K_mm = X_LL' * K_LL * X_LL + M_mm = X_LL' * M_LL * X_LL + K_mB = zeros(nev, length(r_dofs)) + M_mB = X_LL' * (M_LR + Z) + K_Bm = K_mB' + M_Bm = M_mB' + cb.properties.M = M_tot = [M_BB M_Bm; M_mB M_mm] + cb.properties.K = K_tot = [K_BB K_Bm; K_mB K_mm] + dt = round(Base.time() - t0, 2) + info("Reduced system ready in $dt seconds") + + ndofs_reduced = size(K_red, 1) + nr_dofs = length(r_dofs) + nl_dofs = length(l_dofs) + nt_dofs = nr_dofs + nl_dofs + p = round(nr_dofs/nt_dofs*100, 2) + + info("--- CALCULATION SUMMARY ---") + info("Dimension of original system: $nt_dofs dofs") + info("Dimension of reduced system: boundary $nr_dofs dofs + $nev general dofs") + info("Dimension of reduced system is $p % of the original system") + + return M_tot, K_tot +end diff --git a/test/craig_bampton/calc_eigenvalues.jl b/test/craig_bampton/calc_eigenvalues.jl new file mode 100644 index 0000000..c4b803f --- /dev/null +++ b/test/craig_bampton/calc_eigenvalues.jl @@ -0,0 +1,64 @@ +using JuliaFEM +using JuliaFEM.Preprocess +using JuliaFEM.Postprocess +using JuliaFEM.Abaqus: create_surface_elements + +tic() + +# read mesh +mesh = abaqus_read_mesh("model.inp") +info("element sets = ", collect(keys(mesh.element_sets))) +info("surface sets = ", collect(keys(mesh.surface_sets))) + +# create two field problems with different material properties +bracket = Problem(Elasticity, "LDU_Bracket", 3) +bracket.elements = create_elements(mesh, "LDUBracket") +update!(bracket.elements, "youngs modulus", 165.0E3) +update!(bracket.elements, "poissons ratio", 0.275) +update!(bracket.elements, "density", 7.10E-9) +plate = Problem(Elasticity, "AdapterPlate", 3) +plate.elements = create_elements(mesh, "Adapterplate1", "Adapterplate2") +update!(plate.elements, "youngs modulus", 208.0E3) +update!(plate.elements, "poissons ratio", 0.30) +update!(plate.elements, "density", 7.80E-9) + +# create boundary condition from node set + +fixed = Problem(Dirichlet, "fixed", 3, "displacement") +fixed_nodes = mesh.node_sets[:FIXED] +fixed.elements = [Element(Poi1, [nid]) for nid in fixed_nodes] +update!(fixed.elements, "geometry", mesh.nodes) +update!(fixed.elements, "displacement 1", 0.0) +update!(fixed.elements, "displacement 2", 0.0) +update!(fixed.elements, "displacement 3", 0.0) + +""" A helper function to create tie contact. """ +function create_interface(mesh::Mesh, slave::String, master::String) + interface = Problem(Mortar, "tie contact", 3, "displacement") + interface.properties.dual_basis = false + slave_elements = create_surface_elements(mesh, slave) + master_elements = create_surface_elements(mesh, master) + update!(slave_elements, "master elements", master_elements) + interface.elements = [slave_elements; master_elements] + return interface +end + +# call helper function to create tie contacts +tie1 = create_interface(mesh, + "LDUBracketToAdapterplate1", + "Adapterplate1ToLDUBracket") +tie2 = create_interface(mesh, + "LDUBracketToAdapterplate2", + "Adapterplate2ToLDUBracket") + +# add field and boundary problems to solver +solver = Solver(Modal, bracket, plate, fixed, tie1, tie2) +# save results to Xdmf data format ready for ParaView visualization +add_results_writer!(solver, Xdmf("results"; overwrite=true)) +# solve 6 smallest eigenvalues +solver.properties.nev = 10 +solver.properties.which = :SM +solver() +println("Eigenvalues: ", sqrt.(solver.properties.eigvals) / (2*pi)) + +toc() diff --git a/test/test_model_reduction_craig_bampton.jl b/test/test_model_reduction_craig_bampton.jl new file mode 100644 index 0000000..d8318e0 --- /dev/null +++ b/test/test_model_reduction_craig_bampton.jl @@ -0,0 +1,93 @@ +# This file is a part of JuliaFEM. +# License is MIT: see https://github.com/JuliaFEM/ModelReduction.jl/blob/master/LICENSE + +using ModelReduction + +using JuliaFEM +using JuliaFEM.Preprocess +using JuliaFEM.Postprocess +using JuliaFEM.Abaqus: create_surface_elements +using Logging +Logging.configure(level=DEBUG) + +# read mesh +datadir = joinpath(Pkg.dir("ModelReduction"), "test", "test_model_reduction_craig_bampton") +mesh = abaqus_read_mesh(joinpath(datadir, "model.inp")) +info("element sets = ", collect(keys(mesh.element_sets))) +info("surface sets = ", collect(keys(mesh.surface_sets))) +info("node sets = ", collect(keys(mesh.node_sets))) + +#' create two field problems with different material properties + +bracket = Problem(Elasticity, "LDU_Bracket", 3) +bracket.elements = create_elements(mesh, "LDUBracket") +update!(bracket.elements, "youngs modulus", 165.0E3) +update!(bracket.elements, "poissons ratio", 0.275) +update!(bracket.elements, "density", 7.10E-9) +plate = Problem(Elasticity, "AdapterPlate", 3) +plate.elements = create_elements(mesh, "Adapterplate1", "Adapterplate2") +update!(plate.elements, "youngs modulus", 208.0E3) +update!(plate.elements, "poissons ratio", 0.30) +update!(plate.elements, "density", 7.80E-9) + +# create boundary condition from node set + +fixed = Problem(Dirichlet, "fixed", 3, "displacement") +fixed_nodes = mesh.node_sets[:FIXED] +fixed.elements = [Element(Poi1, [nid]) for nid in fixed_nodes] +update!(fixed.elements, "geometry", mesh.nodes) +update!(fixed.elements, "displacement 1", 0.0) +update!(fixed.elements, "displacement 2", 0.0) +update!(fixed.elements, "displacement 3", 0.0) + +""" A helper function to create tie contact. """ +function create_interface(mesh::Mesh, slave::String, master::String) + interface = Problem(Mortar, "tie contact", 3, "displacement") + interface.properties.dual_basis = false + slave_elements = create_surface_elements(mesh, slave) + master_elements = create_surface_elements(mesh, master) + update!(slave_elements, "master elements", master_elements) + interface.elements = [slave_elements; master_elements] + return interface +end + +# call helper function to create tie contacts +tie1 = create_interface(mesh, + "LDUBracketToAdapterplate1", + "Adapterplate1ToLDUBracket") +tie2 = create_interface(mesh, + "LDUBracketToAdapterplate2", + "Adapterplate2ToLDUBracket") + +# solver = Solver(Modal, bracket, plate, fixed, tie1, tie2) +# nev = solver.properties.nev = 10 +# solver.properties.which = :SM +# solver.properties.sigma = 1.0e-9 +# solver.properties.empty_assemblies_before_solution = false +# solver() +# w = sqrt.(real(solver.properties.eigvals)) +# info("Eigenvalues of the original system [Hz]: ", w / (2*pi)) + +# Let's test model reduction + +tic() + +cb = Analysis(CraigBampton) +# Retained nodes = nodes from node set BORDER +cb.properties.r_nodes = r_nodes = collect(mesh.node_sets[:BORDER]) +cb.properties.l_nodes = setdiff(keys(mesh.nodes), r_nodes) +add_problems!(cb, [bracket, plate, fixed, tie1, tie2]) +@time run!(cb) + +info("Calculate 10 smallest eigenvalues of the reduced system") +t0 = Base.time() +w_, X_ = eigs(cb.properties.K, cb.properties.M; which=:SM, nev=10) +w_ = sqrt.(real(w_)) +dt = round(Base.time() - t0, 2) +info("Calculated eigenvalues in $dt seconds") + +w = [111.383, 155.03, 215.399, 358.761, 409.654, 603.515, 714.478, 859.087, 906.233, 1134.16] +info("Eigenvalues of the original system [Hz]: ", w) +info("Eigenvalues of the reduced system [Hz]: ", w_/(2*pi)) + +toc() diff --git a/test/test_model_reduction_craig_bampton.log b/test/test_model_reduction_craig_bampton.log new file mode 100644 index 0000000..fd6521c --- /dev/null +++ b/test/test_model_reduction_craig_bampton.log @@ -0,0 +1,61 @@ +11-May 22:09:04:INFO:root:number of retained nodes = 64 +11-May 22:09:04:INFO:root:number of retained dofs = 192 +11-May 22:09:04:INFO:root:number of other nodes = 97706 +11-May 22:09:04:INFO:root:number of other dofs = 293118 +11-May 22:09:05:INFO:root:Assembling problems ... +11-May 22:11:09:INFO:root:Assembly done! +11-May 22:11:21:INFO:root:Eliminating mesh tie constraint tie contact using static condensation +11-May 22:11:21:INFO:root:# slave dofs = 1344, # master dofs = 105 +11-May 22:11:21:WARNING:root:D is not diagonal, is dual basis used? This might take a long time. +11-May 22:11:29:INFO:root:Eliminating mesh tie constraint tie contact using static condensation +11-May 22:11:29:INFO:root:# slave dofs = 2544, # master dofs = 174 +11-May 22:11:29:WARNING:root:D is not diagonal, is dual basis used? This might take a long time. +11-May 22:11:43:INFO:root:Dimension of K_LL: 293118 +11-May 22:11:43:INFO:root:Number of nonzeros: 286167 +11-May 22:11:43:INFO:root:Zero rows should match to the 3*slave_nodes: 6951 +11-May 22:11:43:INFO:root:1. Static condensation +11-May 22:11:43:INFO:root:Construct recovery matrix B = K_LL^-1 * K_LR +11-May 22:11:43:INFO:root:Factorize K_LL using Cholesky factorization +11-May 22:11:50:INFO:root:Calculate B = K_LL^-1 * K_LR blockwise +11-May 22:11:53:INFO:root:Condensate dof 10 / 192... +11-May 22:11:55:INFO:root:Condensate dof 20 / 192... +11-May 22:11:57:INFO:root:Condensate dof 30 / 192... +11-May 22:11:59:INFO:root:Condensate dof 40 / 192... +11-May 22:12:01:INFO:root:Condensate dof 50 / 192... +11-May 22:12:03:INFO:root:Condensate dof 60 / 192... +11-May 22:12:05:INFO:root:Condensate dof 70 / 192... +11-May 22:12:07:INFO:root:Condensate dof 80 / 192... +11-May 22:12:09:INFO:root:Condensate dof 90 / 192... +11-May 22:12:11:INFO:root:Condensate dof 100 / 192... +11-May 22:12:13:INFO:root:Condensate dof 110 / 192... +11-May 22:12:16:INFO:root:Condensate dof 120 / 192... +11-May 22:12:18:INFO:root:Condensate dof 130 / 192... +11-May 22:12:20:INFO:root:Condensate dof 140 / 192... +11-May 22:12:22:INFO:root:Condensate dof 150 / 192... +11-May 22:12:24:INFO:root:Condensate dof 160 / 192... +11-May 22:12:26:INFO:root:Condensate dof 170 / 192... +11-May 22:12:28:INFO:root:Condensate dof 180 / 192... +11-May 22:12:30:INFO:root:Condensate dof 190 / 192... +11-May 22:12:31:INFO:root:Calculation of recovery matrix B done in 48.0 seconds +11-May 22:12:31:INFO:root:Calculate boundary stiffness matrix K_BB +11-May 22:12:31:INFO:root:Calculating K_BB done in 0.1 seconds +11-May 22:12:31:INFO:root:Calculate boundary mass matrix M_BB +11-May 22:12:31:INFO:root:Calculate Y = M_RL*B +11-May 22:12:36:INFO:root:Calculate M_BB = M_RR + Y + Y' + B'*M_LL*B +11-May 22:12:36:INFO:root:Calculating M_BB done in 5.42 seconds +11-May 22:12:36:INFO:root:2. Dynamic condensation +11-May 22:12:36:INFO:root:Calculate eigenvalue problem K_LL + w^2*M_LL = 0 +11-May 22:12:59:INFO:root:Eigenvalues solved in 23.07 seconds +11-May 22:12:59:INFO:root:Construct the final system M_red*u'' + K_red*u = f_red +11-May 22:13:02:INFO:root:Reduced system ready in 2.8 seconds +11-May 22:13:02:INFO:root:--- CALCULATION SUMMARY --- +11-May 22:13:02:INFO:root:Dimension of original system: 293310 dofs +11-May 22:13:02:INFO:root:Dimension of reduced system: boundary 192 dofs + 10 general dofs +11-May 22:13:02:INFO:root:Dimension of reduced system is 0.07 % of the original system +243.121280 seconds (1.04 G allocations: 74.489 GiB, 18.68% gc time) +11-May 22:13:02:INFO:root:Calculate 10 smallest eigenvalues of the reduced system +11-May 22:13:03:INFO:root:Calculated eigenvalues in 1.16 seconds +11-May 22:13:04:INFO:root:Eigenvalues of the original system [Hz]: [111.383, 155.03, 215.399, 358.761, 409.654, 603.515, 714.48, 859.087, 906.23, 1134.16] +11-May 22:13:04:INFO:root:Eigenvalues of the reduced system [Hz]: [111.377, 155.029, 215.48, 358.793, 410.932, 604.705, 720.236, 861.946, 923.735, 1147.1] + +# difference in pct: [0.0056405, 0.000682211, -0.0376501, -0.00889601, -0.311908, -0.197173, -0.805642, -0.332753, -1.93162, -1.14103]