Reference

noFB

A constant, used by elements' residual or lagrangian as their 3rd output if they do provide any feedback to the solver (for example, on the reduction of the barrier parameter in interior point method).

Example: return L,noFB

See also: noFB

AbstractElement

An abstract data type. An element type MyElement must be declared as a subtype of AbstractElement.

MyELementmust provide a constructor with interface

`eleobj = MyElement(nod::Vector{Node}; kwargs...)`

See also: coord, Node, lagrangian

Acost{Na,ainod,afield,Tcost,Tcostargs} <: AbstractElement

An element to apply a once-off cost on a combination of A-dofs. For costs per unit of time on A-dofs (not recommended), see DofCost.

Named arguments to the constructor

• inod::NTuple{Na,𝕫}=() For each A-dof to enter cost, its element-node number.
• field::NTuple{Na,Symbol}=() For each A-dof to enter cost, its field.
• cost::Functor cost(A,costargs...)→ℝ
• costargs::NTuple=() or NamedTuple of additional arguments splatted when calling cost.

Requestable internal variables

• cost, the value of the cost.

Example

@functor with() acost(A;A0)=(A[1]-A0)^2
ele1 = addelement!(model,Acost,[nod1],inod=(1,),field=(:EI,),
       cost=acost,costargs=(;A0=0.27)

Note that A0 is within a NamedTuple in the call to addelement!, but not in the definition of the Functor acost.

See also: SingleAcost, DofCost, SingleDofCost, ElementCostAndConstraint, addelement!

DirectXUA{OX,OU,IA}

A non-linear direct solver for optimisation FEM.

An analysis is carried out by a call with the following syntax:

initialstate    = initialize!(model)

The solver does not yet support interior point methods.

Parameters

• OX 0 for static analysis 1 for first order problems in time (viscosity, friction, measurement of velocity) 2 for second order problems in time (inertia, measurement of acceleration)
• OU 0 for white noise prior to the unknown load process 2 otherwise
• IA 0 for XU problems (variables of class A will be unchanged) 1 for XUA problems

Named arguments

• dbg=(;) a named tuple to trace the call tree (for debugging).
• verbose=true set to false to suppress printed output (for testing).
• silenterror=false set to true to suppress print out of error (for testing) .
• initialstate an AbstractVector of State: one initial state for each experiment. initialstate must be with zero time derivatives. It does not provide initial conditions for the problem, but an initial guess for the iterative solver.
• time an AbstractVector (of same length as initialstate) of AbstractRange of times at which to compute the steps. Example: 0:0.1:5.
• maxiter=50 maximum number of Newton-Raphson iterations.
• maxΔλ=1e-5 convergence criteria: a norm of the scaled Λ increment.
• maxΔx=1e-5 convergence criteria: a norm of the scaled X increment.
• maxΔu=1e-5 convergence criteria: a norm of the scaled U increment.
• maxΔa=1e-5 convergence criteria: a norm of the scaled A increment.
• saveiter=false set to true so that the output state contains the states at each Newton-Raphson iteration (for debugging non-convergence).

Setting the following flags to true will improve the sparsity of the system. But setting a flag to true when the condition isn't met causes the Hessian to be wrong, which is detrimental for convergence.

• Xwhite=false true if response measurement error is a white noise process.
• XUindep=false true if response measurement error is independant of U
• UAindep=false true if U is independant of A
• XAindep=false true if response measurement error is independant of A

Output

• state, where state[iexp][itime] contains the state of the optimized model at each of these steps, or if saveiter=true then state[iiter][iexp][itime] is a state.

See also: solve, initialize!, SweepX, FreqXU

DofConstraint{λclass,Nλ,Nx,Nu,Na,
              λinod,λfield, xinod,xfield, uinod,ufield, ainod,afield,
              Tg,Tmode} <: AbstractElement

An element to apply physical/optimisation equality/inequality constraints on dofs.

The constraints are holonomic, i.e. they apply to the values, not the time derivatives, of the involved dofs. This element is very general but not very user-friendly to construct, factory functions are provided for better useability. The sign convention is that each gap g≥0 and each Lagrange multiplier λ≥0.

This element can generate three classes of constraints, depending on the input argument λclass.

• λclass=:X Physical constraint. In mechanics, the Lagrange multiplier dof is a generalized force, dual of the gap. The gap Functor must be of the form gap(x,t,gargs...).
• λclass=:U Time varying optimisation constraint. For example: find A-parameters so that at all times, the response does not exceed a given criteria. The gap Functor must be of the form gap(x,u,a,t,gargs...).
• λclass=:A Time invariant optimisation constraint. For example: find A-parameters such that A[1]+A[2]=gargs.somevalue. The gap Functor must be of the form gap(a,gargs...).

All constraints in one DofConstraint element are of the same λclass. The gap Functor must return a SVector{Nλ,𝕫}.

Named arguments to the constructor

• λclass::Symbol The class (:X,:U or :A) of the Lagrange multipliers.
• λinod ::NTuple{Nλ,𝕫 } The element-nodes number of the Lagrange multipliers.
• λfield::NTuple{Nλ,Symbol} The field of the Lagrange multipliers.
• xinod ::NTuple{Nx,𝕫 }=() For each X-dof to be constrained, its element-node number.
• xfield::NTuple{Nx,Symbol}=() For each X-dof to be constrained, its field.
• uinod ::NTuple{Nu,𝕫 }=() For each U-dof to be constrained, its element-node number.
• ufield::NTuple{Nu,Symbol}=() For each U-dof to be constrained, its field.
• ainod ::NTuple{Na,𝕫 }=() For each A-dof to be constrained, its element-node number.
• afield::NTuple{Na,Symbol}=() For each A-dof to be constrained, its field.
• λinod::𝕫 The element-node number of the Lagrange multiplier.
• λclass::Symbol The class (:X,:U or :A) of the Lagrange multiplier. See the explanation above for classes of constraints
• λfield::Symbol The field of the Lagrange multiplier.
• gap::Functor The gap function.
• gargs::NTuple=() or NamedTuple Additional inputs to the gap function.
• mode::Functor where mode(t::ℝ) -> Symbol, with value :equal, :positive or :off at any time. An :off constraint will set the Lagrange multiplier to zero. Applies to all Nλ constraints.

Example TODO

using Muscade,StaticArrays
model           = Model(:TestModel)
n1              = addnode!(model,𝕣[0])
@functor with() gap(x,    t)=SVector(x+.1)
@functor with() res(x,u,a,t)=0.4x.+.08+.5x.^2)
e1              = addelement!(model,DofConstraint,[n1],λclass=:X, 
                              λinod=(1,),λfield=(:λ1,), 
                              xinod=(1,),xfield=(:t1,),
                              gap=gap,
                              mode=positive)
