onsas / onsas.jl Goto Github PK

View Code? Open in Web Editor NEW

6.0 6.0 0.0 16.16 MB

An Open Nonlinear Structural Analysis Solver in Julia

License: MIT License

Julia 99.46% Shell 0.11% Makefile 0.43%

computational-mechanics finite-element-method julia structural-engineering

onsas.jl's Introduction

ONSAS: an Open Nonlinear Structural Analysis Solver for GNU-Octave/Matlab

About ONSAS

ONSAS is a GNU-Octave/Matlab code for static/dynamic and linear/non-linear analysis of structures. The first version was developed for educational purposes and was published in a handbook of the course Análisis no lineal de Estructuras taught at Facultad de Ingeniería, Universidad de la República since 2017.

Some examples

A deployable ring

A propeller model

A truss tower model

Contact

You can publicly post in the discussion section or contact privately sending an e-mail to jorgepz [AT] fing.edu.uy .

License

The code is distributed under a GNU-GPL 3.0 license.

How to use ONSAS

The user should follow these steps to install and run onsas:

Download and install GNU-Octave and Paraview
Download the ONSAS zip source files from the realeases web https://github.com/ONSAS/ONSAS/releases
Open GNU-Octave, move to the examples directory and run one of the examples.

An introduction to using and contributing to ONSAS was presented in 2022. The recording is available at this youtube video.

Contributions

The complete list of authors of code, contributions, affiliations and acknowledgments is available in the documentation.

onsas.jl's People

Contributors

Stargazers

Watchers

onsas.jl's Issues

Update UML draw

Add LazySets as optional dep

Can be useful to add Mesh -> LazySet conversions.

using ONSAS.StaticAnalyses


N = 1000 # Number of elements.
E = 30e6 # Young's modulus.
ν = 0.3
ρ = 7.3e-4 # Density.
L = 200 # L/N is the individual element length.
A = 1 # Cross section.

# -------------------------------
# Materials
# -------------------------------
steel = SVK(E, ν, ρ, "steel")

nodes = [Node(l) for l in LinRange(0, L, N + 1)]
S = Square(sqrt(A))
elems = [Truss(nodes[i], nodes[i+1], S) for i in 1:N]
mesh = Mesh(nodes, elems)
nodes

# Plot con LazySets.jl
using ONSAS: Elements, Meshes

# One-dimensional for now.
function Base.convert(::Type{LazySet}, mesh::Mesh{1})
    arr = []
    for node in Elements.nodes(mesh)
        push!(arr, Singleton(coordinates(node)))
    end
    for elem in Meshes.elements(mesh)
        coords = coordinates(elem)
        push!(arr, Interval(first(coords[1]), first(coords[2])))
    end
    return UnionSetArray{Float64,AbstractHyperrectangle{Float64}}(arr)
end
s = convert(LazySet, mesh)
plot(s)

Split `Geometries` module into `Meshes` and `CrossSections`

Add Interpolation order into `AbstractElement` parameters

Define a new simplified input interface

Thanks to @jorgepz and @mforets a new interface for defining input structs will be used:

## scalar parameters
E = 2e11  # Young modulus in Pa
ν = 0.0  # Poisson's modulus
A = 5e-3  # Cross-section area in m^2
ang = 65 # truss angle in degrees
L = 2 # Length in m 
d = L * cos(deg2rad(65))   # vertical distance in m
h = L * sin(deg2rad(65))
# Fx = 0     # horizontal load in N
Fⱼ = -3e8  # vertical   load in N
# -------------------------------
# Materials
# -------------------------------
steel = SVK(E, ν, "steel")
aluminum = SVK(E / 3, ν, "aluminium")
# -------------------------------
# Geometries
# -------------------------------
## Cross section
a = sqrt(A)
s₁ = Square(a, "S1")
s₂ = Square(2a, "A2")
# -------------------------------
# Create mesh
# -------------------------------
## Nodes
n₁ = Node(0.0, 0.0, 0.0)
n₂ = Node(d, h, 0.0)
n₃ = Node(2d, 0.0, 0.0)
vec_nodes = [n₁, n₂, n₃]
## Elements 
truss₁ = Truss(n₁, n₂, s₁) # [n₁, n₂]
truss₂ = Truss(n₂, n₃, s₂) # [n₂, n₃]
vec_elems = [truss₁, truss₂]
s_mesh = Mesh(vec_nodes, vec_elems)
# -------------------------------
# Materials
# -------------------------------
s_materials = StructuralMaterials(
    dictionary([steel => [truss₁], aluminum => [truss₂]])
)
# -------------------------------
# Boundary conditions
# -------------------------------
bc₁ = PinnedDisplacementBoundaryCondition("fixed")
bc₂ = FⱼLoadBoundaryCondition(Fⱼ, "load in j")
node_bc = dictionary([bc₁ => [n₁, n₃], bc₂ => [n₂]])
elem_bc = dictionary([bc₁ => [truss₁], bc₂ => [truss₂]])
s_boundary_conditions = StructuralBoundaryConditions(node_bc, elem_bc)

for (bc, n) in pairs(node_bc)
    apply!(s, bc, n)
end
# -------------------------------
# Create Structure
# -------------------------------
s = Structure(s_mesh, s_materials, s_boundary_conditions)

# -------------------------------
# Structural Analysis
# -------------------------------
# Final load factor
λ₁ = 10
NSTEPS = 9
s_analysis = StaticAnalysis(s, λ₁, NSTEPS=NSTEPS)
s_state = current_state(s_analysis)
@test norm(displacements(s_state)) == 0
@test norm(external_forces(s_state)) == 0
@test norm(internal_forces(s_state)) == 0
# -------------------------------
# Solve analysis
# -------------------------------
tol_f = 1e-10;
tol_u = 1e-10;
max_iter = 100;
tols = ConvergenceSettings(tol_f, tol_u, max_iter)
nr = NewtonRaphson(tols)
sol = solve(s_analysis, nr)

Add SVK material model 🪨

Include an AbstractMaterial supertype and an SVK concrete type.

Add `NodeIndex` and `ElementIndex` structs

Use meta to define repeatedly methods with symbolic patterns

Create access with @eval :

# Composed accessors 
#TODO:fixme
 for m in methodswith(AbstractBoundaryCondition)
     @eval $m(fbc::FixedDisplacementBoundaryCondition) = $m(fbc.bc)
 end

Instead of:

dofs(pbc::PinnedDisplacementBoundaryCondition) = dofs(pbc.bc)
Base.values(pbc::PinnedDisplacementBoundaryCondition) = values(pbc.bc)
label(pbc::PinnedDisplacementBoundaryCondition) = label(pbc.bc)

Another option for composition is:

using Lazy

@forward ParentObject.attribute method_1, method_2,...

Modules to check

Add constructors inside structs

This issue seek to add inner constructors for the sake of being more neat. Structs allows to define inner constructors into interface.jl like this:

julia> struct OurRational{T<:Integer} <: Real
           num::T
           den::T
           function OurRational{T}(num::T, den::T) where T<:Integer
               if num == 0 && den == 0
                    error("invalid rational: 0//0")
               end
               num = flipsign(num, den)
               den = flipsign(den, den)
               g = gcd(num, den)
               num = div(num, g)
               den = div(den, g)
               new(num, den)
           end
       end

Add analytic solution for von misses truss example

Fix global indices issue

Running elements.jl twice fails, probably due to some global being defined which is not reset.

Reorganize old files and update them to the new interface version

Obviously old files and valuable contributions form @jorgepz and @JoaquinViera will be included from this commit to solve canonical problems with the "new" interface.

Together with @mforets we think that by the moment (for the sake of clarity) is better to include only files that are being used. Are u agree ?

Remove `AbstractType` from concrete type structs

Change this:

