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

ϵ (\epsilon)

an alias for Base.eps(𝕣).

∞ (\infty)

an alias for Base.inf.

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 passed to 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)

See also: SingleAcost, DofCost, SingleDofCost, ElementCost, 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
• 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,Nx,Nu,Na,xinod,xfield,uinod,ufield,ainod,
    afield,λinod,λfield,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 the gap g≥0 and the 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...).

Named arguments to the constructor

• 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 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.

Example

using Muscade
model           = Model(:TestModel)
n1              = addnode!(model,𝕣[0]) 
@functor with() gap(x,t)=x[1]+.1
@functor with() res(x,u,a,t)=0.4x.+.08+.5x.^2)
e1              = addelement!(model,DofConstraint,[n1],xinod=(1,),xfield=(:t1,),
                              λinod=1, λclass=:X, λfield=:λ1,gap=gap,
                              mode=positive)
e2              = addelement!(model,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[1]

See also: Hold, ElementConstraint, 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 passed to 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)

See also: Acost, SingleDofCost, ElementCost, 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::ℝ) → ℝ.

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)=a*t+b
e     = addelement!(model,DofLoad,[node];field=:tx,value=load)

See also: Hold, DofCost

eiginc = solve(EigX{ℝ};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
• 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, geneig

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.
• 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

ElementConstraint{Teleobj,λinod,λfield,Nu,Treq,Tg,Tgargs,Tmode} <: AbstractElement

An element to apply optimisation equality/inequality constraints on the element-results of another "target" element. The target element must not be added separately to the model. Instead, the ElementType, and the named arguments to the target element are provided as input to the ElementConstraint constructor.

This element generates a time varying optimisation constraint. For example: find A-parameters so that at all times, the element-result von-Mises stress does not exceed a given value.

The Lagrangian multiplier introduced by this optimisation constraint is of class :U.

Named arguments to the constructor

• λinod::𝕫 The element-node number of the Lagrange multiplier.

• λfield::Symbol The field of the Lagrange multiplier.

• 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.

• gₛ::𝕣=1. A scale for the gap.

• λₛ::𝕣=1. A scale for the Lagrange multiplier.

• gap A gap function gap(eleres,t,gargs...)→ℝ. eleres is the output of the above-mentionned request to the target element.

• gargs::NTuple 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.

• ElementType The named of the constructor for the relevant element

• elementkwargs A named tuple containing the named arguments of the ElementType constructor.

Requestable internal variables

Not to be confused with the req provided as input to addelement! when adding an ElementConstraint, one can, after the analysis, request results from ElementConstraint

• λ The constraints Lagrange multiplier
• gap The constraints gap Functor
• 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 ElementConstraint. For example @request gap would extract the value of the ElementConstraint's function gap, while @request eleres(gap) refers to the value of a variable called gap in the target element.

Example

@functor with() gap(eleres,t) = eleres.Fh^2
ele1 = addelement!(model,ElementCoonstraint,[nod1];req=@request(Fh),
                   gap,λinod=1,λfield=:λ,mode=equal, 
                   ElementType=AnchorLine,
                   elementkwargs=(Δxₘtop=[5.,0,0], xₘbot=[250.,0],L=290., buoyancy=-5e3))

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

ElementCost{Teleobj,Treq,Tcost,Tcostargs} <: AbstractElement

An element to apply costs on another "target" element's dofs and element-results. The target element must not be added separately to the model. Instead, the ElementType, and the named arguments to the target element are provided as input to the ElementCost constructor.

Named arguments to the constructor

• 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.
• 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 passed to cost.
• ElementType The named of the constructor for the relevant element.
• elementkwargs A named tuple containing the named arguments of the ElementType constructor.

Example

@functor with() cost(eleres,X,U,A,t;Fh0) = (eleres.Fh-Fh0)^2
ele1 = addelement!(model,ElementCost,[nod1];req=@request(Fh),
                   cost=cost,
                   costargs = (;Fh0=0.27),
                   ElementType=AnchorLine,
                   elementkwargs=(Λₘtop=[5.,0,0], xₘbot=[250.,0], L=290., buoyancy=-5e3))

See also: SingleDofCost, DofCost, @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.
• 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  = FunctionFromVector(xs::AbstractRange,ys::AbstractVector)
y  = f(x)

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

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

An "identity vector"

id = IdVec i = 9834987 id[i] == i # true for any i

See also Julia's identity function.

model = Model([ID=:my_model])

Construct a blank model, which will be mutated to create a finite element [optimization] 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

• ElementType 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 ElementType 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

A simple offset vector that only implements "setindex!" and "getindex". Based on Memory, so cannot be extended.

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,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 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() acost(a,three)=(a/three)^2
e     = addelement!(model,SingleAcost,[node];field=:EI,
                    costargs=(3.,),cost=acost)

See also: DofCost, SingleDofCost, Acost, ElementCost

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 tocost``
• 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)