e2              = addelement!(model,Muscade.QuickFix  ,[n1],inod=(1,),field=(:t1,),
                              res=res
initialstate    = initialize!(model)
setdof!(initialstate,1.;field=:λ1)
state           = solve(SweepX{0};initialstate,time=[0.],verbose=false) 
X               = state[1].X[2]

See also: Hold, ElementCostAndConstraint, off, equal, positive

DofCost{Class,Nx,Nu,Na,xinod,xfield,uinod,ufield,ainod,
    afield,Tcost,Tcostargs} <: AbstractElement

An element to apply costs on combinations of dofs. The cost is "per unit of time". For once-off costs on A-dofs, see Acost.

Named arguments to the constructor

• xinod::NTuple{Nx,𝕫}=() For each X-dof to enter cost, its element-node number.
• xfield::NTuple{Nx,Symbol}=() For each X-dof to enter cost, its field.
• uinod::NTuple{Nu,𝕫}=() For each U-dof to enter cost, its element-node number.
• ufield::NTuple{Nu,Symbol}=() For each U-dof to enter cost, its field.
• ainod::NTuple{Na,𝕫}=() For each A-dof to enter cost, its element-node number.
• afield::NTuple{Na,Symbol}=() For each A-dof to enter cost, its field.
• cost::Functor cost(X,U,A,t,costargs...)→ℝ X and U are tuples (derivates of dofs...), and ∂0(X),∂1(X),∂2(X) must be used by cost to access the value and derivatives of X (resp. U)
• costargs::NTuple=() or NamedTuple of additional arguments splatted when calling cost.

Requestable internal variables

• cost, the value of the cost.

Example

@functor with() xcost(X,U,A,t;X0) = (X[1]-X0)^2
ele1 = addelement!(model,DofCost,[nod1],xinod=(1,),xfield=(:tx1,),
       cost=xcost,costargs=(;X0=0.27)

Note that X0 is within a NamedTuple in the call to addelement!, but not in the definition of the Functor xcost.

See also: Acost, SingleDofCost, ElementCostAndConstraint, addelement!

DofLoad{Tvalue,Field} <: AbstractElement

An element to apply a loading term to a single X-dof.

Named arguments to the constructor

• field::Symbol.
• value::Functor, where value(t::ℝ) → ℝ.
• args::NTuple=() or NamedTuple of additional arguments passed to value.

Requestable internal variables

• F, the value of the load.

Examples

using Muscade
model = Model(:TestModel)
node  = addnode!(model,𝕣[0,0])
@functor with(a=3,b=-1) load(t,c)=a*t^c+b
e     = addelement!(model,DofLoad,[node];field=:tx,value=load,args=(2.,))

Note that 2. is within a Tuple in the call to addelement!, but c is not, in the definition of the Functor load.

See also: Hold, DofCost

eiginc = solve(EigX{ℝ     };state=initialstate,nmod)
eiginc = solve(EigX{:fullℝ};state=initialstate,nmod)
eiginc = solve(EigX{ℂ     };state=initialstate,nmod)

Given an initial (typicaly static) state initialstate, computes the lowest nmod eigenmodes of a system. EigX{ℝ} computes real eigenmodes not accounting for damping. EigX{ℝ} computes complex eigenmodes accounting for damping. The data structure eiginc can be passed to increment to obtain dynamic states superimposing mode shapes.

Input

• initialstate - a State, at which the problem is linearized.
• nmod=5 - the number of eigenmodes to identify
• droptol=1e-9 - in the stiffness and mass matrix, the magnitude of a term relative to the largest term in the matrix under which the term is set to zero.
• Further named arguments: see the optional keyword arguments to geneig.

Output

• an object of type EigXℝincrement or EigXℂincrement, for use with increment to create a snapshot of the oscillating system.

See also: solve, initialize!, increment

EigXU{OX,OU}

Study the combinations of load and response that are least detected by sensor systems.

An analysis is carried out by a call with the following syntax:

initialstate    = initialize!(model)
eiginc          = solve(EigXU{OX,OU};Δω, p, nmod,initialstate)

The solver linearises the problem (computes the Hessian of the Lagrangian) at initialstate and solves the ΛXU-eigenvalue problem at frequencies $ωᵢ = Δω*i$ with $i∈{0,...,2ᵖ-1}$.

Parameters

• OX 0 for static analysis 1 for first OX problems in time (viscosity, friction, measurement of velocity) 2 for second OX problems in time (inertia, measurement of acceleration)
• OU 0 for white noise prior to the unknown load process 2 otherwise

Named arguments

• dbg=(;) a named tuple to trace the call tree (for debugging).
• verbose=true set to false to suppress printed output (for testing).
• initialstate a State, at which the problem is linearized.
• nmod the number of eigen-modes to identusy
• Δω frequency step
• p 2^p steps will be analysed.
• droptol=1e-10 set to zero terms in the incremental matrices that are smaller than droptol in absolute value.

Output

• an object of type EigXUincrement for use with increment to create a snapshot of the oscillating system.

See also: increment,EigXU, solve, initialize!, study_singular, SweepX, DirectXUA

ElementCostAndConstraint{TargetElement,λinod,λfield,Nu,Treq,Tg,Tgargs,Tmode} <: AbstractElement

An element to apply a single costs and multiple physical and optimisation equality and inequality constraints on the element-results of another "target" element. The target element must not be added separately to the model. Instead, the TargetElement, and the named arguments to the target element are provided as input to the ElementCostAndConstraint constructor.

The Lagrangian multipliers introduced by physical and optimisation constraints are of class :X and :U respectively.

Named arguments to the constructor

Omit irrelevant arguments if no physical or optimisation constraints, or cost, are required.

• λxinod::NTuple{Nλx,𝕫} The element-nodes number of the physical Lagrange multipliers.
• λxfield::NTuple{Nλx,Symbol} The fields of the physical Lagrange multipliers.
• λuinod::NTuple{Nλu,𝕫} The element-nodes number of the optimisation Lagrange multipliers.
• λufield::NTuple{Nλu,Symbol} The fields of the optimisation Lagrange multipliers.
• req A request for element-results to be extracted from the target element, see @request. The request is formulated as if adressed directly to the target element.
• gap A gap function gap(eleres,t,gargs...) returning a Tuple of ℝ, with Nλx gaps of the physical constraints followed by Nλu gaps of the optimisation constraints. eleres is the output of the above-mentionned request to the target element.
• gargs::NTuple=() or NamedTuple Additional inputs to the gap function.
• mode::Functor where mode(t::ℝ), returning a Tuple of Symbols, with Nλx :off,:equal or :positive for the physical constraints followed by Nλu similar Symbols for the optimisation constraints.
• cost A cost Functor cost(eleres,t,costargs...)→ℝ. eleres is the output of the above-mentionned request to the target element.
• costargs::NTuple=() or NamedTuple of additional arguments splatted when calling cost.
• TargetElement The named of the constructor for the target element. This cannot be a ElementCostAndConstraint.
• elementkwargs A named tuple containing the named arguments of the TargetElement constructor.

Requestable internal variables

These are variables requestable from ElementCostAndConstraint to extract results after the analysis. This is distinct from the req input described above, which relate to element results to be extracted from the target element, to evaluate cost or gaps.

• λ The constraints Lagrange multiplier
• gap The constraints gap Functor
• mode The vector of modes for each constraint
• eleres(...) where ... is the list of requestables wanted from the target element. The "prefix" eleres is there to prevent possible confusion with variables requestable from ElementCostAndConstraint. For example @request gap would extract the value of the ElementCostAndConstraint's varable gap, while @request eleres(gap) refers to the value of a variable called gap in the target element.

See also: DofCost, Hold, DofConstraint, off, equal, positive, @request, @functor

FreqXU{OX,OU}

A linear frequency domain solver for optimisation FEM.

An analysis is carried out by a call with the following syntax:

initialstate    = initialize!(model)
stateXU         = solve(FreqXU{OX,OU};Δt, p, t₀,tᵣ,initialstate)

The solver linearises the problem (computes the Hessian of the Lagrangian) at initialstate with time tᵣ, and solves it at times t=range(start=t₀,step=Δt,length=2^p). The return

Parameters

• OX 0 for static analysis 1 for first order problems in time (viscosity, friction, measurement of velocity) 2 for second order problems in time (inertia, measurement of acceleration)
• OU 0 for white noise prior to the unknown load process 2 otherwise

Named arguments

• dbg=(;) a named tuple to trace the call tree (for debugging).
• verbose=true set to false to suppress printed output (for testing).
• silenterror=false set to true to suppress print out of error (for testing) .
• initialstate a State at which the problem is linearized.
• t₀=0. time of first step.
• Δt time step.
• p 2^p steps will be analysed.
• tᵣ=t₀ reference time for linearisation.
• droptol=1e-10 set to zero terms in the incremental matrices that are smaller than droptol in absolute value.

Output

A vector of length 2^p containing the state of the model at each of these steps.

See also: solve, initialize!, study_singular, SweepX, DirectXUA

f  = Muscade.FunctionFromVector(xs::AbstractRange,ys::AbstractVector)
y  = f(x)

Linear interpolation.  Fails if `x` is outside the range `xs`

function MyElement(nod::Vector{Node}; ... foo::Functor ...)

@functor returns an object of type Functor. Functor is a subtype of the abstract type Function.

In an element constructor, requiring a function-like input to be a Functor forces the user of the element to use the @functor macro to define functions provided as input to an element. This is designed to safeguard the user from some arguably unintuitive behaviour of closures.

See also: @functor

Hold <: AbstractElement

An element to set a single X-dof to zero.

Named arguments to the constructor

• field::Symbol. The field of the X-dof to constraint.
• λfield::Symbol=Symbol(:λ,field). The field of the Lagrange multiplier.

Example

using Muscade
model = Model(:TestModel)
node  = addnode!(model,𝕣[0,0])
e     = addelement!(model,Hold,[node];field=:tx)

See also: DofConstraint, DofLoad, DofCost

v=McLaurin(Ty,Δx)

Ty::∂ℝ has partials to arbitrary order with respect to a variable x. These partials define a McLaurin expansion, which McLaurin evaluates at value Δx.

Thus the elements in v have the same type as Δx

McLaurin handles nested structures of Tuples and SVectors of ∂ℝ, applying the expansion to each element.

McLaurin is a utility function behind chainrule and Taylor

See also: chainrule, Taylor, revariate, apply

model = Model([ID=:my_model])

Construct a blank model, which will be mutated to create a FEM-constrained optimisation problem.

See also: addnode!, addelement!, describe, solve

Muscade.Monitor <: AbstractElement

An element for for monitoring inputs to and outputs from another element, during an analysis.

Instead of adding the element to be monitored directly into the model, add this element with the element to be monitored as argument.

Inputs and outputs are @show'n.

Named arguments to the constructor

• TargetElement The the type of element to be monitored-
• trigger A function that takes dbg as an input and returns a boolean (true) to printout.
• elementkwargs a NamedTuple containing the named arguments of the TargetElement constructor.

Node

The eltype of vectors handed by Muscade as first argument to element constructors.

Example: function SingleDofCost(nod::Vector{Node};class::Symbol, ... )

See also: coord

Muscade.QuickFix <: AbstractElement

An element for creating simple elements with "one line" of code. Elements thus created have several limitations:

• physical elements with only X-dofs.
• only R can be espied.

The element is intended for testing. Muscade-based applications should not include this in their API.

Named arguments to the constructor

• inod::NTuple{Nx,𝕫}. The element-node numbers of the X-dofs.
• field::NTuple{Nx,Symbol}. The fields of the X-dofs.
• res::Functor, where res(X::ℝ1,X′::ℝ1,X″::ℝ1,t::ℝ) → ℝ1, the residual.

Examples

A one-dimensional linear elastic spring with stiffness 2.

using Muscade
model = Model(:TestModel)
node1  = addnode!(model,𝕣[0])
node2  = addnode!(model,𝕣[1])
@functor with() res(x,u,a,t)=0.4x.+.08+.5x.^2) 
e = addelement!(model,Muscade.QuickFix,[node1,node2];inod=(1,2),field=(:tx1,:tx1),res=res)

SingleAcost <: AbstractElement

An element with a single node, for adding a once-off cost to a single A-dof.

Named arguments to the constructor

• field::Symbol.
• cost::Functor, where
  • cost(a::ℝ,[,costargs...]) → ℝ
• costargs::NTuple=() or NamedTuple of additional arguments splatted when calling cost.

Requestable internal variables

• cost, the value of the cost.

Example

using Muscade
model = Model(:TestModel)
node  = addnode!(model,𝕣[0,0])
@functor with() acost(a,three)=(a/three)^2
e     = addelement!(model,SingleAcost,[node];field=:EI,
                    costargs=(3.,),cost=acost)

Note that 3. is within a NamedTuple in the call to addelement!, but three is not, in the definition of the Functor acost.

See also: DofCost, SingleDofCost, Acost, ElementCostAndConstraint

SingleDofCost{Derivative,Class,Field,Tcost} <: AbstractElement

An element with a single node, for adding a cost to a given dof.

Named arguments to the constructor

• class::Symbol, either :X or :U.
• field::Symbol.
• cost::Functor, where
  • cost(x::ℝ,t::ℝ[,costargs...]) → ℝ if class is :X or :U, and
  • cost(x::ℝ, [,costargs...]) → ℝ if class is :A.
• costargs::NTuple=() or NamedTuple of additional arguments passed to cost`
• derivative::Int=0 0, 1 or 2 - which time derivative of the dof enters the cost.

Requestable internal variables

• cost, the value of the cost.

Example

using Muscade
model = Model(:TestModel)
node  = addnode!(model,𝕣[0,0])
@functor with() xcost(x,t,three)=(x/three)^2
e     = addelement!(model,SingleDofCost,[node];class=:X,field=:tx,
                    costargs=(3.,),cost=xcost)

Note that 3. is within a NamedTuple in the call to addelement!, but three is not, in the definition of the Functor xcost.

See also: DofCost, ElementCostAndConstraint

SingleUdof{XField,Ufield,Tcost} <: AbstractElement

An element that creates a Udof, and associates a cost to its value. The value of the Udof is applied as a load to a Xdof on the same node.

Named arguments to the constructor

• Xfield::Symbol.
• Ufield::Symbol.
• cost::Functor, where cost(u::ℝ,t::ℝ[,costargs...]) → ℝ
• costargs::NTuple=() or NamedTuple of additional arguments passed to cost.

Requestable internal variables

• cost, the value of the cost.

Example

using Muscade
model = Model(:TestModel)
node  = addnode!(model,𝕣[0,0])
@functor with() ucost(u,t,three)->(u/three)^2
e     = addelement!(model,SingleUdof,[node];Xfield=:tx,Ufield=:utx,
                    costargs=(3.,),cost=ucost)

Note that 3. is within a NamedTuple in the call to addelement!, but three is not, in the definition of the Functor ucost.

See also: DofCost, ElementCostAndConstraint

axis = Muscade.SpyAxis()

Spoof a GLMakie.jl Axis/Axis3 object so that calls like

lines!(  axis,args...;kwargs...)

result in args and kwargs being stored in axis, allowing to test functions that generate plots. Results are accessed by for example

axis.call[3].fun        
axis.call[3].args[2]

To get the name of the 3rd GLMakie.jl function that was called, and the 2nd input argument of this call.

Only lines!, scatter! and mesh! logging functions are implemented for now, but more functions can easily be added.

SweepX{OX}

A non-linear, time domain solver, that solves the problem time-step by time-step. Only the X-dofs of the model are solved for, while U-dofs and A-dofs are unchanged.

• SweepX{0} is Newton-Raphson.
• SweepX{1} is a first order variant of Newmark-β with Newton-Raphson iterations.
• SweepX{2} is Newmark-β, with Newton-Raphson iterations.

Important:

Muscade does not allow elements to have state variables, for example, plastic strain, or shear-free position for dry friction. Where the element implements such physics, this is implemented by introducing the state as a degree of freedom of the element, and solving for its evolution, even in a quasi-static problem, requires the use of OX≥1.

An analysis is carried out by a call with the following syntax:

initialstate    = initialize!(model)
setdof!(initialstate,1.;class=:U,field=:λcsr)
states           = solve(SweepX{2};initialstate=initialstate,time=0:10)

Named arguments to solve:

• dbg=(;) a named tuple to trace the call tree (for debugging)
• verbose=true set to false to suppress printed output (for testing)
• silenterror=false set to true to suppress print out of error (for testing)
• initialstate a State, obtain from ìnitialize! or SweepX, describing the initial conditions of the problem.
• time maximum number of Newton-Raphson iterations
• β=1/4,γ=1/2 parameters to the Newmark-β algorithm. β is dummy if OX<2. γ is dummy if OX<1.
• maxiter=50 maximum number of equilibrium iterations at each step.
• maxΔx=1e-5 convergence criteria: norm of X.
• maxLλ=∞ convergence criteria: norm of the residual.
• saveiter=false set to true so that output states contains the state at the iteration of the last step analysed. Useful to study a step that fails to converge.

Output

A vector of length equal to that of the named input argument time containing the states at the time steps.

See also: solve, initialize!, findlastassigned, study_singular, DirectXUA, FreqXU

SweepXA{OX}

A non-linear, time domain solver, that solves the problem time-step by time-step. Only the X-dofs of the model are solved for, while U-dofs and A-dofs are unchanged.

• SweepXA{0} is Newton-Raphson.
• SweepXA{1} is implicit Euler.
• SweepXA{2} is Newmark-β, with Newton-Raphson iterations.

IMPORTANT NOTE: Muscade does not allow elements to have state variables, for example, plastic strain, or shear-free position for dry friction. Where the element implements such physics, this is implemented by introducing the state as a degree of freedom of the element, and solving for its evolution, even in a quasi-static problem, requires the use of ORDER≥1.

An analysis is carried out by a call with the following syntax:

initialstate    = initialize!(model)
setdof!(initialstate,1.;class=:X,field=:tx1,order=1) 
states           = solve(SweepXA{2};initialstate=initialstate,time=0:10)

Named arguments to solve:

• dbg=(;) a named tuple to trace the call tree (for debugging)
• verbose=true set to false to suppress printed output (for testing)
• silenterror=false set to true to suppress print out of error (for testing)
• initialstate a State, obtain from ìnitialize! or SweepXA.
• time maximum number of Newton-Raphson iterations
• β=1/4,γ=1/2 parameters to the Newmark-β algorithm. Dummy if OX<2
• maxXiter=50 maximum number of equilibrium iterations at each step.
• maxΔx=1e-5 convergence criteria: norm of X.
• maxLλ=∞ convergence criteria: norm of the residual.
• maxAiter=50 maximum number of A-iterations
• maxΔa=1e-5 convergence criteria: norm of A.
• maxLa=∞ convergence criteria: norm of the residual.

Output

A vector of length equal to that of the named input argument time containing the states at the time steps.

See also: solve, initialize!, findlastassigned, study_singular, DirectXUA, FreqXU

y,... = apply{:chainrule}(f,x)
y,... = apply{:direct   }(f,x)

In the context of forward automatic differentiation using ∂ℝ, apply{:chainrule} accelerates the evaluation of y,...= f(x) if the length of x is smaller than the length of its partials.

apply{:direct} simply executes f(x) (no chain rule is applied)

Also works where x is a nested structure of Tuples and NamedTuples where the leaves are ℝ or SArray{S,R} where {S,R<:ℝ}.

Wrapper function of revariate and chainrule

chainrule_Jacobian{P}(Ty,X_)

Given Ty obtained using revariate, and X_, obtained using motion{P}(X) where X is a tuple of SVectors and P=precedence(X), compute y, a tuple of length ND of AbstractArrays of same eltype as vectors in X, andy∂X₀, the Jacobian of∂0(y)with respect to∂0(X)`.

See also revariate, motion, motion⁻¹

value      = default{:fieldname}(namedtuple,defaultvalue)
namedtuple = default(inputnamedtuple,defaultnamedtuple)

The first syntax attempts to access field fieldname from namedtuple. If namedtuple does not have such a field - or is not a NamedTuple, return defaultvalue.

The second syntax creates namedtuple from inputnamedtuple, supplementing with fields and values from defaultnamedtuple where there is no corresponding field in inputnamedtuple. This a thin wrapper of Julia's Base.merge.

Muscade.diffed_lagrangian{P}(eleobj;Λ,X,U,A,t=0.,SP=nothing)

Compute the Lagrangian, its gradients and Hessian, and the memory of an element. For element debugging and testing.

P, the order of differentiation must be 1 or 2.

The output is a NamedTuple with fields Λ, X, U, A, t, SP echoing the inputs and fields

• ∇L of format ∇L[iclass][ider]so that for example ∇L[2][3] contains the gradient of the Lagrangian wrt to the acceleration. iclass is 1,2,3 and 4 for Λ, X, U and A respectively.
• if P==2: HL of format HL[iclass,jclass][ider,jder]so that for example HL[1,2][1,3] contains the mass matrix.
• FB as returned by lagrangian

See also: diffed_residual, print_element_array

state = increment{OX}(initialstate,eiginc,imod,A)

Starting from initalstate for which an EigX analysis has been carried out, and using the output eiginc of that analysis, construct new States representing the instantaneous state of the vibrating structure

Input

• OX the number of time derivatives to be computed. increment(initialstate,eiginc,imod,A) defaults to OX=2
• initialstate the same initial State provided to EigX to compute eiginc
• eiginc obtained from EigX
• imod, an AbstractVector of integer mode numbers
• A, an AbstractVector of same length as imod, containing real or complex amplitudes associated to the modes

Output

• state a snapshot of the vibrating system

See also: EigX

state = increment{OX}(initialstate,eiginc,iω,imod,A)

Starting from initalstate for which an EigX analysis has been carried out, and using the output eiginc of that analysis, construct new States representing the instantaneous state of the vibrating structure

Input

• OX the number of time derivatives to be computed. increment(initialstate,eiginc,imod,A) defaults to OX=2
• initialstate the same initial State provided to EigXU to compute eiginc
• eiginc obtained from EigXU
• iω, the number of the frequency to consider. ω=iω*Δω where Δω is an input to EigXU.
• imod, an AbstractVector of integer mode numbers
• A, an AbstractVector of same length as imod, containing real or complex amplitudes associated to the modes

Output

• state a snapshot of the vibrating system

See also: EigXU

P,ND,X_  = motion(X)
Y_       = f(Y_)    
Y        = motion⁻¹{P,ND}(Y_)

motion⁻¹ transforms Y_ into a NTuple Y. The values P and `ND must be as returned by motion. Time derivatives can be accessed as

Y₀ = ∂0(Y)
Y₁ = ∂1(Y)
Y₂ = ∂2(Y)

One can also obtain the i-th time derivative directly as

Yᵢ = motion⁻¹{P,ND,i}(Y_)

See also motion

V_ = revariate{P}(V)

The variable V is a nested structure of NamedTuples, Tuples and SArrays of Reals (possibly: ∂ℝs).

V is stripped of its partials, an revariated to order P.

V_ = revariate(V)

revariates to the order precedence(V). TV has the same structure as V. One typical usage is with a Tuple:

a_,b_,c_ = revariate((a,b,c))

revariate, in conjunction with chainrule can be used to improve performance when the length of V is smaller than the length of its partials.

The output(s) of revariate must not be combined in computation with variables obtained by variate or revariate (including the inputs Λ, X, U, A and t) to lagrangian and residual (or variables computed from these). To prevent bugs, mark variables "touched" directly or indirectly by revariate's output with e.g. an underscore to keep track of the separation.

See [toolbox.EulerBeam3D] for an example of usage.

See also: chainrule, variate_indices, variate

y = value{P}(Y)

Extract the value of an automatic differentiation object, or SArray of such objects.

See also: precedence, Muscade.variate_, δ_, ∂, VALUE, value_∂

y,yₓ = value_∂{P,N}(Y)
y,y′ = value_∂{P  }(Y)

Get value and derivative in one operation.

See also: precedence, Muscade.variate_, Muscade.δ_, value, ∂, VALUE

V = variate{O}(v[;context,scale])

variate v to the O-th order.

V = variate(v[;context,scale])

variates v to the first order.

The variable v is a nested structure of NamedTuples, Tuples and SArrays of Reals (possibly: ∂ℝs). V has the same structure as v.

The output of variate must not be combined in computation with variables obtained by motion, revariate or a separate call to variate (including the inputs Λ, X, U, A and t) to lagrangian and residual (or variables computed from these).

If output from variate (or variables computed from this ouput), are combined in computation with variables obtained by revariate, motion or a separate call to variate (including the inputs Λ, X, U, A and t) to lagrangian and residual (or variables computed from these), then the later variables must be included in context. context can be a variable or a `tuple of variables.

V = variate(∂1{X};context=(Λ,X,U,t)) # A left out
W = V*t+∂0(X)                        # safe
Y = V*A[1]                           # not safe!!!

Not for use within elements: scale has the same structure as v except that SArrays and NTuples can be replaced by AllElement(scaling) to apply scaling to all elements in theSArrayorNTupleinv`.

See also: motion, revariate, chainrule, variate_indices

V = variate0{O}(v[;context,scale])

The same as variate, but while all derivatives are as from variate, all values are set to 0.

ℂ (\bbC)

an alias for abstract type Complex{<:Real}. For use in dispatching. ℂ1... ℂ4 are AbstractArrays of dimensions 1 to 4. ℂ11 is an AbstractVector of AbstractVector.

ℝ (\bbR)

an alias for abstract type Real. For use in dispatching. ℝ1... ℝ4 are AbstractArrays of dimensions 1 to 4. ℝ11 is an AbstractVector of AbstractVector.

ℤ (\bbZ)

an alias for abstract type Integer. For use in dispatching. ℤ1... ℤ4 are AbstractArrays of dimensions 1 to 4. ℤ11 is an AbstractVector of AbstractVector.

yₓ = ∂{P,N}(Y)

Extract the gradient of an automatic differentiation object. If Y is a SArray, the index of the partial derivative is appended to the indices of Y.

y′ = ∂{P}(Y)

Extract the derivative of an automatic differentiation object (or SArray of such), where the variation was created by the syntax Muscade.variate_{P}.

See also: precedence, Muscade.variate_, δ_, value, VALUE, value_∂

𝔹 (\bbB)

an alias for Bool. For use in dispatching. 𝔹1... 𝔹4 are AbstractArrays of dimensions 1 to 4.

𝕓 (\bbb)

an alias for Bool. For use in struct definitions. 𝕓1... 𝕓4 are Arrays of dimensions 1 to 4.

𝕣 (\bbr)

an alias for Float64. For use in struct definitions. 𝕣1... 𝕣4 are Arrays of dimensions 1 to 4. 𝕣11 is a Vector of Vector.

𝕫 (\bbz)

an alias for Int64. For use in struct definitions. 𝕫1... 𝕫4 are Arrays of dimensions 1 to 4. 𝕫11 is a Vector of Vector.

get(model::Model,dofID::dofID)

where dofID is found in the output of get(model,eleID) and get(model,nodID), returns a datastructure:

• class, dofclass :X, :U or :A
• field, field of the dof
• nodID identifying the nod of the model bearing this dof
• scale scaling factor
• eles: an array with, for each element sharing this dof
  • eleID identifying the element within the model
  • eletyp the datatype of the element
  • ielnod the element's node number

get(model::Model,eleID::EleID)

where eleID is returned by addelement!, returns a datastructure:

• eletyp: the datatype of the element
• dofs: an array with, for each node of the element
  • dofID identifying the dof within the model
  • nodID identifying the nod of the model bearing this dof
  • class, dofclass :X, :U or :A
  • field, field of the dof
  • scale scaling factor
• nods:
  • nodID identifying the nod of the model bearing this dof
  • coord the coordinate of the node

get(model::Model,nodID::nodID)

where nodID is returned by addnode!, returns a datastructure:

• coord the coordinate of the node
• eles: an array with, for each element connected to this node
  • eleID identifying the element within the model
  • eletyp the datatype of the element
• nods:
  • dofID identifying the dof of the model
  • class, dofclass :X, :U or :A
  • field, field of the dof
  • scale scaling factor

c = a∘₁b

Compute the single-dot product of two arrays, so that cᵢⱼ=Σₖ aᵢₖ bₖⱼ where i and j can be multiple indices.

See also: ⊗,∘₂

c = a∘₂b

Compute the double-dot product of two arrays, so that cᵢⱼ=Σₖₗ aᵢₖₗ bₖₗⱼ where i and j can be multiple indices.

See also: ∘₁,⊗

c = a⊗b

Compute the exterior product of two arrays, so that cᵢⱼ=aᵢ bⱼ where i and j can be multiple indices. A synonym is ∘₀.

See also: ∘₁,∘₂

Muscade.GUI(eiginc,initialstate;[draw_shadow=true],[shadow=...],[model=...])

Taking the output eiginc obtained from an EigXU, and the state initstate provided to EigXU, provide a GUI to explore the results.

Optional keyword arguements are

• draw_shadow whether to superimpose a drawing of initstate
• shadow a NamedTuple with any arguments to be passed to draw! initstate
• model a NamedTuple with any arguments to be passed to draw! the EigXU modes.

See also EigXU

Muscade.GUI(state,refstate=state[1];dim=3,kwargs...)

Taking state, a Vector of States output by various solvers, provide a GUI to explore the results.

This assumes that elements' drawing methods are writen for GLMakie.

The GUI allows to intereactively amplify responses (Λ,X,U and A-dofs) to make then easier to visualise. For X-dofs, it is the difference from the refstate (by default: state[1]) that is amplified.

Optional keyword arguements are

• dim, 2 or 3 depending on whether elements assume Axis or Axis3
• kwargs keywords argument, that will be passed to draw!

See also EigXU

Taylor(Ty,x₀,x)

Ty::∂ℝ has partials to arbitrary order evaluated at x₀. These partials define a Taylor expansion, which Taylor evaluates at value x

Taylor handles nested structures of Tuples and SVectors of ∂ℝ, applying the expansion to each element.

See also: chainrule, McLaurin, revariate, apply

VALUE(Y)

Completely strip Y of partial derivatives. Use only for debugging purpose.

See also: precedence, Muscade.variate_, δ_, value, ∂, value_∂

eleid = addelement!(model,ElType,nodid;kwargs...)

Add one or several elements to model, connecting them to the nodes specified by nodid.

If nodid is an AbstractVector: add a single element to the model. eleid is then a single element identifier.

If nodid is an AbstractMatrix: add multiple elements to the model. Each row of nodid identifies the node of a single element. eleid is then a vector of element identifiers.

For each element, addelement! will call eleobj = ElType(nodes;kwargs...) where nodes is a vector of nodes of the element.

See also: addnode!, describe, coord

nodid = addnode!(model,coord)

If coord is an AbstractVector of Real: add a single node to the model. Muscade does not prescribe what coordinate system to use. Muscade will handle to each element the coord of the nodes of the element, and the element constructor must be able to make sense of it. nodid is a node identifier, that is used as input to addelement!.

If coord is an AbstractMatrix, its rows are treated as vectors of coordinates. nodid is then a vector of node identifiers.

See also: addelement!, coord , describe

mut,opt = Muscade.allocate_drawing(axis,eleobjs;kwargs...)

Elements that are to be displayed in graphical output must implement a method for Muscade.allocate_drawing.

The method is to allocate opt, a NamedTuple of data that will not be mutated from frame to frame, but are usefull in Muscade.update_drawing or Muscade.display_drawing!.

The method is also to allocate mut, a NamedTuple of data that will be mutated from frame to frame. When implementing graphics with GLMakie.jl, the fields of mut must be exactly the updatable inputs provided to GLMakie.jl's drawing primitives: in Muscade.display_drawing!

lines!(axis,mut.x,mut.y)

is acceptable, but

lines!(axis,mut.x[:,s],mut.y[:,s])
lines!(axis,mut.a.x,mut.a.y)

are not.

The content of Arrays in opt and mut can be undef-ined.

Inputs are:

• axis the "canvas" to draw on, typicaly a GLMakie.jl Axis.
• eleobjs an AbstractVector of element objects, of length nel.
• kwargs a NamedTuple containing the keyword arguments provided by the user. See default.

See also: Muscade.update_drawing, Muscade.display_drawing!

chainrule(Ty,x)

Apply a chain rule in automatic differentiation. For example

Tx = revariate(x)
Ty = f(Tx)
y  = chainrule(Ty,x)

is faster than

y  = f(x)

if the length of x is smaller than the length of its partials.

See also: revariate, apply

c = coord(node)

Used by element constructors to obtain the coordinates of a vector of Nodes handed by Muscade to the constructor. c is accessed as

c[inod][icoord]

where inod is the element-node number and icoord an index into a vector of coordinates.

Note that c[inod] points at the same memory as nod[inod].coord: do not mutate c[inod]!

See also: addnode!, addelement!, describe, solve

describe(model,spec)

Print out information about model. spec can be

• an EleID to describe an element, additional optional keyword argument:

  • [depth=0] Set to value>0 to show the data in the element object up to depth.

• a DofID to describe a dof.

• a NodID to describe a node,

• :doftyp to obtain a list of doftypes,

• :dof to obtain a list of dofs or

• :eletyp for a list of element types.

  describe(state,[class=:all])

Provide a description of the dofs stored in state. class can be either :all, :Λ, :ΛX, :X, :U, :A or :scale.

If an element type requires a customized printout of its data by describe, implement str = Muscade.sdump(eleobj::MyElementType;maxdepth=depth). describe will add the initial indentation of 6 spaces to the begining of each line.

See also: get to obtain data instead of printing to screen.

Muscade.diffed_residual(eleobj;X,U,A,t=0.,SP=nothing)

Compute the residual, its gradients, and the memory of an element. For element debugging and testing.

The output is a NamedTuple with fields X, U, A, t, SP echoing the inputs and fields

• R containing the residual
• ∇R of format ∇R[iclass][ider]so that for example ∇R[2][3] contains the mass matrix. iclass is 2,3 and 4 for X, U and A respectively.
• FB as returned by residual

See also: diffed_lagrangian, print_element_array

Muscade.display_drawing!(axis,MyElement,mut,opt)

Elements that are to be displayed in graphical output must implement a method for Muscade.display_drawing!.

Inputs are:

• axis the "canvas" to draw on, typicaly a GLMakie.jl Axis.
• MyElement used for dispatching to the right method.
• mut is a NamedTuple, as output by Muscade.update_drawing. More specificaly, if implementing graphics with GLMakie.jl, mut has been
• opt is as returned by Muscade.allocate_drawing

When implementing graphics with GLMakie.jl, the fields of mut must be exactly the updatable inputs provided to GLMakie.jl's drawing primitives: in Muscade.display_drawing!

lines!(axis,mut.x,mut.y)

is acceptable, but

lines!(axis,mut.x[:,s],mut.y[:,s])

is not. The reason is that when doing graphics with GLMakie.jl, Muscade will recursively wrap each field of mut into an Observable before calling the elements' methods display_drawing!. This allows Muscade to update the graphics by just calling Muscade.update_drawing for each element.

See also: Muscade.allocate_drawing, Muscade.update_drawing

Muscade.doflist(::Type{E<:AbstractElement})

Elements must provide a method for Muscade.doflist.

The method must take the element type as only input, and return a NamedTuple with fieldnames inod,class and field. The tuple-fields are NTuples of the same length. For example

Muscade.doflist( ::Type{<:Turbine}) = (inod =(1   ,1   ,2        ,2        ),
                                       class=(:X  ,:X  ,:A       ,:A       ),
                                       field=(:tx1,:tx2,:Δseadrag,:Δskydrag))

In Λ, X, U and A handed by Muscade to residual or lagrangian, the dofs in the vectors will follow the order in the doflist. Element developers are free to number their dofs by node, by field, or in any other way.

See also: Muscade.lagrangian, Muscade.residual, Muscade.no_second_order

c = Muscade.dots(a,b,Val(N))

Compute the N-dot product of two arrays. N=0, N=1 and N=2 respectively correspond to ⊗, ∘₁ and ∘₂.

graphic = draw!(axis   ,state[,els];kwargs...)
          draw!(graphic,state[,els];kwargs...)

Plot all or part of a Model.

Currently, only GLMakie.jl is supported and tested, but Muscade is designed to allow application developers to chose other graphic system, including exporting data to Paraview. GLMakie.jl is thus not a dependency of Muscade, and must be installed and invoked (using) separately to run demos provided with Muscade.

Application developers can implement methods Muscade.allocate_drawing, Muscade.update_drawing and Muscade.display_drawing! to make their element "drawable".

axis a GLMakie.jl Axis, a Muscade.SpyAxis (for automated testing of graphic generation), and in the future a HDF5/VTK file handle for export of data to Paraview. state a single State. els specifies which elements to draw and can be either

• a vector of EleIDs (obtained from addelement!`), all corresponding to the same concrete element type
• a concrete element type (see describe).
• omitted: all the element of the model are drawn.

