-
Notifications
You must be signed in to change notification settings - Fork 8
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
Implementation of the Craig-Bampton method
Implemented. Test model is not included to the repository.
- Loading branch information
Showing
5 changed files
with
363 additions
and
0 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -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 |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -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() |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -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() |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -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] |