See also: DofCost, ElementCost

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 tocost``

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)

See also: DofCost, ElementCost

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 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 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.
• 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 unchangeδ

• SweepXA{0} is Newton-Raphson, with feasibility line-search, to handle inequality constraints.
• SweepXA{1} is implicit Euler, with feasibility line-search.
• SweepXA{2} is Newmark-β, with Newton-Raphson iterations and feasibility line search

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=:U,field=:λcsr)
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.
• maxLineIter=50 Maximum number of iteration in the feasibility line search. set to 0 to skip the line search (not recommended for models with inequality constraints).
• sfac=0.5 Parameter in the line search for a feasible point. If a tentative result is not feasible, backtrack by a factor sfac. If still not feasible, backtrack what is left by a factor sfac, and so forth, up to maxLineIter times.
• γfac=0.5 Parameter for feasibility. For an inequality constraint g(X) with reaction force λ, require g(X)*λ==γ, and multiply γ *= γfac at each iteration.

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

addin!(asm,global,block,ibr,ibc,factor=1.)

Add a sparse block into a large out sparse matrix, at block-row and -column ibr and ibc. Use prepare to allocate memory for global and build the assembler asm.

addin!(asm,outvec,blockvec,ibr)

Add a full block vector into a large outvec full vector. at block-row ibr. Use prepare to create asm.

See also: prepare

composeJacobian{P}(Ty,X_)

