From 4e6ecf19d63e08e722450fd51cc39b4711a4f4ff Mon Sep 17 00:00:00 2001 From: Marcelo Forets Date: Sat, 17 Jun 2023 15:51:10 -0300 Subject: [PATCH 1/2] use single field in FixedDof --- ...bench_linear_cylinder_internal_pressure.jl | 12 ++-- .../bench_linear_extension.jl | 12 ++-- .../bench_uniaxial_compression.jl | 12 ++-- .../bench_uniaxial_extension.jl | 8 +-- examples/clamped_truss/clamped_truss.jl | 2 +- .../cylinder_internal_pressure.jl | 26 ++++----- examples/linear_extension/linear_extension.jl | 14 ++--- .../uniaxial_compression.jl | 56 +++++++++---------- .../uniaxial_extension/uniaxial_extension.jl | 6 +- examples/von_misses_truss/von_misses_truss.jl | 20 +++---- .../FixedDofBoundaryConditions.jl | 37 ++++++------ src/Entities/Nodes.jl | 2 +- test/boundary_conditions.jl | 41 ++++++++------ test/static_analyses.jl | 5 +- test/structural_model.jl | 6 +- 15 files changed, 133 insertions(+), 126 deletions(-) diff --git a/benchmark/linear_cylinder_internal_pressure/bench_linear_cylinder_internal_pressure.jl b/benchmark/linear_cylinder_internal_pressure/bench_linear_cylinder_internal_pressure.jl index 19f01ab54..b50ef379c 100644 --- a/benchmark/linear_cylinder_internal_pressure/bench_linear_cylinder_internal_pressure.jl +++ b/benchmark/linear_cylinder_internal_pressure/bench_linear_cylinder_internal_pressure.jl @@ -66,11 +66,11 @@ function linear_cylinder_structure(; ms::Real=0.5) # ------------------------------- # Boundary conditions # ------------------------------- - # Dirichlet boundary conditions - bc₁ = FixedDof(; components=[1], name=bc₁_label) - bc₂ = FixedDof(; components=[2], name=bc₂_label) - bc₃ = FixedDof(; components=[3], name=bc₃_label) - # Neumann boundary conditions + # Dirichlet boundary conditions + bc₁ = FixedDof(:u, [1], bc₁_label) + bc₂ = FixedDof(:u, [2], bc₂_label) + bc₃ = FixedDof(:u, [3], bc₃_label) + # Neumann boundary conditions bc₄ = Pressure(; values=pressure, name=bc₄_label) boundary_conditions = StructuralBoundaryCondition(bc₁, bc₂, bc₃, bc₄) # Assign boundary conditions to the ones defined in the mesh @@ -98,7 +98,7 @@ function point_eval_handler(structure::Structure; ## scalar parameters (dimensions in mm an MPa) Lₖ = 30.0 # cylinder length in 𝐞ₖ mm Rₑ = 200.0 # outer radius in mm - Lᵢ = Lⱼ = 2.25Rₑ # hyper rectangle origin in 𝐞ᵢ,𝐞ⱼ and 𝐞ₖ in mm + Lᵢ = Lⱼ = 2.25Rₑ # hyper rectangle origin in 𝐞ᵢ,𝐞ⱼ and 𝐞ₖ in mm O = (x=-Lᵢ / 2, y=-Lⱼ / 2, z=0.0) # hyper rectangle origin in 𝐞ᵢ,𝐞ⱼ and 𝐞ₖ in mm # Create an hyper rectangle Lᵢ x Lⱼ x Lₖ diff --git a/benchmark/linear_extension/bench_linear_extension.jl b/benchmark/linear_extension/bench_linear_extension.jl index 85c75fbf6..27df82ba7 100644 --- a/benchmark/linear_extension/bench_linear_extension.jl +++ b/benchmark/linear_extension/bench_linear_extension.jl @@ -9,7 +9,7 @@ Linear extension GMSH mesh and `IsotropicLinearElastic` material. function linear_extension_structure(; ms=0.5) # x, y and z dimensions of the box in the mesh respectively. - Lᵢ = 2.0 # Dimension in x of the box in m + Lᵢ = 2.0 # Dimension in x of the box in m Lⱼ = 1.0 # Dimension in y of the box in m Lₖ = 1.0 # Dimension in z of the box in m @@ -25,16 +25,16 @@ function linear_extension_structure(; ms=0.5) # Material types without assigned elements. materials = StructuralMaterial(mat) - # Dirichlet boundary conditions + # Dirichlet boundary conditions bc₁_label = "fixed-ux" bc₂_label = "fixed-uj" bc₃_label = "fixed-uk" bc₄_label = "tension" - bc₁ = FixedDof([:u], [1], bc₁_label) - bc₂ = FixedDof([:u], [2], bc₂_label) - bc₃ = FixedDof([:u], [3], bc₃_label) + bc₁ = FixedDof(:u, [1], bc₁_label) + bc₂ = FixedDof(:u, [2], bc₂_label) + bc₃ = FixedDof(:u, [3], bc₃_label) - # Neumann boundary conditions + # Neumann boundary conditions bc₄ = GlobalLoad([:u], t -> [p * t, 0, 0], bc₄_label) boundary_conditions = StructuralBoundaryCondition(bc₁, bc₂, bc₃, bc₄) diff --git a/benchmark/uniaxial_compression/bench_uniaxial_compression.jl b/benchmark/uniaxial_compression/bench_uniaxial_compression.jl index 67c0130da..521afb3b3 100644 --- a/benchmark/uniaxial_compression/bench_uniaxial_compression.jl +++ b/benchmark/uniaxial_compression/bench_uniaxial_compression.jl @@ -7,7 +7,7 @@ function strain_energy_neo(𝔼::AbstractMatrix, K::Real, μ::Real) J = sqrt(det(ℂ)) # First invariant I₁ = tr(ℂ) - # Strain energy function + # Strain energy function return Ψ = μ / 2 * (I₁ - 2 * log(J)) + K / 2 * (J - 1)^2 end @@ -39,16 +39,16 @@ function uniaxial_compression_structure(; ms=0.5) # Material types without assigned elements. materials = StructuralMaterial(neo_hookean_hyper) - # Dirichlet boundary conditions + # Dirichlet boundary conditions bc₁_label = "fixed-ux" bc₂_label = "fixed-uj" bc₃_label = "fixed-uk" bc₄_label = "tension" - bc₁ = FixedDof([:u], [1], bc₁_label) - bc₂ = FixedDof([:u], [2], bc₂_label) - bc₃ = FixedDof([:u], [3], bc₃_label) + bc₁ = FixedDof(:u, [1], bc₁_label) + bc₂ = FixedDof(:u, [2], bc₂_label) + bc₃ = FixedDof(:u, [3], bc₃_label) - # Neumann boundary conditions + # Neumann boundary conditions bc₄ = Pressure([:u], t -> p * t, bc₄_label) bc_labels = [bc₁_label, bc₂_label, bc₃_label, bc₄_label] diff --git a/benchmark/uniaxial_extension/bench_uniaxial_extension.jl b/benchmark/uniaxial_extension/bench_uniaxial_extension.jl index ba2429ae4..362f524ce 100644 --- a/benchmark/uniaxial_extension/bench_uniaxial_extension.jl +++ b/benchmark/uniaxial_extension/bench_uniaxial_extension.jl @@ -28,14 +28,14 @@ function uniaxial_extension_structure(; ms=0.5) # Material types without assigned elements. materials = StructuralMaterial(svk) - # Dirichlet boundary conditions + # Dirichlet boundary conditions bc₁_label = "fixed-ux" bc₂_label = "fixed-uj" bc₃_label = "fixed-uk" bc₄_label = "tension" - bc₁ = FixedDof([:u], [1], bc₁_label) - bc₂ = FixedDof([:u], [2], bc₂_label) - bc₃ = FixedDof([:u], [3], bc₃_label) + bc₁ = FixedDof(:u, [1], bc₁_label) + bc₂ = FixedDof(:u, [2], bc₂_label) + bc₃ = FixedDof(:u, [3], bc₃_label) # Neumann boundary conditions the tension is apliad in the negative direction of the z local axis. bc₄ = Pressure(; values=t -> -p * t, name=bc₄_label) diff --git a/examples/clamped_truss/clamped_truss.jl b/examples/clamped_truss/clamped_truss.jl index 927770ed2..554e637a4 100644 --- a/examples/clamped_truss/clamped_truss.jl +++ b/examples/clamped_truss/clamped_truss.jl @@ -41,7 +41,7 @@ function run_clamped_truss_example() # Boundary conditions # ------------------------------- # Fixed dofs - bc₁ = FixedDof(; components=[1], name="fixed_uₓ") + bc₁ = FixedDof(:u, [1], "fixed_uₓ") # Load bc₂ = GlobalLoad(; values=t -> [F * t], name="load in j") # Apply bcs to the nodes diff --git a/examples/cylinder_internal_pressure/cylinder_internal_pressure.jl b/examples/cylinder_internal_pressure/cylinder_internal_pressure.jl index 09e7e4af9..67e0d4c97 100644 --- a/examples/cylinder_internal_pressure/cylinder_internal_pressure.jl +++ b/examples/cylinder_internal_pressure/cylinder_internal_pressure.jl @@ -1,5 +1,5 @@ # -------------------------------------------------- -# Cylinder submitted to an Internal Pressure Example +# Cylinder submitted to an Internal Pressure Example #---------------------------------------------------- using LinearAlgebra, Test, Suppressor using ONSAS @@ -17,7 +17,7 @@ function run_cylinder_internal_pressure_example() E = 210.0 # Young modulus in MPa ν = 0.3 # Poisson ratio pressure(t::Real) = p * t - ## number of steps + ## number of steps NSTEPS = 9 ## tolerances for testing ATOL = 1e-2 * (Rₑ - Rᵢ) @@ -76,11 +76,11 @@ function run_cylinder_internal_pressure_example() # ------------------------------- # Boundary conditions # ------------------------------- - # Dirichlet boundary conditions - bc₁ = FixedDof(; components=[1], name=bc₁_label) - bc₂ = FixedDof(; components=[2], name=bc₂_label) - bc₃ = FixedDof(; components=[3], name=bc₃_label) - # Neumann boundary conditions + # Dirichlet boundary conditions + bc₁ = FixedDof(:u, [1], bc₁_label) + bc₂ = FixedDof(:u, [2], bc₂_label) + bc₃ = FixedDof(:u, [3], bc₃_label) + # Neumann boundary conditions bc₄ = Pressure(; values=pressure, name=bc₄_label) boundary_conditions = StructuralBoundaryCondition(bc₁, bc₂, bc₃, bc₄) # Assign boundary conditions to the ones defined in the mesh @@ -118,9 +118,9 @@ function run_cylinder_internal_pressure_example() # Numerical solution # ------------------------------- states_lin_sol = solve!(linear_analysis) - # get time vector or load factors + # get time vector or load factors λᵥ = load_factors(linear_analysis) - # Get the solution at a random point + # Get the solution at a random point "Return a rand point in the cylinder (R, θ, L)." function rand_point_cylinder(Rᵢ::Real=Rᵢ, Rₑ::Real=Rₑ, Lₖ::Real=Lₖ) [rand() * (Rₑ - Rᵢ) + Rᵢ, rand() * 2 * π, rand() * Lₖ] @@ -142,7 +142,7 @@ function run_cylinder_internal_pressure_example() uᵢ_numeric_p_rand = displacements(states_lin_sol, point_evaluator, 1) uⱼ_numeric_p_rand = displacements(states_lin_sol, point_evaluator, 2) uₖ_numeric_p_rand = displacements(states_lin_sol, point_evaluator, 3) - # + # uᵣ_numeric_p_rand = sqrt.(@. uᵢ_numeric_p_rand^2 + uⱼ_numeric_p_rand^2) # ------------------------------- # Analytic solution @@ -173,7 +173,7 @@ function run_cylinder_internal_pressure_example() atol::Real=ATOL, atolr=ATOLR, Rᵢ::Real=Rᵢ, Rₑ::Real=Rₑ, Lₖ::Real=Lₖ) structure = ONSAS.structure(analysis(sol)) - # Generic surface s at z = Lₖ + # Generic surface s at z = Lₖ rand_R, rand_θ₁, Lₖ = rand_point_cylinder(Rᵢ, Rₑ, Lₖ) # Set by force Lₖ rand_θ₂ = rand() * 2 * π @@ -250,7 +250,7 @@ function run_cylinder_internal_pressure_example() @test zero_uⱼ_axis_x_nonlinear end #----------------------------- - # Plot & plots + # Plot & plots #----------------------------- PLOT_RESULTS && plot_results(λᵥ, nonlinear_analysis, uᵣ_numeric_nᵢ, uᵣ_numeric_nₑ, uᵣ_analytic_nᵢ, uᵣ_analytic_nₑ) @@ -272,7 +272,7 @@ function plot_results(λᵥ, nonlinear_analysis, uᵣ_numeric_nᵢ, uᵣ_numeric plot!(fig, vec_p, uᵣ_analytic_nₑ; label="analytic linear uᵣ(Rₑ)", legend=:topleft, color=:black, lw=2, ls=:solid) - # Plot comparing linear and non linear solutions + # Plot comparing linear and non linear solutions plot!(fig, vec_p_non_in, uᵣ_numeric_nonlinear_nᵢ; label="non-linear uᵣ(0, Rᵢ, 0)", color=:red, lw=2, marker=:circle, markersize=3) diff --git a/examples/linear_extension/linear_extension.jl b/examples/linear_extension/linear_extension.jl index 1a6c96d6e..e1e28a033 100644 --- a/examples/linear_extension/linear_extension.jl +++ b/examples/linear_extension/linear_extension.jl @@ -1,4 +1,4 @@ -# ---------------------------------------------------------------- +# ---------------------------------------------------------------- # Uniaxial Extension Example 1 from (Zerpa et. Al., 2019, CMAME). # ---------------------------------------------------------------- using Test, LinearAlgebra, Suppressor @@ -17,7 +17,7 @@ function run_linear_extension_example() ν = 0.4 # Poisson's ratio p = 3 # Tension load in Pa tension(t) = p * t # Tension load function - Lᵢ = 2.0 # Dimension in x of the box in m + Lᵢ = 2.0 # Dimension in x of the box in m Lⱼ = 1.0 # Dimension in y of the box in m Lₖ = 1.0 # Dimension in z of the box in m RTOL = 1e-4 # Relative tolerance for tests @@ -34,11 +34,11 @@ function run_linear_extension_example() # ------------------------------- # Fixed dofs bc₁_label = "fixed-ux" - bc₁ = FixedDof([:u], [1], bc₁_label) + bc₁ = FixedDof(:u, [1], bc₁_label) bc₂_label = "fixed-uj" - bc₂ = FixedDof([:u], [2], bc₂_label) + bc₂ = FixedDof(:u, [2], bc₂_label) bc₃_label = "fixed-uk" - bc₃ = FixedDof([:u], [3], bc₃_label) + bc₃ = FixedDof(:u, [3], bc₃_label) # Load bc₄_label = "tension" bc₄ = GlobalLoad([:u], t -> [tension(t), 0, 0], bc₄_label) @@ -172,7 +172,7 @@ function run_linear_extension_example() σⱼ_analytic_p_rand_e = σ_analytic_p_rand_e[2] σₖ_analytic_p_rand_e = σ_analytic_p_rand_e[3] #----------------------------- - # Test boolean for CI + # Test boolean for CI #----------------------------- @testset "Linear Extension example" begin # Displacements @@ -186,7 +186,7 @@ function run_linear_extension_example() @test ϵᵢ_numeric_e_rand ≈ ϵᵢ_analytic_p_rand_e rtol = RTOL @test norm(ϵⱼ_numeric_e_rand) ≈ 0 atol = RTOL @test norm(ϵₖ_numeric_e_rand) ≈ 0 atol = RTOL - # Stresses + # Stresses @test σᵢ_analytic_p_rand_e ≈ σᵢ_analytic_p_rand_e rtol = RTOL @test σⱼ_analytic_p_rand_e ≈ σⱼ_analytic_p_rand_e atol = RTOL @test σₖ_analytic_p_rand_e ≈ σₖ_analytic_p_rand_e atol = RTOL diff --git a/examples/uniaxial_compression/uniaxial_compression.jl b/examples/uniaxial_compression/uniaxial_compression.jl index 0f1cafa38..1b41bee74 100644 --- a/examples/uniaxial_compression/uniaxial_compression.jl +++ b/examples/uniaxial_compression/uniaxial_compression.jl @@ -1,4 +1,4 @@ -# ----------------------------- +# ----------------------------- # Uniaxial Compression Example # ----------------------------- using Test, LinearAlgebra, Suppressor @@ -12,10 +12,10 @@ function run_uniaxial_compression() ## scalar parameters E = 1.0 # Young modulus in Pa ν = 0.3 # Poisson's ratio - μ = G = E / (2 * (1 + ν)) # Second Lamé parameter + μ = G = E / (2 * (1 + ν)) # Second Lamé parameter K = E / (3 * (1 - 2 * ν)) # Bulk modulus p = 1 # Pressure load in Pa - Lᵢ = 2.0 # Dimension in x of the box in m + Lᵢ = 2.0 # Dimension in x of the box in m Lⱼ = 1.0 # Dimension in y of the box in m Lₖ = 1.0 # Dimension in z of the box in m ms = 0.5 # Refinement factor for the mesh @@ -39,7 +39,7 @@ function run_uniaxial_compression() n₈ = Node(Lᵢ, Lⱼ, 0.0) vec_nodes = [n₁, n₂, n₃, n₄, n₅, n₆, n₇, n₈] s₁_mesh = Mesh(; nodes=vec_nodes) - ## Faces + ## Faces f₁ = TriangularFace(n₅, n₈, n₆, "loaded_face_1") f₂ = TriangularFace(n₆, n₈, n₇, "loaded_face_2") f₃ = TriangularFace(n₄, n₁, n₂, "x=0_face_1") @@ -50,7 +50,7 @@ function run_uniaxial_compression() f₈ = TriangularFace(n₄, n₈, n₅, "z=0_face_2") vec_faces = [f₁, f₂, f₃, f₄, f₅, f₆, f₇, f₈] append!(faces(s₁_mesh), vec_faces) - ## Entities + ## Entities t₁ = Tetrahedron(n₁, n₄, n₂, n₆, "tetra_1") t₂ = Tetrahedron(n₆, n₂, n₃, n₄, "tetra_2") t₃ = Tetrahedron(n₄, n₃, n₆, n₇, "tetra_3") @@ -76,15 +76,15 @@ function run_uniaxial_compression() # ------------------------------- # Fixed dofs bc₁_label = "fixed-ux" - bc₁ = FixedDof(; components=[1], name=bc₁_label) + bc₁ = FixedDof(:u, [1], bc₁_label) bc₂_label = "fixed-uj" - bc₂ = FixedDof(; components=[2], name=bc₂_label) + bc₂ = FixedDof(:u, [2], bc₂_label) bc₃_label = "fixed-uk" - bc₃ = FixedDof(; components=[3], name=bc₃_label) + bc₃ = FixedDof(:u, [3], bc₃_label) # Load bc₄_label = "compression" bc₄ = GlobalLoad(; values=t -> [-p * t, 0, 0], name=bc₄_label) - # Assign this to faces + # Assign this to faces face_bc = dictionary([bc₁ => [f₃, f₄], bc₂ => [f₅, f₆], bc₃ => [f₇, f₈], bc₄ => [f₁, f₂]]) # Crete boundary conditions struct s₁_boundary_conditions = StructuralBoundaryCondition(; face_bcs=face_bc) @@ -136,17 +136,17 @@ function run_uniaxial_compression() e = rand(elements(s₁)) # Cosserat or second Piola-Kirchhoff stress tensor ℙ_numeric_case₁ = stress(states_sol_case₁, e) - # ℙᵢᵢ component: + # ℙᵢᵢ component: ℙᵢᵢ_numeric_case₁ = getindex.(ℙ_numeric_case₁, 1, 1) - # ℙⱼⱼ component: + # ℙⱼⱼ component: ℙⱼⱼ_numeric_case₁ = getindex.(ℙ_numeric_case₁, 2, 2) - # ℙₖₖ component: + # ℙₖₖ component: ℙₖₖ_numeric_case₁ = getindex.(ℙ_numeric_case₁, 3, 3) - # Get the Right hand Cauchy strain tensor ℂ at a random state + # Get the Right hand Cauchy strain tensor ℂ at a random state ℂ_rand_numeric_case₁ = rand(strain(states_sol_case₁, e)) - # Get the Second Piola Kirchhoff stress tensor ℙ at a random state + # Get the Second Piola Kirchhoff stress tensor ℙ at a random state ℙ_rand_numeric_case₁ = rand(stress(states_sol_case₁, e)) - # Load factors + # Load factors load_factors_case₁ = load_factors(sa₁) # ----------------------------------------------- # Case 2 - GMSH mesh and `HyperElastic` material @@ -163,7 +163,7 @@ function run_uniaxial_compression() J = sqrt(det(ℂ)) # First invariant I₁ = tr(ℂ) - # Strain energy function + # Strain energy function Ψ = μ / 2 * (I₁ - 2 * log(J)) + K / 2 * (J - 1)^2 end params = [K, μ] # The order must be the same defined in the strain energy (splatting) @@ -174,7 +174,7 @@ function run_uniaxial_compression() # ------------------------------- # Boundary Conditions # ------------------------------- - # Redefine the load boundary condition + # Redefine the load boundary condition bc₄ = Pressure(; values=t -> p * t, name=bc₄_label) # BoundaryConditions types without assigned node, feces and elements s_boundary_conditions = StructuralBoundaryCondition(bc₁, bc₂, bc₃, bc₄) @@ -226,20 +226,20 @@ function run_uniaxial_compression() numeric_α_case₂, numeric_β_case₂, numeric_γ_case₂, numeric_uᵢ_case₂, _, _ = αβγ_numeric(states_sol_case₂) # Cosserat or second Piola-Kirchhoff stress tensor ℙ_numeric_case₂ = stress(states_sol_case₂, e) - # ℙᵢᵢ component: + # ℙᵢᵢ component: ℙᵢᵢ_numeric_case₂ = getindex.(ℙ_numeric_case₂, 1, 1) - # ℙⱼⱼ component: + # ℙⱼⱼ component: ℙⱼⱼ_numeric_case₂ = getindex.(ℙ_numeric_case₂, 2, 2) - # ℙₖₖ component: + # ℙₖₖ component: ℙₖₖ_numeric_case₂ = getindex.(ℙ_numeric_case₂, 3, 3) - # Get the Right hand Cauchy strain tensor ℂ at a random state + # Get the Right hand Cauchy strain tensor ℂ at a random state ℂ_rand_numeric_case₂ = rand(strain(states_sol_case₂, e)) - # Get the Second Piola Kirchhoff stress tensor ℙ at a random state + # Get the Second Piola Kirchhoff stress tensor ℙ at a random state ℙ_rand_numeric_case₂ = rand(stress(states_sol_case₂, e)) - # Load factors + # Load factors load_factors_case₂ = load_factors(sa₂) #----------------------------- - # Analytic solution + # Analytic solution #----------------------------- "Computes displacements numeric solution uᵢ, uⱼ and uₖ for analytic validation." function u_ijk_numeric(numerical_α::Vector{<:Real}, numerical_β::Vector{<:Real}, @@ -261,11 +261,11 @@ function run_uniaxial_compression() analytic_ℙⱼⱼ(α, β, μ, K) end # Compute the analytic Second Piola-Kirchoff stress tensor `ℙ` for the numeric vectors α and β - # Case 1 + # Case 1 ℙᵢᵢ_analytic_case₁ = analytic_ℙᵢᵢ(numeric_α_case₁, numeric_β_case₁) ℙⱼⱼ_analytic_case₁ = analytic_ℙⱼⱼ(numeric_α_case₁, numeric_β_case₁) ℙₖₖ_analytic_case₁ = analytic_ℙₖₖ(numeric_α_case₁, numeric_β_case₁) - # Case 2 + # Case 2 ℙᵢᵢ_analytic_case₂ = analytic_ℙᵢᵢ(numeric_α_case₂, numeric_β_case₂) ℙⱼⱼ_analytic_case₂ = analytic_ℙⱼⱼ(numeric_α_case₂, numeric_β_case₂) ℙₖₖ_analytic_case₂ = analytic_ℙₖₖ(numeric_α_case₂, numeric_β_case₂) @@ -274,7 +274,7 @@ function run_uniaxial_compression() #-------------------------------- rand_point = [[rand() * Lᵢ, rand() * Lⱼ, rand() * Lₖ]] eval_handler_rand = PointEvalHandler(mesh(s₂), rand_point) - # Compute analytic solution at a random point + # Compute analytic solution at a random point uᵢ_case₂, uⱼ_case₂, uₖ_case₂ = u_ijk_numeric(numeric_α_case₂, numeric_β_case₂, numeric_γ_case₂, rand_point[]...) rand_point_uᵢ = displacements(states_sol_case₂, eval_handler_rand, 1) @@ -282,7 +282,7 @@ function run_uniaxial_compression() rand_point_uₖ = displacements(states_sol_case₂, eval_handler_rand, 3) stress_point = stress(states_sol_case₂, eval_handler_rand)[] #----------------------------- - # Test boolean for CI + # Test boolean for CI #----------------------------- @testset "Case 1 Uniaxial Compression Example" begin @test ℙᵢᵢ_analytic_case₁ ≈ ℙᵢᵢ_numeric_case₁ rtol = RTOL diff --git a/examples/uniaxial_extension/uniaxial_extension.jl b/examples/uniaxial_extension/uniaxial_extension.jl index 2fc7d2db7..8bd148f92 100644 --- a/examples/uniaxial_extension/uniaxial_extension.jl +++ b/examples/uniaxial_extension/uniaxial_extension.jl @@ -71,11 +71,11 @@ function run_uniaxial_extension() # ------------------------------- # Fixed dofs bc₁_label = "fixed-ux" - bc₁ = FixedDof([:u], [1], bc₁_label) + bc₁ = FixedDof(:u, [1], bc₁_label) bc₂_label = "fixed-uj" - bc₂ = FixedDof([:u], [2], bc₂_label) + bc₂ = FixedDof(:u, [2], bc₂_label) bc₃_label = "fixed-uk" - bc₃ = FixedDof([:u], [3], bc₃_label) + bc₃ = FixedDof(:u, [3], bc₃_label) # Load bc₄_label = "tension" bc₄ = GlobalLoad([:u], t -> [p * t, 0, 0], bc₄_label) diff --git a/examples/von_misses_truss/von_misses_truss.jl b/examples/von_misses_truss/von_misses_truss.jl index 00d71b854..d826be856 100644 --- a/examples/von_misses_truss/von_misses_truss.jl +++ b/examples/von_misses_truss/von_misses_truss.jl @@ -1,4 +1,4 @@ -# ------------------------------------------------------------- +# ------------------------------------------------------------- # Von Misses Truss Example from (Zerpa, Bazzano 2017 ) - 2.5.4 # ------------------------------------------------------------- using Test, LinearAlgebra @@ -12,8 +12,8 @@ function run_von_misses_truss_example() ν = 0.0 # Poisson's modulus A₀ = 2.5e-3 # Cross-section area in m² ANG = 65 # truss angle in degrees - L = 2 # Length in m - V = L * cos(deg2rad(ANG)) # vertical distance in m + L = 2 # Length in m + V = L * cos(deg2rad(ANG)) # vertical distance in m H = L * sin(deg2rad(ANG)) # horizontal distance in m Fₖ = -3e8 # Vertical load in N RTOL = 1e-4 # Relative tolerance for tests @@ -31,7 +31,7 @@ function run_von_misses_truss_example() s₁ = Circle(d) a = sqrt(A₀) s₂ = Square(a) - ## Entities + ## Entities truss₁ = Truss(n₁, n₂, s₁, strain_model, "left_truss") # [n₁, n₂] truss₂ = Truss(n₂, n₃, s₂, strain_model, "right_truss") # [n₂, n₃] elements = [truss₁, truss₂] @@ -52,9 +52,9 @@ function run_von_misses_truss_example() # Boundary conditions # ------------------------------- # Fixed dofs - bc₁ = FixedDof([:u], [1, 2, 3], "fixed_uₓ_uⱼ_uₖ") - bc₂ = FixedDof([:u], [2], "fixed_uⱼ") - # Load + bc₁ = FixedDof(:u, [1, 2, 3], "fixed_uₓ_uⱼ_uₖ") + bc₂ = FixedDof(:u, [2], "fixed_uⱼ") + # Load bc₃ = GlobalLoad([:u], t -> [0, 0, Fₖ * t], "load in j") node_bc = dictionary([bc₁ => [n₁, n₃], bc₂ => [n₂], bc₃ => [n₂]]) s_boundary_conditions = StructuralBoundaryCondition(; node_bcs=node_bc) @@ -89,14 +89,14 @@ function run_von_misses_truss_example() @test norm(numerical_uⱼ) ≤ RTOL numerical_λᵥ = -load_factors(sa) * Fₖ # TODO: fix me - # Test stress and strains + # Test stress and strains σ_truss₂ = stress(states_sol, truss₂) ϵ_truss₂ = strain(states_sol, truss₂) # if strain_model == RotatedEngineeringStrain # @test σ_truss₂ == E * ϵ_truss₂ # end #----------------------------- - # Analytic solution + # Analytic solution #----------------------------- "Analytic load factor solution for the displacement `uₖ` towards z axis at node `n₂` `RotatedEngineeringStrain` ." function load_factors_analytic(uₖ::Real, ::Type{RotatedEngineeringStrain}, @@ -115,7 +115,7 @@ function run_von_misses_truss_example() end analytics_λᵥ = load_factors_analytic.(numerical_uₖ, strain_model) #----------------------------- - # Test boolean for CI + # Test boolean for CI #----------------------------- @test analytics_λᵥ ≈ numerical_λᵥ rtol = RTOL end diff --git a/src/BoundaryConditions/FixedDofBoundaryConditions.jl b/src/BoundaryConditions/FixedDofBoundaryConditions.jl index ac3a6990a..a7c312321 100644 --- a/src/BoundaryConditions/FixedDofBoundaryConditions.jl +++ b/src/BoundaryConditions/FixedDofBoundaryConditions.jl @@ -16,38 +16,35 @@ using ..Utils export FixedDof, components """ -Fixed displacement boundary condition. +Fixed boundary condition. -This is a particular instance of the struct `DisplacementBoundaryCondition` -considering null dof value at an specific component of the dof displacements. +Considers null dof values at specific component(s) of the given field. """ -Base.@kwdef struct FixedDof <: AbstractDirichletBoundaryCondition - "Symbols where the where the boundary condition is subscribed." - dofs::Vector{Field} = [:u] - "Vectors of integer indicating the fixed degree of freedom component." - components::Vector{Dof} +struct FixedDof <: AbstractDirichletBoundaryCondition + "Field where the boundary condition is subscribed." + field::Field + "Components of the field which are fixed." + components::Vector{Int64} "Label of the boundary condition." - name::Label = NO_LABEL + name::Label + function FixedDof(dofs::Field, components::Vector{Int64}, name::Label=NO_LABEL) + new(dofs, components, name) + end end "Return the fixed components of the `Dof`s defined in the boundary condition `bc`." components(bc::FixedDof) = bc.components "Return fixed `Dof`s of an `AbstractNode` imposed in the `FixedDof` `fbc`." -function apply(fbc::FixedDof, n::AbstractNode) - fbc_dofs_symbols = dofs(fbc) - dofs_to_delete = Dof[] - for dof_symbol in fbc_dofs_symbols - push!(dofs_to_delete, getindex(dofs(n), dof_symbol)[components(fbc)]...) - end - dofs_to_delete +function apply(bc::FixedDof, n::AbstractNode) + # TODO Rename method to fixed_dofs ? + dofs(n, bc.field)[bc.components] end "Return fixed `Dof`s of an `AbstractFace` or `AbstractElement` imposed in the `FixedDof` `fbc`." -function apply(fbc::FixedDof, e::E) where {E<:Union{AbstractFace,AbstractElement}} - dofs_to_delete = Dof[] - [push!(dofs_to_delete, apply(fbc, n)...) for n in nodes(e)] - dofs_to_delete +function apply(bc::FixedDof, e::E) where {E<:Union{AbstractFace,AbstractElement}} + # TODO Rename method to fixed_dofs ? + reduce(vcat, apply(bc, n) for n in nodes(e)) end end diff --git a/src/Entities/Nodes.jl b/src/Entities/Nodes.jl index 826c662c8..0008cce91 100644 --- a/src/Entities/Nodes.jl +++ b/src/Entities/Nodes.jl @@ -61,7 +61,7 @@ Base.getindex(n::AbstractNode, i::Int) = n.x[i] "Return node degrees of freedom." dofs(n::AbstractNode) = n.dofs -dofs(vn::Vector{<:AbstractNode}) = vcat(dofs.(vn)...) +dofs(vn::Vector{<:AbstractNode}) = vcat(dofs.(vn)...) # mapreduce(dofs, vcat, vn) dofs(n::AbstractNode, s::Field) = n.dofs[s] "Sets a vector of dofs `vd` to the node assigned to the field `s`." diff --git a/test/boundary_conditions.jl b/test/boundary_conditions.jl index 977a2457b..88e5554bf 100644 --- a/test/boundary_conditions.jl +++ b/test/boundary_conditions.jl @@ -11,7 +11,7 @@ using ONSAS.Nodes using ONSAS.TriangularFaces using ONSAS.Tetrahedrons -# Entities +# Entities n₁ = Node(0, 0, 0, dictionary([:u => [Dof(1), Dof(2), Dof(3)], :θ => [Dof(13), Dof(14), Dof(15)], :T => [Dof(25)]])) @@ -29,29 +29,38 @@ t_face = TriangularFace(n₁, n₂, n₃) tetra = Tetrahedron(n₁, n₂, n₃, n₄) @testset "ONSAS.BoundaryConditions.FixedDof" begin + # Generic labeled boundary condition. + generic_fixed_dofs = :u + fixed_components = [1, 3] + generic_bc_label = :fixed_bc_generic + fixed_bc = FixedDof(generic_fixed_dofs, fixed_components, generic_bc_label) + + @test components(fixed_bc) == fixed_components + @test label(fixed_bc) == generic_bc_label + + @test apply(fixed_bc, n₃) == [Dof(7), Dof(9)] + @test apply(fixed_bc, t_face) == [Dof(1), Dof(3), Dof(4), Dof(6), Dof(7), Dof(9)] + @test apply(fixed_bc, tetra) == + [Dof(1), Dof(3), Dof(4), Dof(6), Dof(7), Dof(9), Dof(10), Dof(12)] - # Generic labeled boundary condition - generic_fixed_dofs = [:u, :θ] - fixed_components = [1, 3] # Fixes first, second and third component of the corresponding dofs + generic_fixed_dofs = :θ + fixed_components = [1, 3] generic_bc_label = :fixed_bc_generic fixed_bc = FixedDof(generic_fixed_dofs, fixed_components, generic_bc_label) - @test dofs(fixed_bc) == generic_fixed_dofs @test components(fixed_bc) == fixed_components @test label(fixed_bc) == generic_bc_label - @test apply(fixed_bc, n₃) == [Dof(7), Dof(9), Dof(19), Dof(21)] - @test apply(fixed_bc, t_face) == [Dof(1), Dof(3), Dof(13), Dof(15), Dof(4), - Dof(6), Dof(16), Dof(18), Dof(7), Dof(9), Dof(19), Dof(21)] - @test apply(fixed_bc, tetra) == [Dof(1), Dof(3), Dof(13), Dof(15), Dof(4), Dof(6), Dof(16), - Dof(18), Dof(7), Dof(9), Dof(19), Dof(21), Dof(10), Dof(12), Dof(22), Dof(24)] + @test apply(fixed_bc, n₃) == [Dof(19), Dof(21)] + @test apply(fixed_bc, t_face) == [Dof(13), Dof(15), Dof(16), Dof(18), Dof(19), Dof(21)] + @test apply(fixed_bc, tetra) == + [Dof(13), Dof(15), Dof(16), Dof(18), Dof(19), Dof(21), Dof(22), Dof(24)] end t_to_test = 2.0 @testset "ONSAS.BoundaryConditions.GlobalLoad" begin - - # Generic labeled global load boundary condition + # Generic labeled global load boundary condition. dofs_toapply_bc = [:u, :θ] load_fact_generic(t) = [sin(t), t, t^2] generic_values = t -> load_fact_generic(t) .* [1, 1, 1] @@ -63,17 +72,17 @@ t_to_test = 2.0 @test label(generic_bc) == generic_bc_label @test generic_bc(t_to_test) == values(generic_bc)(t_to_test) - # Node force computation + # Node force computation loaded_dofs, f_vec = apply(generic_bc, n₁, t_to_test) @test loaded_dofs == [Dof(1), Dof(2), Dof(3), Dof(13), Dof(14), Dof(15)] @test f_vec == repeat(generic_bc(t_to_test), 2) - # Face tension computation + # Face tension computation loaded_dofs, p_vec = apply(generic_bc, t_face, t_to_test) @test loaded_dofs == vcat(Dof.(1:9), Dof.(13:21)) @test p_vec == repeat(generic_bc(t_to_test) * area(t_face) / 3, 6) - # Volume tension computation + # Volume tension computation loaded_dofs, b_vec = apply(generic_bc, tetra, t_to_test) @test loaded_dofs == vcat(Dof.(1:24)) @@ -94,7 +103,7 @@ end @test generic_bc(t_to_test) == values(generic_bc)(t_to_test) @test label(generic_bc) == generic_bc_label - # Face tension computation + # Face tension computation loaded_dofs, p_vec = apply(generic_bc, t_face, t_to_test) @test loaded_dofs == vcat(Dof.(1:9)) @test p_vec == -repeat([generic_values(t_to_test) * area(t_face) / 3, 0, 0], 3) diff --git a/test/static_analyses.jl b/test/static_analyses.jl index edf09fca6..66789b078 100644 --- a/test/static_analyses.jl +++ b/test/static_analyses.jl @@ -31,6 +31,7 @@ L = 2 # Length in m d = L * cos(deg2rad(65)) # vertical distance in m h = L * sin(deg2rad(65)) Fₖ = -3e8 # vertical load in N + # ------------------------------- # Materials # ------------------------------- @@ -69,8 +70,8 @@ s_materials = StructuralMaterial(mat_dict) # Boundary conditions # ------------------------------- # Fixed dofs -bc₁ = FixedDof([:u], [1, 2, 3], "fixed_uₓ_uⱼ_uₖ") -bc₂ = FixedDof([:u], [2], "fixed_uⱼ") +bc₁ = FixedDof(:u, [1, 2, 3], "fixed_uₓ_uⱼ_uₖ") +bc₂ = FixedDof(:u, [2], "fixed_uⱼ") # Load bc₃ = GlobalLoad([:u], t -> [0, 0, Fₖ * t], "load in j") node_bcs = dictionary([bc₁ => [n₁, n₃], bc₂ => [n₂], bc₃ => [n₂]]) diff --git a/test/structural_model.jl b/test/structural_model.jl index e3bbceb33..cd5118ca1 100644 --- a/test/structural_model.jl +++ b/test/structural_model.jl @@ -54,11 +54,11 @@ s_materials = StructuralMaterial(mat_dict) Fⱼ = 20.0 Fᵢ = 10.0 dof_dim = 3 -bc₁ = FixedDof([:u], collect(1:dof_dim), "fixed_uₓ_uⱼ_uₖ") -bc₂ = FixedDof([:u], [2], "fixed_uⱼ") +bc₁ = FixedDof(:u, collect(1:dof_dim), "fixed_uₓ_uⱼ_uₖ") +bc₂ = FixedDof(:u, [2], "fixed_uⱼ") bc₃ = GlobalLoad([:u], t -> [0, Fⱼ * t, 0], "load in j") bc₄ = GlobalLoad([:u], t -> [Fᵢ * sin(t), 0, 0], "load in i") -bc₅ = FixedDof([:T], [1], "fixed_T") +bc₅ = FixedDof(:T, [1], "fixed_T") node_bc = dictionary([bc₁ => [n₁, n₃], bc₂ => [n₂], bc₃ => [n₂, n₁]]) face_bc = dictionary([bc₃ => [face₁], bc₅ => [face₁]]) elem_bc = dictionary([bc₄ => [truss₁, truss₂]]) From 196f31a4180ce26cbe2e160e75392ecb53838c1d Mon Sep 17 00:00:00 2001 From: Mauricio Vanzulli <50339940+mvanzulli@users.noreply.github.com> Date: Sun, 18 Jun 2023 12:10:01 +0200 Subject: [PATCH 2/2] Update benchmark/uniaxial_compression/bench_uniaxial_compression.jl --- benchmark/uniaxial_compression/bench_uniaxial_compression.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/benchmark/uniaxial_compression/bench_uniaxial_compression.jl b/benchmark/uniaxial_compression/bench_uniaxial_compression.jl index 521afb3b3..210fcbd7a 100644 --- a/benchmark/uniaxial_compression/bench_uniaxial_compression.jl +++ b/benchmark/uniaxial_compression/bench_uniaxial_compression.jl @@ -8,7 +8,7 @@ function strain_energy_neo(𝔼::AbstractMatrix, K::Real, μ::Real) # First invariant I₁ = tr(ℂ) # Strain energy function - return Ψ = μ / 2 * (I₁ - 2 * log(J)) + K / 2 * (J - 1)^2 + Ψ = μ / 2 * (I₁ - 2 * log(J)) + K / 2 * (J - 1)^2 end # Include `create_mesh` function.