struct Truss{dim,G<:AbstractCrossSection,T<:Real} <: AbstractElement{dim,T}
    n₁::AbstractNode{dim,T}
    n₂::AbstractNode{dim,T}
    cross_section::G
    label::Symbol
    function Truss(n₁::AbstractNode{dim,T}, n₂::AbstractNode{dim,T}, g::G, label=:no_labelled_elem) where
    {dim,G<:AbstractCrossSection,T<:Real}
        new{dim,G,T}(n₁, n₂, g, Symbol(label))
    end
end

for

struct Truss{dim,N<:AbstractNode{dim},G<:AbstractCrossSection,T<:Real} <: AbstractElement{dim,T}
    n₁::N
    n₂::N
    cross_section::G
    label::Symbol
    function Truss(n₁::N, n₂::N, g::G, label=:no_labelled_elem) where
    {dim,T<:Real,N<:AbstractNode{dim,T},G<:AbstractCrossSection}
        new{dim,N,G,T}(n₁, n₂, g, Symbol(label))
    end
end

This is more idiomatic for dispatch. However, in cases that the filed type is changed during the execution (for instance using singletons) an abstract type will be used (or ScalarWrapper) trick.

Delete symbol from `Dof` attribute

Apply displacements and load boundary conditions

Update julia version in to `1` (lastest) in the C.I 📟

Fix strain and stresses methods from `StructuralState`

Remove intermediate iterations form strain and state fields

julia> strain(states_sol[1])
8-element Vector{Any}:
  0.0
  0.0
 -0.03141352193657888
 -0.03141352193657888
 -0.031756065363709
 -0.031756065363709
 -0.03175610984877188
 -0.03175610984877188

Run examples in CI

Consider to add `AbsrtactEntity` for common methods along `Node`s, `Element`s and `Face`s

Maybe common methods such as:

dimension(::AbstractEntity)
coordinates(::AbstractEntity)
dofs(::AbstractEntity)
nodes(::AbstractEntity)

Can be defined in an AbstractEntity interface, and the a AbstractNode, AbstractFace and AbstractElement can use abstract methods by inheritance.

Test `Geometries` module ⚪ 🟦

Geometries to be added:

An example of a rectangle cross-section:

""" Rectangle cross-section.
### Fields:
- `width_y` -- width in `y` axis.
- `width_z` -- width in `z` axis.
"""
struct Rectangle <: AbstractCrossSection
    width_y::Number
    width_z::Number
end

area(r::Rectangle) = cs.width_y * cs.width_z
function Ixx(r::Rectangle)
    #  torsional constant from table 10.1 from Roark's Formulas for Stress and Strain 7th ed.
    a = 0.5 * max([r.width_y, r.width_z])
    b = 0.5 * min([r.width_y, r.width_z])

    Ixx = a * b^3 * (16 / 3 - 3.36 * b / a * (1 - b^4 / (12 * a^4)))
    return Ixx
end
Iyy(r::Rectangle) = r.width_z^3 * r.width_y / 12
Izz(r::Rectangle) = r.width_y^3 * r.width_z / 12
Ixy(r) = 0.0
Ixz(r) = 0.0
Iyz(r) = 0.0

Add a mesh implementation

AbstractMesh
Mesh

Implement a `TriangularFace` element

According to the discussion with @jorgepz this element will set with no material only for processing boundary conditions on a surface.

Add linear dynamic analysis

Should assemble the matrices M, C, K given a problem formulation like that in #97

Can be composed with https://github.com/ONSAS/StructuralDynamicsODESolvers.jl for the dynamic analysis.

Add progress meter

May be useful in ONSAS.jl @jorgepz

Progress meter

Add `StructuralAnalysisCase` and `StructuralSolver` interface

Update `README.md`

Implement `Triangle` element

The triangle element should be comparable with the Triangle defined in ONSAS considering a material field.

Element nodes storage

Thanks to @mforets comment in #78 we can use SVector{<:AbstractNode} to store the nodes of each element.

Change this:

"""
A `Truss` represents an element composed by two `Node`s that transmits axial force only.
### Fields:
- `n₁`             -- stores first truss node.
- `n₂`             -- stores second truss node.
- `cross_sections` -- stores the truss cross-section properties.
- `label`          -- stores the truss label.
"""
struct Truss{dim,N<:AbstractNode{dim},G<:AbstractCrossSection,T<:Real} <: AbstractElement{dim,T}
    n₁::N
    n₂::N
    cross_section::G
    label::Symbol
    function Truss(n₁::N, n₂::N, g::G, label=:no_labelled_element) where
    {dim,T<:Real,N<:AbstractNode{dim,T},G<:AbstractCrossSection}
        new{dim,N,G,T}(n₁, n₂, g, Symbol(label))
    end
end

For this:

"""
A `Truss` represents an element composed by two `Node`s that transmits axial force only.
### Fields:
- `n₁`             -- stores first truss node.
- `n₂`             -- stores second truss node.
- `cross_sections` -- stores the truss cross-section properties.
- `label`          -- stores the truss label.
"""
struct Truss{dim,N<:AbstractNode{dim},G<:AbstractCrossSection,T<:Real} <: AbstractElement{dim,T}
    nodes::SVector{N}
    cross_section::G
    label::Symbol
    function Truss(n₁::N, n₂::N, g::G, label=:no_labelled_element) where
    {dim,T<:Real,N<:AbstractNode{dim,T},G<:AbstractCrossSection}
        new{dim,N,G,T}(n₁, n₂, g, Symbol(label))
    end
end
# Node acccesor:
nodes(t::Truss) = t.nodes

Implement the same function to compute forces vectors and tangent matrices

This violates single responsibility principle but leverages performance.

Features by example

Von Misses Truss example:

Materials: SVK
CrossSections: Circle, Square
Elements: Truss, Node, Dof
BoundaryConditions: GlobalLoadBoundaryCondition FixedDofBoundaryCondition
Meshes: Mesh
StructuralModel: StructuralMaterials, StructuralBoundaryConditions (only applied into nodes)
StructuralSolvers: NewtonRaphson
StructuralAnalyses: StaticAnalysis
Analytic solution

Uniaxial extension:

Materials: SVK
Materials: NeoHookean
Elements: Tetrahedron, Node, Dof, TriangularFace
BoundaryConditions: GlobalLoadBoundaryCondition FixedDofBoundaryCondition
Meshes: Mesh
Meshes: GMSHMesh
StructuralSolvers: NewtonRaphson
StructuralAnalyses: StaticAnalysis

Create `BoundaryConditions` module and interface

The following boundary conditions structs are included:

Fix overflow issue in x86 SO for random tests

ONSAS.StructuralAnalyses.StaticAnalyses.StaticState: Test Failed at /home/runner/work/ONSAS.jl/ONSAS.jl/test/static_analysis.jl:139
  Expression: ≈((internal_forces(default_s))[1:6], fᵢₙₜ_e_1, rtol = RTOL)
   Evaluated: [0.7863760384118906, -4.7362539917743816e166, -7.009402958251185e166, -7.009042741363477e166, -6.023890074008387e208, -1.0110027777808359e235] ≈ [0.7863760384118906, 0.0996784825685848, 0.7665694221672829, 0.13643254313390907, 0.817823788943396, 0.32092303565108105] (rtol=0.05)

Stacktrace:
 [1] macro expansion
   @ /opt/hostedtoolcache/julia/nightly/x86/share/julia/stdlib/v1.10/Test/src/Test.jl:478 [inlined]
 [2] macro expansion
   @ ~/work/ONSAS.jl/ONSAS.jl/test/static_analysis.jl:139 [inlined]
 [3] macro expansion
   @ /opt/hostedtoolcache/julia/nightly/x86/share/julia/stdlib/v1.10/Test/src/Test.jl:1504 [inlined]
 [4] top-level scope
   @ ~/work/ONSAS.jl/ONSAS.jl/test/static_analysis.jl:80