Given Ty obtained using revariate, and X_, obtained using motion{P}(X) where X is a tuple of SVectors and P=constants(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⁻¹, composevalue

composevalue{P,ND}(Ty,X_)

Given Ty obtained using revariate, and X_, obtained using motion{P}(X) where X is a tuple of length ND and P=constants(X), compute y, a tuple of length ND of AbstractArrays of same eltype as vectors in `X.

See also revariate, motion, motion⁻¹, composeJacobian

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

λ,v,ncv = geneig{ALGO}(A,B,neig=5)

Solves (A-λ*B)*v=0, finding the neig lowest eigenvalues λ (in absolute value) and the corresponding eigenvectors v (a Vector{Vector})

Input

ALGO can be

• :SDP if A is symmetric definite positive and B is symmetric. Will return real λ and v.
• :Hermitian if A is symmetric indefinite and B is symmetric. Will return real λ and v.
• :Complex otherwise, will return complex λ and v.

Optional keyword arguments:

• maxiter = 300
• verbosity = 0 ∈ {0,1,2,3}
• krylovdim = 2neig+6

Uses KrylovKit.jl. Freely based on VibrationGEPHelpers.jl and input from PetrKryslUCSD and stevengj. See GIThub-blame for bug-credits.

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  = constant(X)
X_ = motion{P}(X)

Transform a NTuple of SVectors, for example the vector X provided as an input to residual or Lagrangian into a SVector of ∂ℝ. This can be used by an element to compute time derivatives, for example Euler, Coriolis and centrifugal accelerations, or strain rates.

Some principles of safe automatic differentiation must be adhered to:

• the function that uses motion must also 'unpack' : no variable that is touched by the output of motion must be returned by the function without having been unpacked by motion⁻¹. Touched variables can for example be marked with an underscore.
• The precendence P must be calculated using constants with all variables that are input to the function and may be differentiated.
• If other levels of automatic differentiation are introduced within the function, unpack in reverse order of packing.

See motion⁻¹

P  = constants(X,U,A,t)
ND = length(X)
X_  = motion{P,ND}(X)
Y_  = f(Y_)    
Y₀ = motion⁻¹{P,ND,0}(Y_)
Y₁ = motion⁻¹{P,ND,1}(Y_)
Y₂ = motion⁻¹{P,ND,2}(Y_)
Y  = motion⁻¹{P,ND  }(Y)

Extract the value and time derivatives from a variable that is a function of the output of motion. In the above Y is a tuple of length ND. One can use ∂0,∂1 and ∂2 to unpack Y.

See also motion

TX = revariate{P}(V)

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

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

TV = revariate(VX)

revariates to the order precedence(V).

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

Be extremely careful never to mix any variable that is a function of V with any other variables containing ∂ℝs but not produced by the same revariate.

A special version of revariate

V = (;X,U,A)
S = (X=scale.X,U=scale.U,A=scale.A)
TV = revariate{P}(V,S)

allows to introduce scaling in automatic differentiation. For this method, S and V have the same structure, with the important exception that Tuples in V must be ntuples (have elements of identical type T), and, to a Tuple in V corresponds a variable of type T in V. Put simply: the same scale scale.X will be applied to ∂0[X], ∂1[X] and ∂2[X].

See also: chainrule

to_order{P}(V)

Decrease (lossy) or increase (pad partials with zeros) the order of differentiation of V. V is a nested structure of NamedTuple, Tuple, SArray, and the components of V must be of type ∂ℝ (otherwise, the number of partials would be undefined).

y = value{P}(Y)

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

See also: constants, variate, δ, ∂, VALUE, value_∂

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

Get value and derivative in one operation.

See also: constants, variate, δ, value, ∂, VALUE

X = variate{P,N}(x)

where typeof(x)<:SVector{N}, create a SVector of automatic differentiation objects of precedence P.

X = variate{P}(x)

where typeof(x)<:Real, create an object of precedence P.

See also: constants, δ, value, ∂, VALUE, value_∂

X = δ{P,N,R}()

create a SVector of automatic differentiation objects of precedence P and value zero.

X = δ{P}()

Create automatic differentiation object of precedence P and value zero.

See also: constants, variate, value, ∂, VALUE, value_∂

ℂ (\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.

ℕ (\bbN)

an alias for UInt64. For use in dispatching. ℕ1... ℕ4 are AbstractArrays of dimensions 1 to 4.

ℝ (\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 variate{P}.

See also: constants, 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.

𝕔 (\bbc)

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

𝕟 (\bbn)

an alias for UInt64. 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.

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.

See also: ∘₁,∘₂

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

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

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, as if Ty had been computed at 0.

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, fast

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, fast

@show VALUE(Y)

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

See also: constants, 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, fast

colnormalize(a)

Euclidian-normalize the columns of an SMatrix

columnmatrix(v)

Reshape a vector into a matrix of size (length(v),1)

P = constants(a,b,c)

Generate a precedence P that is higher than the precedence of the arguments.

See also: variate, δ, value, ∂, VALUE, value_∂

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,

• 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.

See also: addelement!, addnode!

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

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 eletyp).
• 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

et = eletyp(model)

Return a vector of the concrete types of elements in the model.

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, ElementConstraint, off, positive

y,... = fast(f,x)

In the context of forward automatic differentiation using ∂ℝ, accelerate the evaluation of y,...= f(x) if the length of x is smaller than the length of its partials.

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

Be extremely careful with closures, making sure that f does not capture variables of type ∂ℝ.

Wrapper function of revariate and McLaurin

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 asum 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 eletyp).

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, eletyp

rotations = getsomedofs(X,[3,6])

Used by elements' residual or lagrangian to some degrees of freedom, and their time derivatives, from the variables X and U.

See also: ∂0,∂1,∂2

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

`n` is the length of the time series

δω=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)

"Outer" elements like ElementCost and ElementConstraint 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 ElementCost's constructor and lagrange method for an example.

See also: ElementCost, @request, getresult

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[

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.

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(false)

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, ElementConstraint, 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, ElementConstraint, off, equal

bigmat,bigmatasm,bigvecasm,bigvecdis = prepare(pattern)

Prepare for the assembly of sparse blocks into a large sparse matrix. bigmat is allocated, with the correct sparsity structure, but its nzval undef'ed. Where some blocks share the same sparsity structure, blocks in pattern can have === elements.

pattern is a SparseMatrixCSC{<:SparseMatrixCSC}, where empty blocks are structuraly zero

See also: addin!

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

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

rowmatrix(v)

Reshape a vector into a matrix of size (1,length(v))

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!

sparser!(S::SparseMatrixCSC,keep::Function)

Eliminate terms that do not satisfy a criteria from the storage of a sparse matrix. S will be mutated, and the size of its internal storage modified. keep is a Function which to an index into S.nzval associate true if storage is to be kept for this term and false otherwise. Alternatively, keep can be a Vector{Bool}

Examples

sparser!([T],S,i->abs(S.nzval[i])>tol)
sparser!([T],S,keep::Vector{Boolean})

If T is provided (initialised as T = copy(S)), then result is set in T, S is unchanged, other wise S is mutated in place. The input keep::Vector{Boolean} must be of length nnz(S).

sparser!([S1,S2,...],rtol=1e-9)

Operates in place, reducing the sparses S1, S2 etc... to a common sparsity pattern.

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.

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!

zero!(a)

Set to zero all elements of an arrays. If a is sparse, the vector nzval of values is set to zero and the sparsity structure is unchanged.

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 = 𝔉(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 = 𝔉⁻¹(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,

𝕫log2(i::𝕫)

Compute the integer log2 of an integer, fast. Fails if i is not a power of two.

@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.
• 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

@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

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.

Functor is a subtype of the abstract type Function: functions that accept a arg::Function as an input will accept arg to be a Functor. Functions that require arg:Functor will not accept a classical Function, thus enforcing capture by value etc.

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)

• 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])
elnod           = [model.nod[n.inod] for n∈[node1,node2]]
mat             = AxisymmetricBarCrossSection(EA=EA,μ=μ)
bar             = Bar3D(elnod;mat)

See also: Bar3D, EulerBeam3D

Bar3D <: AbstractElement

A three-dimensional bar element, with two nodes, six X-dofs and three U-dofs

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])
elnod           = [model.nod[n.inod] for n∈[node1,node2]]
mat             = AxisymmetricBarCrossSection(EA=EA,μ=μ)
bar             = Bar3D(elnod;mat)

See also: AxisymmetricBarCrossSection, EulerBeam3D

EulerBeam3D

An Euler beam element

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 ElementCost.

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
• 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

StrainGaugeOnEulerBeam3D

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. StrainGaugeOnEulerBeam3D provides requestables allowing to model strain gauges. The strain gauges are placed halfway along the EulerBeam3D. In an inverse analysis, the StrainGaugeOnEulerBeam3D is itself wraped by an ElementCost 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.
• ElementType=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 StrainGaugeOnEulerBeam3D elements are given as in this example:

draw!(axis,state>;StrainGaugeOnEulerBeam3D=(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

Rodrigues(v::SVector{3})

Transform a rotation vector v into the rotation matrix M.

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

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.

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⁻¹.

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⁻¹.

scac(x)

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

See also Toolbox.sinc1

sinc1(x)

sinc1(x) = sin(x)/x - but sinc1(0.) = 1.. The function can be differentiated to the fourth order.

This differs from Julia's sinc(x) = sin(π*x)/(π*x).

See also Toolbox.scac

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⁻¹.

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⁻¹.

trace(v::SMatrix{3,3})

Computes the trace of a matrix.

Drawing a Bar3D.

draw!(axis,state)

Optional arguments (and their default values) are

• line_color = :black color of the line
• Udof (true iff element has Udofs) wether to draw U-forces.
• Uscale = 1. How many meter is a Newton per meter?

Drawing a EulerBeam3D.

draw!(axis,state)

draw!(axis,state;EulerBeam3D=(;style=:shape))

α      = 2π*(0:19)/20
circle = 0.1*[cos.(α) sin.(α)]'
draw!(axis,state;EulerBeam3D=(;style=:solid,section = circle))

style=:shape shows the deformed neutral axis of the element. It has optional arguments frame=true (draws the element's corotated frame of reference) and nseg=10 (number of points to show the deflected shape of each element).

style=:solid shows the deformed shape of the element. It requires the input section=... to be given a matrix of size (2,nsec) describing nsec points around the cross section of the element (no need to close the circumference by repeating the first point at the end). It has optional arguments nseg=10 as above, marking=true to draw a longitudinal marking and solid_color=:yellow.

Other optional arguments (and their default values) are

• Udof (true iff element has Udofs) wether to draw U-forces.
• draw_frame = false wether to draw the local reference frame of each element
• draw_marking = true wether to draw "longitudinal marking" along the element. Will only draw if style=:solid.
• nseg = 1 number of segments to display the shape of a deformed element
• solid_color = :yellow color of the surface if style=:solid
• line_color = :black color of the line if style=:sshape
• Uscale = 1. How many meter is a Newton per meter?