kwargs... is any additional key words arguments that will be passed to the draw method of each element, for example to specify colors, etc. See the elements' documentation.

When a plot of the Model is first generated, axis must be provided, and draw! returns graphic. graphic can then be provided for further calls to draw! to update the graphic.

See also: getdof, @request, @espy, addelement!, solve

equal(t) → :equal

A function which for any value t returns the symbol equal. Useful for specifying the keyword argument mode=equal in adding an element of type `DofConstraint to a Model.

See also: DofConstraint, ElementCostAndConstraint, off, positive

ilast = findlastassigned(state)

Find the index ilast of the element before the first non assigment element in a vector state.

In multistep analyses, solve returns a vector state of length equal to the number of steps requested by the user. If the analysis is aborted, solve still returns any available results at the begining of state, and the vector state[1:ilast] is fully assigned.

See also: solve

dofres = getdof(state;[class=:X],field=:somefield,nodID=[nodids...],[order=0])

Obtain the value of dofs of the same class and field, at various nodes and for various states.

If state is a vector, the output dofres has size (ndof,nstate). If state is a scalar, the output dofres has size (ndof,).

See also: getresult, addnode!, solve

getndof(model|Element)
getndof(model|Element,class)
getndof(model|Element,(class1,class2,[,...]))

where class can be any of :X, :U, :A: get the number of dofs of each specified dof-classes for the variable model or the type Element. If no class is specified getndof return the sum of the number of dofs of all classes.

See also: describe

eleres = getresult(state,req,els)

Obtain an array of nested NamedTuples and NTuples of element results. req is a request defined using @request. state a vector of States or a single State. els can be either

• a vector of EleIDs (obtained from addelement!) all corresponding to the same concrete element type
• a concrete element type (see describe).

If state is a vector, the output dofres has size (nele,nstate). If state is a scalar, the output dofres has size (nele).

See also: getdof, @request, @espy, addelement!, solve, describe

t3        = getsomedofs(X,3)
rotations = getsomedofs(X,SVector(4:6))

Used by elements' residual, lagrangian or update_drawing to obtain some degrees of freedom.

X and U are provided as NTuples of SVector (of Matrix in the case of update_drawing](@ref)). getsomedofs forms a NTuple containing the selected array components. Where e.g. X is a Matrix, the indexing is applied to the first dimension (corresponding to the dofs)

See also: ∂0,∂1,∂2,motion

Muscade.getδt(n,δω) = 2π/(n*δω)

`n` is the length of the time series

δω=Muscade.getδω(n,δt) = 2π/(n*δt)

`n` is the length of the time series

initialstate = initialize!(model)

Return an initial State for the model with all dofs set to zero

Modifying a model (invoquing addnode! and addelement! after initialize! will result in an error)

Optional keyword arguments: nΛder=1, nXder=1, nUder=1 to specify the number of "derivatives" to store in the State. note that "n⋅der==order+1", that is, for a dynamic problem (with accelerations), nXder=3 so that order==2. Setting n⋅der is only required if setdof! will be used. A dynamic solver handles a "static" initial state perfectly well. time=-∞ the time associated to the initial state. Note that for example setting time=0., and then calling a solver with a first time step also at 0. causes an error.

See also: setdof!, Model, addnode!, addelement!, solve

@espy function Muscade.lagrangian(eleobj::MyElement,Λ,X,U,A,t,SP,dbg)
    ...
    return L,FB
end

Elements must implement a method for Muscade.lagrangian or Muscade.residual.

Inputs

• eleobj an element object
• Λ a SVector{nXdof,R} where{R<:Real}, Lagrange multipliers (aka δX virtual displacements).
• X a NTuple of SVector{nXdof,R} where{R<:Real}, containing the Xdofs and, depending on the solver, their time derivatives. Use x=∂0(X), v=∂1(X) and a=∂2(X) to safely obtain vectors of zeros where the solver leaves time derivatives undefined.
• U a NTuple of SVector{nUdof,R} where{R<:Real}, containing the Udofs and, depending on the solver, their time derivatives. Use u=∂0(U), ̇u=∂1(U) and ̈u=∂2(U) to safely obtain vectors of zeros where the solver leaves time derivatives undefined.
• A a SVector{nAdof,R} where{R<:Real}.
• t a `Real containing the time.
• SP solver parameters (for example: the barrier parameter γ for interior point methods).
• dbg a NamedTuple to be used only for debugging purposes.

