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

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)
stateXUA        = solve(DirectXUA{OX,OU,IA};initialstate,time=0:1.:5)

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 a State.
• time an 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 is a vector (over the iterations) of vectors (over the steps) of States of the model (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

A vector of length equal to that of time containing the state of the optimized model at each of these steps.

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 function 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 function 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 function 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::Function The gap function.
• gargs::NTuple Additional inputs to the gap function.
• mode::Function 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]) 
e1              = addelement!(model,DofConstraint,[n1],xinod=(1,),xfield=(:t1,),
                              λinod=1, λclass=:X, λfield=:λ1,gap=(x,t)->x[1]+.1,
                              mode=positive)
e2              = addelement!(model,QuickFix  ,[n1],inod=(1,),field=(:t1,),
                              res=(x,u,a,t)->0.4x.+.08+.5x.^2)
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.

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.
• class:Symbol :A for cost on A-dofs only, :I ("instant") otherwise.
• cost::Function if class==:I, cost(X,U,A,t,costargs...)→ℝ if class==:A, cost(A,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 tocost``

Requestable internal variables

• cost, the value of the cost.

Example

ele1 = addelement!(model,DofCost,[nod1],xinod=(1,),field=(:tx1,),
       class=:I,cost=(X,U,A,t;X0)->(X[1]-X0)^2,costargs=(;X0=0.27)

See also: 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::Function, where value(t::ℝ) → ℝ.

Requestable internal variables

• F, the value of the load.

Examples

using Muscade
model = Model(:TestModel)
node  = addnode!(model,𝕣[0,0])
e     = addelement!(model,DofLoad,[node];field=:tx,value=t->3t-1)

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!, studysingular, 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,X,U,A,t,gargs...)→ℝ. eleres is the output of the above-mentionned request to the target element. 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). X, U and A are the degrees of freedom of the element ElementType.

• gargs::NTuple Additional inputs to the gap function.

• mode::Function 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 function
• 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

@once gap(eleres,X,U,A,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, @once

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 function cost(eleres,X,U,A,t,costargs...)→ℝ. eleres is the output of the above-mentionned request to the target element. 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). X, U and A are the degrees of freedom of the element ElementType.
• [costargs::NTuple=() or NamedTuple] 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.

Requestable internal variables

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

• cost The value of cost(eleres,X,U,A,t,costargs...).
• 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 cost would extract the value of the ElementConstraint's function cost, while @request eleres(cost) refers to the value of a variable called cost in the target element.

Example

@once 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, @once

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!, studysingular, 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

Monitor <: AbstractElement

Debbuging tool: 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

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::Function, 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])
e = addelement!(model,QuickFix,[node1,node2];inod=(1,2),field=(:tx1,:tx1),
                res=(X,X′,X″,t)->Svector{2}(2*(X[1]-X[2]),2*(X[2]-X[1])) )

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, :U or :A.
• field::Symbol.
• cost::Function, 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])
e     = addelement!(model,SingleDofCost,[node];class=:X,field=:tx,
                    costargs=(3.,),cost=(x,t,three)->(x/three)^2)

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::Function, 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])
e     = addelement!(model,SingleUdof,[node];Xfield=:tx,Ufield=:utx,
                    costargs=(3.,),cost=(x,t,three)->(x/three)^2)

See also: DofCost, ElementCost

axis = Muscade.SpyAxis()

Spoof a GLMakie Axis 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 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{ORDER}

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, with feasibility line-search, to handle inequality constraints.
• SweepX{1} is implicit Euler, with feasibility line-search.
• SweepX{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 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(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. Dummy if ORDER<2
• maxiter=50 maximum number of equilibrium iterations at each step.
• maxΔx=1e-5 convergence criteria: norm of X. D
• 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.
• 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, studysingular, DirectXUA, FreqXU

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

default{:fieldname}(namedtuple,defval)

attempt to get a field fieldname from a NamedTuple. If namedtuple does not have such a field - or is not a NamedTuple, return defval.

λ,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{O}(X)

The vector X of Reals (possibly: ∂ℝs) is stripped of its partials, an revariated to the order precedence(X)+O.

TX = revariate(X)

revariates to the order precedence(X).

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

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

See also: compose

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

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 compose and Taylor

See also: compose, Taylor, revariate, fast

f  = QuadraticFunction(μ,σ)

μ and σ are 𝕣 (Float64)

y  = f(x) # == 1/2*((x-μ)/σ)^2
f  = QuadraticFunction(μ,σ)

Alternatively, μ can be a Function of time, in which case

y  = f(x,t) # == 1/2*((x-μ(t))/σ)^2

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: compose, 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

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

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!

colnormalize(a)

Euclidian-normalize the columns of an SMatrix

columnmatrix(v)

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

compose(Ty,x)

Compose automatic differentiation. For example Tx = revariate(x) Ty = f(Tx) y = compose(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

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.

See also: addelement!, addnode!

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: solve

diffed_lagrangian(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.

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

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

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

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.

nosecondorder(::Type{E<:AbstractElement})

Elements that define residual which would give excessive compilation and/or execution time if differentiated to the second order can implement a method after the below pattern to limit differentiation to first order:

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

off(t) → :off

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

See also: DofConstraint, ElementConstraint, equal, positive

positive(t) → :positive

A function which for any value t returns the symbol positive. Usefull 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!

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

@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.nosecondorder, 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, studyscale

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. 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!(S,i->abs(S.nzval[i])>tol)
sparser!(S,keep)
sparser!([S1,S2],1e-9)
sparser!(S1,S2],1e-9)

scale = studyscale(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 = studysingular(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.

To do so, 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!

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

@once tag f(x)= x^2

do not parse the definition of function f again if not modified. Using in a script, this prevents recompilations in Muscade or applications based on it when they receive such functions as argument.

tag must be a legal variable name, and unique to this invocation of @once

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 = @typeof(foo(args...[;kwargs...]))

Determine the inferred type and the returned type of the output[s] returned by the relevant method-instance of foo.