ONSAS.StructuralAnalyses.StaticAnalyses.StaticState: Test Failed at /home/runner/work/ONSAS.jl/ONSAS.jl/test/static_analysis.jl:141
  Expression: ≈((internal_forces(default_s))[1:6], fᵢₙₜ_e_1, rtol = RTOL)
   Evaluated: [0.7863760384118906, -4.7362539917743816e166, -7.009402958251185e166, -7.009042741363477e166, -6.023890074008387e208, -1.0110027777808359e235] ≈ [0.7863760384118906, 0.0996784825685848, 0.7665694221672829, 0.13643254313390907, 0.817823788943396, 0.32092303565108105] (rtol=0.05)

Stacktrace:
 [1] macro expansion
   @ /opt/hostedtoolcache/julia/nightly/x86/share/julia/stdlib/v1.10/Test/src/Test.jl:478 [inlined]
 [2] macro expansion
   @ ~/work/ONSAS.jl/ONSAS.jl/test/static_analysis.jl:141 [inlined]
 [3] macro expansion
   @ /opt/hostedtoolcache/julia/nightly/x86/share/julia/stdlib/v1.10/Test/src/Test.jl:1504 [inlined]
 [4] top-level scope
   @ ~/work/ONSAS.jl/ONSAS.jl/test/static_analysis.jl:80