Outputs

• L the lagrangian
• FB feedback from the element to the solver (for example: can γ be reduced?). Return noFB of the element has no feedback to provide.

See also: Muscade.residual, Muscade.doflist, @espy, ∂0, ∂1, ∂2, noFB,

req = mergerequest(o.req)

"Wrapping" elements like ElementCostAndConstraint use requests to apply a cost or a constraint to requestables from another "target" element. These outer elements must be coded carefully so that getresult can be used to extracted both requestable internal results from the outer and from the target element.

mergerequest is used to merge the requests for the request needed to enforce a cost or constraint, and the user's request for element to be obtained from the analysis. The call to mergerequest, to be inserted in the code of lagrange for the outer element will be modified by @espy to something like req = mergerequest(o.req,req), to merge o.req of the outer element to any requests req transmitted by the user to extract results (or by an outer element to the outer element).

See the code of ElementCostAndConstraint's constructor and lagrange method for an example.

See also: ElementCostAndConstraint, @request, getresult

Muscade.mod_onebased(i,n) = mod(i-1,n)+1

For i::ℤ, returns a value in {1,...n}. This differs from mod which return a value in [0,n[

P,ND,X_ = motion(X [;context])

Input variables X and U provided toresidual or Lagrangian are NTuples which can have a form like (X₀,X₁,X₂), with the zeroth, first and second time derivative. motion transform such a NTuple into a SVector of ∂ℝ, containing derivatives with respect to time.

Using a scalar instead of SVector to provide an idea, (motion actually only operates on NTuple of SVectors) this transforms (x₀,x₁,x₂) into x₀+∂₁⟨x₁⟩+∂₂⟨x₁+∂₁⟨x₂⟩⟩.

This enables an element to use automatic differentiation to use time derivatives of various results, such as Coriolis and centrifugal accelerations, or strain rates.

If the third output from motion (or variables computed from this 3rd ouput), are combined in computation with variables obtained by variate or revariate (including the inputs Λ, X, U, A and t) to lagrangian and residual (or variables computed from these), then the variables obtained by variate or revariate must be included in context. context can be a variable or a `tuple of variables. To prevent bugs, mark variables "touched" directly or indirectly by motion's 3rd output with e.g. an underscore to keep track of the separation.

P,ND,X_ = motion(X ,context=(A,t)])
b       = f(A)         # b is "touched by A", which may be an addifed input to `residual`
Y_      = g(X_,b,t)    # Y_ is touched by t, and indirectly by A.  This is OK, they where 
                       # declared in the context of X_

See [toolbox.EulerBeam3D] for an example of usage.

!!! Warning Use ∂0, ∂1 and ∂2 to extract time derivatives from X and U, do not use indexing X[1] as the NTuple may be truncated in e.g. static analyses.

See motion⁻¹, variate

muscadeerror([[dbg,]msg])

Throw a MuscadeException, where

• dbg is a NamedTuple that contains "location information"

(for example: solver, step, iteration, element, quadrature point) that will be displayed with the error message.

• msg is a String describing the problem.

Muscade.no_second_order(::Type{E<:AbstractElement})

Elements that define residual are normaly mostly differentiated only to the first order, to avoid excessive compilation and/or execution time. To allow differentiation to the second order (for elements with few dofs), implement a method after the below pattern:

Muscade.no_second_order(::Type{<:MyElementType}) = Val(true)

Muscade.npartial(a::∂ℝ{P,N,R}) → N
Muscade.npartial(typeof(a)) → N

Also handle static arrays and tuples.

See also: precedence

off(t) → :off

A function which for any value t returns the symbol off. Useful for specifying the keyword argument mode=off in adding an element of type `DofConstraint to a Model.

See also: DofConstraint, ElementCostAndConstraint, equal, positive

Muscade.plot_block_matrix_sparsity(M)

Specialised tool to visualise the sparsity pattern of a matrix produced by DirectXUA. M is either a Matrix or a SparseMatrixCSC (the block structure), whose entries are themselve SparseMatrixCSC (the structure of each block).

Optional inputs:

• size=500 Size in pizel of the figure window.
• markersize=3 Size of dots for non-zero elements.

Muscade.plot_matrix_sparsity(M)

Opens a GLMakie figure and plots the sparsity pattern of M::SparseMatrixCSC.

Actual non zero-elements are plotted in green. Remaining structuraly non-zero elements are plotted in red.

Optional inputs:

• size=500 Size in pizel of the figure window.
• title=nothing Title of the figure window.
• markersize=3 Size of dots for non-zero elements.
• atol=1e-9 Tolerance for actual non-zero elements.

positive(t) → :positive

A function which for any value t returns the symbol positive. Useful for specifying the keyword argument mode=positive in adding an element of type `DofConstraint to a Model.

See also: DofConstraint, ElementCostAndConstraint, off, equal

precedence(a::∂ℝ{P,N,R}) → P
precedence(typeof(a)) → P

Also handles static arrays and tuples.

See also: npartial

Muscade.print_element_array(eleobj,class,V)

Show a vector (or a matrix) V, the rows of V being described as corresponding to eleobj dof of class class (:X, :U or :A). This can be used to print degrees of freedom, residuals, their derivatives, or gradients and Hessian of the Lagrangian.

See also: diffed_residual, diffed_lagrangian

Muscade.print_nz(S::SparseMatrixCSC)

List the structuraly non-zero entries of the sparse matrix.

@espy function Muscade.residual(eleobj::MyElement,X,U,A,t,SP,dbg) ... return R,FB end

Elements must implement a method for Muscade.residual or Muscade.lagrangian.

Inputs

• eleobj an element object
• X a NTuple of SVector{nXdof,R} where{R<:Real}, containing the Xdofs and, depending on the solver, their time derivatives. Use x=∂0(X), v=∂1(X) and a=∂2(X) to safely obtain vectors of zeros where the solver leaves time derivatives undefined.
• U a NTuple of SVector{nUdof,R} where{R<:Real}, containing the Udofs and, depending on the solver, their time derivatives. Use u=∂0(U), ̇u=∂1(U) and ̈u=∂2(U) to safely obtain vectors of zeros where the solver leaves time derivatives undefined.
• A a SVector{nAdof,R} where{R<:Real}.
• t a `Real containing the time.
• SP solver parameters (for example: the barrier parameter γ for interior point methods).
• dbg a NamedTuple to be used only for debugging purposes.

Outputs

• R the residual
• FB feedback from the element to the solver (for example: can γ be reduced?). Return noFB of the element has no feedback to provide.

See also: Muscade.lagrangian, Muscade.no_second_order, Muscade.doflist, @espy, ∂0, ∂1, ∂2, noFB

str = sdump(complex(2.,3.))

Same as dump, but returns a string

state = setdof!(state,value        ;[class=:X],field=:somefield,                  [order=0])
state = setdof!(state,value::Vector;[class=:X],field=:somefield,nodID=[nodids...],[order=0])

Set the value of dofs of the same class and field, at various nodes and for various states. There are two methods:

1. A single value is applied to all relevant nodes in the model
2. value and nodID are vectors of the same lengths, and each element in value is applied to the corresponding node.

setdof! is peculiar in that it modifies its input state variable, but must be used as a function. A call like stateout = setdof!(statein,value;class=:X,field=:somefield,order=1) can turn out in two ways: If the state already stores derivatives in X to order 1, then statein is mutated and statein===stateout. Otherwise, statein is unchanged, stateout is a new object, sharing as much memory as possible with statein. To avoid confusion, always use the syntax shown above.

See also: getresult, addnode!, solve

setscale!(model;scale=nothing,Λscale=nothing)

Provide scale value for each type (class and field) of dof in the model. This is usued to improve the conditioning of the incremental problems and for convergence criteria. scale is a NamedTuple with fieldnames within X, U and A. Each field is itself a NamedTuple with fieldnames being dof fields, and value being the expected order of magnitude.

For example scale = (X=(tx=10,rx=1),A=(drag=3.)) should be read as: X-dofs of field :tx are expected to be of the order of magnitude of 10m in the solution, :rx to be of the order of 1 radian, and A-dofs of field drag of the order of 3. All other degrees of freedom are of the order of 1.

Determining scaling coefficients that improve the condition number of incremental matrices is a hard problem.

Λscale is a scalar. The scale of a Λ-dof will be deemed to be the scale of the corresponding X-dof, times Λscale.

See also: addnode!, describe, coord, Muscade.study_scale

solve(Solver;dbg=(;),verbose=true,silenterror=false,kwargs...)

Execute an analysis using solver Solver (e.g. SweepX, DirectXUA...), and safeguard partial results in the case of error.

Named arguments

• dbg=(;) a named tuple to trace the call tree (for debugging)
• verbose=true set to false to suppress printed output (for testing)
• silenterror=false set to true to suppress print out of error (for testing)
• kwargs... Further arguments passed on to the method solve provided by the solver

See also: SweepX, DirectXUA, initialize!

scale = Muscade.study_scale(state;[SP=nothing],[verbose=false],[dbg=(;)])

Returns a named tuple of named tuples for scaling the model, accessed as scaled.myclass.myfield, for example scale.X.tx1.

If verbose=true, prints out a report of the analysis underlying the proposed scale. The proposed scaling depends on the state passed as input - as it is computed for a given incremental matrix.

See also: setscale!

matrix = Muscade.study_singular(state;SP,[iclasses=(Λ,:X,:U,:A)],[jclasses=iclasses],[verbose::𝕓=true],[dbg=(;)])

Generates an incremental matrix for state (no time derivatives) corresponding to the classes required, and report on the null space of the matrix.

In teh present implementation, the incremental matrix is converted to full format, limiting the applicability to small models.

The function returns the incremental matrix.

Muscade.toggle(condition,a,b)

Typestable equivalent of condition ? a : b. Returns a value converted to promote_type(typeof(a),typeof(b))

mut = Muscade.update_drawing(  axis,::AbstractVector{E},oldmut,opt, Λ,X,U,A,t,SP,dbg)

Elements that are to be displayed in graphical output must implement a method for Muscade.allocate_drawing.

For parametric element types

Muscade.update_drawing(axis,o::AbstractVector{Teleobj}, Λ,X,U,A,t,SP,dbg;kwargs...) 
    where{Teleobj<:MyElement}

For non-parametric element types, one can simplify the above to:

Muscade.update_drawing(axis,o::AbstractVector{MyElement}, Λ,X,U,A,t,SP,dbg;kwargs...)

Inputs are:

• axis the "canvas" to draw on, typicaly a GLMakie.jl Axis.
• eleobjs an AbstractVector of element objects, of length nel.
• oldmut the output mut of Muscade.allocate_drawing or of a previous call to Muscade.update_drawing.
• opt the output opt of Muscade.allocate_drawing
• Λ a matrix of size (nXdof,nel)
• X a NTuple (over the derivatives) of matrices of size (nXdof,nel)
• U a NTuple (over the derivatives) of matrices of size (nUdof,nel)
• A a matrix of size (nAdof,nel)
• t time
• SP solver parameters
• dbg debuging information
• kwargs a NamedTuple containing the keyword arguments provided by the user. See default.

See also: Muscade.allocate_drawing, Muscade.display_drawing!

position = ∂0(X)

Used by elements' residual or lagrangian to extract the zero-th order time derivative from the variables X and U.

See also: ∂1,∂2,getsomedofs

velocity = ∂1(X)

Used by elements' residual or lagrangian to extract the first order time derivative from the variables X and U. Where the solver does not provide this derivative (e.g. a static solver), the output is a vector of zeros.

See also: ∂0,∂2,getsomedofs

position = ∂2(X)

Used by elements' residual or lagrangian to extract the zero-th order time derivative from the variables X and U. Where the solver does not provide this derivative (e.g. a static solver), the output is a vector of zeros.

See also: ∂0,∂1,getsomedofs

X = Muscade.𝔉(x,δt)  # typeset with \mfrakF\Bbbr

Fourrier transform of a real time series x stored at time steps δt and length 2N = 2*2^p into a complex spectre X stored at frequency intervals δω=getδω(2N,δt)=2π/(2N*δt). The length of the spectre is N: only positive frequencies are stored (the Fourrier transform of real functions are Hermitian).

This provides a discretization of the unitary Fourrier transform,

G(ω) = 𝔉(g)(ω) = 1/√(2π) ∫exp(-𝑖ωt) g(t) dt

𝔉 is unitary, in the sense that

sum(abs2.(x))*δt ≈ 2*(sum(abs2.(X)) - abs2.(X[1])/2)*δω

(since the discrete spectre is provided for ω≥0, it contains only half the energy)

Arguments

• x a vector of real numbers representing a time series. Its length must be a power of two.
• δt the time step of the time series

Example

X   = 𝔉(x,δt) 
δω  = getδω(length(x),δt)
x′  = 𝔉⁻¹(X,δω) # ≈ x

See also: 𝔉⁻¹, getδω, getδt,

x = Muscade.𝔉⁻¹(X,δω)  # typeset with \mfrakF\^-\^1

See 𝔉

Arguments

• X a vector of complex numbers representing one side of a spectra. Its length must be a power of two.
• δω, the angular frequency step of spectra

Example

X   = 𝔉(x,δt) 
δω  = getδω(length(x),δt)
x′  = 𝔉⁻¹(X,δω) # ≈ x

See also: 𝔉⁻¹, getδω, getδt,

@espy function ... end

From an anotated function code, generate - "clean" code, in which the anotations have been deleted, and with the call syntax argout... = foo(argin...) - "espying" code, with added input and ouput arguments argout...,res = foo(argin...,req) where req has been generated using @request and res is a nested structure of NamedTuples and NTuples containing the requested data.

The macro is not general: it is designed for residual and lagrangian, which for performance have to be programmed in "immutable" style: they must never mutate variables (this implies in particular, no adding into an array in a loop over Gauss points). So @espy only supports the specific programming constructs needed in this style.

The following is an example of anotated code:

@espy function residual(x::Vector{R},y) where{R<:Real}
    ngp=2
    accum = ntuple(ngp) do igp
        ☼z = x[igp]+y[igp]
        ☼s,☼t  = ☼material(z)
        ♢square = s^2
        @named(s) 
    end
    r = sum(i->accum[i].s,ngp)
    return r,nothing,nothing
end

• The keyword function is preceded by the macro-call @espy.
• The name of requestable variables is preceded by ☼ (\sun). Such anotation must always appear on the left of an assigment.
• If the name of a variable is preceded by ♢ (\diamond), then the variable is evaluated only if requested. Such a notation can only be used if there is only one variable left of the assignement.
• The name of a function being called must be preceded by ☼ if the definition of the function is itself preceeded by the macro-call @espy.
• for-loops are not supported. do-loops must be used: to be efficient, residual and lagrangian must not allocate and thus use immutables.
• One-line function definition is not supported.
• The keyword return must be explicitly used, and if must be followed the a comma separated list of output variables. Syntaxes like return if... are not supported.

See also: @request, @espydbg, getresult

Muscade.@espydbg function ... end

Generate the same code as @espy and print it (for debug purposes).

See also: @espy, @request

a = 3
@functor with(a,e=2) function f(x::Real)
    return a*x^e
end

or

a = 3
@functor with(a,e=2)  f(x::Real)=a*x^e
e = 1
@functor with(a,e)    f(x::Real)=a*x^e
@functor with()       f(x::Real)=x^2

Creates a function-like object, of type Functor.

This is roughly equivalent to a closure defined as

f(x::Real)=a*x^e

Functors are meant to facilitate the definition of "functions" in a Muscade input script, while avoiding several of the issues associated with defining a function (and in particular a closure) in a script:

• A closure captures a variable "by reference", while @functor captures it by value, which might be more intuitive.
• To ensure type stability, the variables captured by a closure would have to be declared const - forbidding to update the input value in a script without restarting Julia.
• If the code of the function is not changed, the function is not parsed and compiled again, accelerating the re-analysis.

It is not possible to associate multiple methods to a functor.

See also: Functor

req = @request expr

Create a request of internal results wanted from a function. Considering the function presented as example for @espy, examples of possible syntax include

req       = @request gp(s,z,material(a,b))
req       = @request gp(s)
req       = @request gp(material(a))

The first expression can be read as follows: "In the function, there is a do loop over variable igp, and within this loop a call to a function material. Results s and z are wanted from within the loop, and results a and b from within material.

The corresponding datastructure containing the results for each element is a nesting of NTuples and NamedTuples, and can be accessed as out.gp[igp].material.a and so forth.

See also: @espy, @espydbg

inftyp,rettyp = Muscade.@typeof(foo(args...[;kwargs...]))

Determine the inferred type and the returned type of the output[s] returned by the relevant method-instance of foo.
Useful to study type-stability in `lagrangian`, `residual` and more.
This does not work on all operating systems, and should thus only be used for debugging.  
In tests, use `Test.@inferred`.

AxisymmetricBarCrossSection

Data structure containing the cross section material properties, for example to a Bar3D

Arguments to the constructor

• EA :: 𝕣 is the axial stiffness [N]
• μ :: 𝕣 is the mass per unit length [kg/m]

Optional argument to the constructor (all set to zero by default)

• w :: 𝕣 is the weight per unit length [N/m]
• g̃ :: Functor describes the gravity field divided by acceleration of gravity [-], function of time, set to SVector(0.,0.,1.) by default
• Caₜ :: 𝕣 is the tangential added mass per unit length [kg/m]
• Clₜ :: 𝕣 is the tangential linear damping coefficient per unit length [N/m/(m/s)]
• Cqₜ :: 𝕣 is the tangential quadratic damping coefficient per unit length [N/m/(m/s)^2], for example from drag
• Caₙ :: 𝕣 is the normal added mass per unit length [kg/m]
• Clₙ :: 𝕣 is the normal linear damping coefficient per unit length [N/m/(m/s)]
• Cqₙ :: 𝕣 is the normal quadratic damping coefficient per unit length [N/m/(m/s)^2]

Example

EA = 10.
L₀ =  2.
μ = 1. 
model           = Model(:TestModel)
node1           = addnode!(model,𝕣[0,0,0])
node2           = addnode!(model,𝕣[L₀,0,0])
mat             = AxisymmetricBarCrossSection(EA=EA,μ=μ)

See also: Bar3D, EulerBeam3D

Bar3D <: AbstractElement

A three-dimensional bar element that comes in two flavors: with two nodes and six X-dofs, or with an additional third node that carries three U-dofs.

Description of the degrees of freedom

The node coordinates provided to the constructor of Bar3D define the "as-meshed" coordinates of the nodes. The three X-dofs carried by each node, with field names :t1, :t2, and :t3, describe the displacement vector of that node with respect to the as-meshed coordinates in a cartesian global coordinate system.

Three additional U-dofs can be added by calling addelement! with Bar3D{true} instead of Bar3D. These U-dofs, carried by a third node, and with field names :t1, :t2, and :t3, represent the three components of an unknown uniformly distributed load (unit will be force per unit length) on the bar, expressed in the global coordinate system.

Arguments to the constructor

• nod :: Vector{Node} contains the element's nodes
• mat :: Mat contains the material properties (AxisymmetricBarCrossSection, for example)

Optional argument to the constructor

• ϵₛ ::𝕣 is such that the stress-free length of the element is (1-ϵₛ) times the as-meshed length of the element.

Providing ϵₛ is optional and set to machine precision by default. A non-zero ϵₛ means that the bar element exhibits some strain in the as-meshed configuration, and hence has some transverse stiffness, which facilitates convergence in static analyses.

Example

EA = 10.
L₀ =  2.
μ = 1. 
model           = Model(:TestModel)
node1           = addnode!(model,𝕣[0,0,0])
node2           = addnode!(model,𝕣[L₀,0,0])
mat             = AxisymmetricBarCrossSection(EA=EA,μ=μ)
addelement!(model,Bar3D,[node1 node2];mat)

See also: AxisymmetricBarCrossSection, EulerBeam3D

BeamCrossSection

Data structure containing the cross section material properties, for example to a EulerBeam3D

Arguments to the constructor

• EA :: 𝕣 is the axial stiffness [N]
• EI₂ :: 𝕣 is the bending stiffness [Nm/(1/m)] about second axis
• EI₃ :: 𝕣 is the bending stiffness [Nm/(1/m)] about third axis
• GJ :: 𝕣 is the torsional stiffness [Nm/(rad/m)] about longitudinal axis
• μ :: 𝕣 is the mass per unit length [kg/m]
• ι₁ :: 𝕣 is the (mass) moment of inertia about longitudial axis per unit length [kgm²/m]

Optional argument to the constructor (all set to zero by default)

• w :: 𝕣 is the weight per unit length [N/m]
• g̃ :: Functor describes the gravity field divided by acceleration of gravity [-], function of time, set to SVector(0.,0.,1.) by default
• Ca₁ :: 𝕣 is the tangential added mass per unit length [kg/m]
• Cl₁ :: 𝕣 is the tangential linear damping coefficient per unit length [N/m/(m/s)]
• Cq₁ :: 𝕣 is the tangential quadratic damping coefficient per unit length [N/m/(m/s)^2], for example from drag
• Ca₂ :: 𝕣 is the transverse added mass per unit length [kg/m] for motions along second axis
• Cl₂ :: 𝕣 is the transverse linear damping coefficient per unit length [N/m/(m/s)] for motions along second axis
• Cq₂ :: 𝕣 is the transverse quadratic damping coefficient per unit length [N/m/(m/s)^2], for motions along second axis
• Ca₃ :: 𝕣 is the transverse added mass per unit length [kg/am] for motions along third axis
• Cl₃ :: 𝕣 is the transverse linear damping coefficient per unit length [N/m/(m/s)] for motions along third axis
• Cq₃ :: 𝕣 is the transverse quadratic damping coefficient per unit length [N/m/(m/s)^2], for motions along third axis

Example

EA = 10.
EI₂ = 1.
EI₃ = 1.
GJ = 1.
μ = 1. 
ι₁ = 0.1
mat             = BeamCrossSection(EA=EA,EI₂=EI₂,EI₃=EI₃,GJ=GJ,μ=μ,ι₁=ι₁)

See also: EulerBeam3D, AxisymmetricBarCrossSection

EulerBeam3D <: AbstractElement

A three-dimensional Euler beam element, with two nodes and twelve X-dofs. An additonal third node carrying three U-dofs can optionally be added.

Description of the degrees of freedom

The node coordinates provided to the constructor of an EulerBeam3D define the "as-meshed" coordinates of the nodes. The three X-dofs carried by each node with field names :t1, :t2, and :t3 describe the displacement vector of that node with respect to the as-meshed coordinates in a cartesian global coordinate system. The three X-dofs carried by each node with field names :r1, :r2, and :r3, describe the rotation of this node, using a Rodrigues representation: the direction of (:r1,:r2,:r3) defines the axis of rotation, and the magnitude of this vector defines the angle of rotation (in gradian).

Three additional U-dofs can be added by calling addelement! with EulerBeam3D{true} instead of EulerBeam3D. These U-dofs, carried by a third node, and with field names :t1, :t2, and :t3, represent the three components of an unknown uniformly distributed load (unit will be force per unit length) on the beam, expressed in the global coordinate system.

Arguments to the constructor

• nod :: Vector{Node} contains the element's nodes
• mat :: Mat contains the material properties (BeamCrossSection, for example)

Optional argument to the constructor

• orient2 :: SVector{3,𝕣} defines the direction of the first bending axis in the global coordinate system.

Default is SVector(0.,1.,0.).

Example

EA = 10.
EI₂ = 3.
EI₃ = 3.
GJ = 4.
μ = 1.
ι₁ = 1.0
L = 5
model = Model(:TestModel)
node1 = addnode!(model,𝕣[0,0,0])
node2 = addnode!(model,𝕣[L,0,0])
mat = BeamCrossSection(EA=EA,EI₂=EI₂,EI₃=EI₃,GJ=GJ,μ=μ,ι₁=ι₁)
addelement!(model,EulerBeam3D,[node1 node2];mat)

See also: BeamCrossSection, Bar3D

EulerBeam3DwithStrainGauge

An element designed to wrap around an EulerBeam3D. The element makes no contribution to residual or Lagrangian besides the contribution made by the element it wraps. EulerBeam3DwithStrainGauge provides requestables allowing to model strain gauges. The strain gauges are placed halfway along the EulerBeam3D. In an inverse analysis, the EulerBeam3DwithStrainGauge is itself wraped by an ElementCostAndConstraint element.

Keyword arguments when adding elements:

• P SMatrix{3,Nsensor,𝕣}, giving the offset between the point on the axis of the element halfway along its length. P[1,:] must be zero.
• D SMatrix{3,Nsensor,𝕣}, the orientation of a strain gauge. If multiple strain gauges are present at the same position, one can repeat the columns of P.
• TargetElement=EulerBeam3D the constructor to the wrapped element. This is typicaly an EulerBeam3D but could be another element wrapping an EulerBeam3D.
• elementkwargs a NamedTuple with the keyword arguments to the wrapped element. See EulerBeam3D

Requestables:

• εₐₓ (scalar), the axial strain at the middle of the element
• κ a vector with elements: κ[1] the torsion at the middle of the element. κ[2] the component of te curvature in EulerBeam3D's direction 2. κ[3] the component of te curvature in EulerBeam3D's direction 3. Although a tri-vector that includes torsion may erroneously suggest a rotation rate vector, the κ[2:3] points to the inside of the curvature and the unit is in 1/L (not rad/L) where L is the unit of length of the model.
• ε a vector of length Nsensor containing the strain values defined by P and D.

#Keyword arguments when drawing elements:

When calling draw!, instruction to EulerBeam3DwithStrainGauge elements are given as in this example:

draw!(axis,state>;EulerBeam3DwithStrainGauge=(L= @SVector [0.1,0,0],point_size=10,accelerometer_color=:orange))`

The optional keword arguments and their default values are

• gauge_color = :blue

• expand = 1.02 a multiplicative factor applied to P in the drawing only (not in strain calculation) to ensure the gauge is not hidden by a patch-drawing of the beam itself.

• L = SVector{Nsensor,𝕣}(norm(P)/5 for i=1:Nsensor) the length of the strain gauge in the drawing TODO

ExcentricRigidConnection <: AbstractElement

An element to constrain a rigid body motion between a master node with fields (:t1,t2,:r3) and an emissary node with fields (:t1,:t2)

Example

using Muscade
model = Model(:TestModel)
MasterNode  =   addnode!(model,𝕣[0.,    0.,     0.])
EmissaryNode  = addnode!(model,𝕣[1/√2,  1/√2,   0.])
e     = addelement!(model,ExcentricRigidConnection,[MasterNode EmissaryNode])

See also: Hold, DofConstraint

Position3D

An 3D single-node element, to be connected to a node with both translation and rotation dofs. The element makes zero contribution to residual or Lagrangian. It only provides requestables allowing to model accelerometers and optical position measurements. For an inverse analysis, this is done by including the element in an ElementCostAndConstraint.

Keyword arguments when adding elements:

• P SMatrix{3,Nsensor,𝕣}, giving the offset between the nodal position and the point(s) at which position(s) and/or accelerations will be measured.
• D SMatrix{3,Nsensor,𝕣}, the orientation of an accelerometer. If no accelerometer is present at a given postion P, use NaNs. If multiple accelerometers are present at the same position, one can repeat the columns of P.

Requestables:

• a a[isensor] is the XXXXXXXXXXXX TODO
• x x[oder+1][:,isensor] contains the position, velocity and acceleration (oder=0,1,2) of the sensor which position was described in the isensor-th column of P.
• rᵢ rₑ[oder+1] is a vector containing a zero vector, the intrinsic rotation rate vector and its time derivative (oder=0,1,2), for the element's node.
• rₑ rₑ[oder+1] is a vector containing the rotation vector, the extrinsic rotation rate vector and its time derivative (oder=0,1,2), for the element's node.
• R R[oder+1] is a matrix containing the rotation matrix, spin matrix and its time derivative (oder=0,1,2), for the element's node.
• xₙ x[oder+1] is a vector containing the position, velocity and acceleration (oder=0,1,2) of the element's node.

#Keyword arguments when drawing elements:

When calling draw!, instruction to Position3D elements are given as in this example:

draw!(axis,state>;Position3D=(L= @SVector [0.1,0,0],point_size=10,accelerometer_color=:orange))`

The opional keword arguments and their default values are

• L = 0. (thus if not given, directions of accelerometers are not shown)
• point_size = 6
• point_color = :black
• stalk_color = :grey
• stalk_width = 1
• accelerometer_color = :teal
• accelerometer_width = 2

MeshLine!(model,topNode, azimuth, eltype, xSection, segLength, nel)

Create the mesh of a multi-segment line in a Muscade model using beam or bar elements. Nodes and elements are generated along a straight line, contained in the (:t1,:t2) plane, and oriented based on the given azimuth angle in radian. Note that the numbering of added nodes and elements starts from the bottom node.

Arguments

• model: The Muscade model
• topNode: The top node of the meshed line.
• azimuth: The azimuth angle (in radians) defining the direction of the line.
• eltype: The type of element to use for the mesh, for example Bar3D or EulerBeam3D.
• xSection: A vector of cross-section properties for each segment, for example BeamCrossSection or AxisymmetricBarCrossSection
• segLength: A vector of lengths for each segment.
• nel: A vector specifying the number of elements in each segment.

Keyword Arguments

• kwargsToElement...: Additional keyword arguments passed to the element constructor (e.g., orientation parameters for beams).

Returns

• nodeList: A vector of vectors containing the node IDs for each segment.
• elementList: A vector of element IDs for the entire mesh.
• nodeCoord: A vector of matrices containing the coordinates of nodes for each segment.

Notes

Assumes at least 2 segments.

Toolbox.Rodrigues(v::SVector{3})

Transform a rotation vector v into the rotation matrix M.

See also Toolbox.spin, Toolbox.spin⁻¹, Toolbox.Rodrigues⁻¹, Toolbox.adjust.

Toolbox.Rodrigues⁻¹(v::SVector{3})

Transform a rotation matrix M into the rotation vector v, such that |v| < π. Undefined for rotations of angle π

See also Toolbox.spin, Toolbox.spin⁻¹, Toolbox.Rodrigues, Toolbox.adjust.

Toolbox.adjust(u::SVector{3},v::SVector{3})

Compute the matrix of the rotation with smallest angle that transforms u into a vector colinear with v. Fails if |u|=0, |v|=0 or if the angle of the rotation is π.

See also Toolbox.spin, Toolbox.spin⁻¹, Toolbox.Rodrigues, Toolbox.Rodrigues⁻¹.

clutch(x,x₁,x₂,y₁,y₂,γ)

Interpolation between x₁ and x₂ in ℝ, to (possibly vector-valued) y₁ and y₂. The parameter γ enables controlling how fast clutch goes from y₁ and y₂. A typical application of clutch are static structural anlayses, where different load components should be applied with a specific sequence/magnitude to ensure convergence of the analysis

Inputs

• x ∈ ℝ is the point at which clutch should be evaluated
• x₁ ∈ ℝ is the point at which clutch starts. The output is equal to y₁ for x<x₁.
• x₂ ∈ ℝ is the point at which clutch stops. The output is equal to y₂ for x>x₂.
• y₁ and y₂ are the values of clutch for x<x₁ and x>x₂, respectively.
• γ>0 controls the shape of the interpolation. If γ=1, linear interpolation between x₁ and x₂ If γ<1, faster progression towards y₂ first, and then slower If γ>1, slower progression towards y₂ first, and then faster

Output

• The interpolated value at x, possibly vector-valued.

Example

x = -1:0.01:2
γ_ = [0.1,0.5,1.,2,10.]
using GLMakie
fig = Figure(size = (1000,1000))
ax = Axis(fig[1, 1],ylabel="clutch() value")
[lines!(ax,x,clutch(x,0,1,0,1,γ),label="γ="*string(γ)) for γ∈γ_]; 
axislegend(); 
xlims!(ax,-1,2); display(fig)

Toolbox.intrinsicrotationrates(rₑ::NTuple{ND,SMatrix{3,3}}) where{ND}

Transform a NTuple containing a rotation matrix and its extrinsic time derivatives, into a NTuple containing a (zero) rotation vector and its intrinsic time derivatives.

See also Toolbox.spin, Toolbox.spin⁻¹, Toolbox.Rodrigues, Toolbox.Rodrigues⁻¹.

Toolbox.scac(x)

scac(x) = sinc1(acos(x)), The function can be differentiated to the fourth order over ]-1,1] .

See also Toolbox.sinc1

Toolbox.spin(v::SVector{3})

Transform a rotation vector v into the cross product matrix M, such that M ∘₁ a = v × a.

See also Toolbox.spin⁻¹, Toolbox.Rodrigues, Toolbox.Rodrigues⁻¹.

toolbox.spin⁻¹(M::SMatrix{3,3})

Transform a cross product matrix M into the rotation vector v, such that v × a = M ∘₁ a.

See also Toolbox.spin, Toolbox.Rodrigues, Toolbox.Rodrigues⁻¹.