Reference manual

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,noχ,noFB

See also: noFB

noχ

A constant, used by elements' residual or lagrangian as their 2nd output if they do not export any material memory (for example, plastic strain).

Example: return L,noχ,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

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.
• gₛ::𝕣=1. A scale for the gap.
• λₛ::𝕣=1. A scale for 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)
state           = solve(staticX;model,time=[0.],verbose=false) 
X               = state[1].X[1]

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

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)->X[1]^2

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

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

An element to apply costs on another element's element-results.

Named arguments to the constructor

• req a request for element-results for ElementType, resulting in the output eleres
• cost a cost function cost(eleres,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=(;) A named tuple of additional arguments to the cost function
• ElementType The named of the constructor for the relevant element
• elementkwargs... Additional named arguments to the ElementCost constructor are passed on to the ElementType constructor.

Requestable internal variables

• cost, the value of the cost.

Example

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

See also: SingleDofCost, DofCost, @request

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

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

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:

• no internal state.
• no initialisation.
• 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...]) → ℝ.
• costargs::NTuple
• derivative::Int 0, 1 or 2 - which 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

StaticX

A non-linear static solver for forward (not inverse, optimisation) FEM. The current implementation does not handle element memory.

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

initialstate    = initialize!(model)
state           = solve(StaticX;initialstate=initialstate,time=[0.,1.])

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 single state - obtained from a call to initialize!, or from a previous analysis
• time an AbstractVector vector of the times at which to compute the
• maxiter=50 maximum number of Newton-Raphson iteration at any given step
• maxΔx=1e-5 convergence criteria: a norm on the scaled X increment
• maxincrement=∞ convergence criteria: a norm on the scaled residual
• saveiter=false set to true so that the output state is a vector describing the states of the model at the last iteration (for debugging non-convergence)
• γ0=1. an initial value of the barrier coefficient for the handling of contact using an interior point method
• γfac1=0.5 at each iteration, the barrier parameter γ is multiplied
• γfac2=100. by γfac1*exp(-min(αᵢ)/γfac2)^2), where αᵢ is computed by the i-th interior point savvy element as αᵢ=abs(λ-g)/γ

Output

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

See also: solve, StaticXUA, initialize!

StaticXUA

A non-linear static solver for optimisation FEM. The current algorithm does not handle element memory.

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

initialstate    = initialize!(model)
stateX          = solve(StaticX  ;initialstate=initialstate,time=[0.,1.])
stateXUA        = solve(StaticXUA;initialstate=stateX)

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 vector of states, one for each load case in the optimization problem, obtained from one or several previous StaticX analyses
• maxAiter=50 maximum number of "outer" Newton-Raphson iterations over A
• maxΔa=1e-5 "outer" convergence criteria: a norm on the scaled A increment
• maxLa=∞ "outer" convergence criteria: a norm on the scaled La residual
• maxYiter=0 maximum number of "inner" Newton-Raphson iterations over X and U for every value of A. Experience so far is that these inner iterations do not increase performance, so the default is "no inner iterations".
• maxΔy=1e-5 "inner" convergence criteria: a norm on the scaled Y=[XU] increment
• maxLy=∞ "inner" convergence criteria: a norm on the scaled Ly=[Lx,Lu] residual
• saveiter=false set to true so that the output state is a vector describing the states of the model at the last iteration (for debugging non-convergence)
• γ0=1. an initial value of the barrier coefficient for the handling of contact using an interior point method
• γfac1=0.5 at each iteration, the barrier parameter γ is multiplied
• γfac2=100. by γfac1*exp(-min(αᵢ)/γfac2)^2), where αᵢ is computed by the i-th interior point savvy element as αᵢ=abs(λ-g)/γ

Output

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

See also: solve, StaticX

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.

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

ℤ (\bbℤ)

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.

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

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

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

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.

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

describe(model,spec)

Print out information about model. spec can be an EleID, a DofID, a NodID to describe an element, a dof or a node. spec can be :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 or :A

See also: solve

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

Elements must overload Muscade's doflist function. 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: lagrangian, residual

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

ilast = findlastassigned(state)

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

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

See also: solve

dofres = getdof(state;[class=:X],field=:somefield,nodID=[nodids...],[iders=0|1|2])

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,nder+1,nstate). If state is a scalar, the output dofres has size (ndof,nder+1).

See also: getresult, addnode!, solve

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

where class can be any of :X, :U, :A: get the number of dofs in the specified classes (or in the whole model).

model can be replaced with a concrete element type. The number of dofs is then per element.

See also: describe

eleres = getresult(state,req,eleids)

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. eleids can be either

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

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

initialstate = initialize!(model)

Finalize a model (invoquing addnode! and addelement! after initialize! will result in an error) and return an initial State (with all dofs set to zero, as starting point for an analysis.)

See also: addnode!, addelement!, solve

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

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.
• χ the element memory
• χcv a function used to built the updated element memory.
• 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
• χ an updated element memory. Return noχ if the element has no memory.
• 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: residual, doflist, @espy, ∂0, ∂1, ∂2, noχ, noFB

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.

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`](@ref), [`equal`](@ref), [`positive`](@ref)

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

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

The inputs and outputs to residual are the same as for lagrangian with two exceptions:

• residual does no have input Λ.
• the first output argument of residual is not a Lagrangian L but R,

a SVector{nXdof,R} where{R<:Real}, containing the element's contribution to the residual of the equations to be solved.

See lagrangian for the rest of the inputs and outputs.

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

Provides an order of magnitude for each type 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 in tens of meters, :rx are in radians, and A-dofs of field drag are of the order of 3. All other degrees of freedom are expected to be roughly of the order of 1.

Λ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

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

Execute an analysis using Solver, 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

This will call the method solve provided by the solver with solve(Solver,pstate,verbose,(dbg...,solver=Symbol(Solver));kwargs...)

See also: StaticX, StaticXUA, initialize!

toggle(condition,a,b)

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

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

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

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

@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 of interest 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.

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

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