ONSAS.StructuralAnalyses.StaticAnalyses.StaticState: Test Failed at /home/runner/work/ONSAS.jl/ONSAS.jl/test/static_analysis.jl:163
  Expression: ≈(internal_forces(default_s), Fᵢₙₜ, rtol = RTOL)
   Evaluated: [0.7863760384118906, -4.7362539917743816e166, -7.009402958251185e166, -7.009042741363477e166, -6.023890074008387e208, -1.0110027777808359e235, NaN, 0.7[153](https://github.com/ONSAS/ONSAS.jl/actions/runs/4244459372/jobs/7378604333#step:6:156)354026002874, -1.8870917777[161](https://github.com/ONSAS/ONSAS.jl/actions/runs/4244459372/jobs/7378604333#step:6:164)854e170] ≈ [0.7863760384118906, 0.0996784825685848, 0.766569422[167](https://github.com/ONSAS/ONSAS.jl/actions/runs/4244459372/jobs/7378604333#step:6:170)2829, 0.30965578128694404, 0.8678596494771985, 0.32133320780145946, 0.9850382101024203, 0.7153354026002874, 0.8333443798517911] (rtol=0.05)

Efficient way to compute element dofs

dofs(e::AbstractElement) = mergewith(vcat, dofs.(nodes(e))...)
"Returns the `AbstractElement` `e` dofs."
function dofs(e::AbstractElement)
    vecdfs = dofs.(nodes(e))

Create a cache with element dofs once the Mesh is created.
Loop over nodes(e)

Create a Model Struct module defining interfaces:

New structs to define the old MEBI are used:

StructuralMaterials
StructuralElements
StructuralBoundaryConditions

Add Strain model as a parametric type into `Truss` elements

This can be done like extending the current Truss element:

"""
A `Truss` represents an element composed by two `Node`s that transmits axial force only.
### Fields:
- `n₁`             -- stores first truss node.
- `n₂`             -- stores second truss node.
- `cross_sections` -- stores the truss cross-section properties.
- `label`          -- stores the truss label.
"""
struct Truss{dim,G<:AbstractCrossSection,T<:Real} <: AbstractElement{dim,T}
    n₁::AbstractNode{dim,T}
    n₂::AbstractNode{dim,T}
    cross_section::G
    label::Symbol
    function Truss(n₁::AbstractNode{dim,T}, n₂::AbstractNode{dim,T}, g::G, label=:no_labelled_elem) where
    {dim,G<:AbstractCrossSection,T<:Real}
        new{dim,G,T}(n₁, n₂, g, Symbol(label))
    end
end

struct Truss{E::AbstractStrain,dim,G<:AbstractCrossSection,T<:Real} <: AbstractElement{dim,T}
    n₁::AbstractNode{dim,T}
    n₂::AbstractNode{dim,T}
    cross_section::G
    label::Symbol
    function Truss(n₁::AbstractNode{dim,T}, n₂::AbstractNode{dim,T}, g::G, label=:no_labelled_elem) where
    {dim,G<:AbstractCrossSection,T<:Real}
        new{dim,G,T}(Green(),n₁, n₂, g, Symbol(label))
    end
end

and then compute the strain with

_strain(::Green, l_ini::Number, l_def::Number) = (l_def^2 - l_ini^2) / (l_ini * (l_ini + l_def))

Document `Element`s formulations implemented

Include a `ONSAS.Utils` module :hammer:

This module is used to define first methods (which do not belong to a specific module) and then overload them. For instance the label method should be defined in:

module Utils

" Function that returns the label of an object "
function label() end 

end #  module

Then inside the Geometry module it would be used:

module BoundaryCondition

@reexport import ..Utils: label

label(bc::AbstractBoundaryCondition) = bc.label 
...

The same principle can be applied to Gauss integrations points. This is based on the I of SOLID principle

Add Clamped example (1D)

#=
This model is taken from [1]. The theoretical derivations of the analytic solution can be found in [2].
The implementation of stiffness and mass matrices in Julia can be found in [3].

[1] Malakiyeh, Mohammad Mahdi, Saeed Shojaee, and Klaus-Jürgen Bathe. "The Bathe time integration method revisited for prescribing desired numerical dissipation." Computers & Structures 212 (2019): 289-298.

[2] Mechanical Vibrations, Gerardin et al, page 250-251.

[3] https://github.com/JuliaReach/SetPropagation-FEM-Examples/blob/main/examples/Clamped/Clamped_Model.jl
=#

using ONSAS.StaticAnalyses

N = 1000 # Number of elements.
E = 30e6 # Young's modulus.
ν = 0.3
ρ = 7.3e-4 # Density.
L = 200 # L/N is the individual element length.
A = 1 # Cross section.

# -------------------------------
# Materials
# -------------------------------
steel = SVK(E, ν, ρ, "steel")


nodes = [Node(l) for l in LinRange(0, L, N + 1)]
S = Square(sqrt(A))
elems = [Truss(nodes[i], nodes[i+1], S) for i in 1:N]
mesh = Mesh(nodes, elems)
nodes

# -------------------------------
# Dofs
#--------------------------------
dof_dim = 1
add!(mesh, :u, dof_dim)

Generate Docs with `Documenter.jl`

Add a windows S.O into a continuous integration

Consider to remove 32 bits S.O

Consider to use `Symmetric` from LinearAlgebra

julia> B = [1.5 2 -4; 2 -1 -3; -4 -3 5]
3×3 Matrix{Float64}:
  1.5   2.0  -4.0
  2.0  -1.0  -3.0
 -4.0  -3.0   5.0

julia> sB = Symmetric(B)
3×3 Symmetric{Float64, Matrix{Float64}}:
  1.5   2.0  -4.0
  2.0  -1.0  -3.0
 -4.0  -3.0   5.0

Element matrices
Global structure matrices

Add UML diagram

Add uniaxial extension test

Set author/contributor policies

Discuss and define a author/contributor policy in a md file.

Create an interface and module for elements features

Move `apply` method into `BoundaryConditions` module and dispatch using type of `Element`s

Document and add tests for `StructuralSolvers` module

Return an`AbstractConvergenceCriterion` from `isconverged` function

Change isconverged!(current_iteration(sa), tolerances(alg)) for isconverged(current_iteration(sa), alg) isa NotConvergedYet`

Implement NewtonRaphson algorithm

Validate the method with:

Von misses truss
Uniaxial extension

Implement a `Tetrahedron` element

Delete `DofIndex` struct

As we discussed with @mforets and @jorgepz, for the sake of simplicity, degree of freedom indexes will be just Integers. This is essential to develop a faster assembler. In contrast, element and nodes will keep being indexed with ElementIndex and NodeIndex objects. This is important to define boundary conditions and also to extract information about the element/node from the mesh struct:

mesh[ElementIndex(1)] = elem_1

s_boundary_conditions = StructuralBoundaryConditions(bcs_dict)
bc_sets = Dict(
    "fixed" => Set([NodeIndex(1), NodeIndex(3)]),
    "load" => Set([NodeIndex(2)])
)