doc of sign convention in elements. draw is graphic-engine agnostic

docs/src/Elements.md

@@ -19,11 +19,49 @@ In `Muscade`, the broad view is taken that anything that contributes to the Lagr
 
 Because "everything" is an element in `Muscade`, app developers can express a wide range of ideas through `Muscade`'s element API.
 
+## No internal variables
+
+In classical finite element formulations, plastic strain is implemented by letting an element have plastic strain and a hardening parameter as an *internal variable* at each quadrature point of the element.  Internal variables are not degrees of freedom.  Instead, they are a memory of the converged state of the element at the previous load or time step, used to affect the residual computed at the present time step.  Internal variables are also used when modeling friction and damage processes.
+
+`Muscade` does not allow elements to have internal variables. The reason is that problem involing dofs of class `U` are not causal: our estimation of the state of a system a step `i` also depends on on measurements taken at steps `j` with `j>i`.  Hence "sweep" procedures, that is, procedures that solve a problem one load or time step at a time, are not applicable to such problems.  Solvers must hence solve for all dofs and steps at once.  When using Newton-Raphson iterations to solve problems of this class, internal variables make the Hessian matrix full: the value of a stress at a given stress depends (through the internal variable) on the strain at all preceeding steps. Or more formaly: Internal variables transforms the problem from *differential* to *integral*. A full matrix quickly leads to impossibly heavy computations.
+
+To model phenomena usualy treated using internal variable, it is necessary in `Muscade` to make the "internal" variable into a degree of freedom, and describe the equation of evolution of this degree of freedom. See [`examples/DecayAnalysis.jl`](DecayAnalysis.md) for an example of implementation.  
+!!! warning
+    Because the equations of evolution involves first order time derivative, one can not use a static solver in combination with such elements.
+
+## Sign convention in elements
+
+Starting with matrix methods in structural analysis, the traditional convention is that in an equation of the form
+
+```math
+K \cdot ΔX = R\\
+X \leftarrow X + ΔX
+```
+
+``K`` is (typicaly) symmetric positive definite, ``ΔX`` are incremental nodal displacements, and ``R`` are *external* loads applied to the structure (a positive load tends to induce a positive displacement).  As a consequence, when an element is implemented within this convention, the element must return its stiffness ``K`` and its "internal reaction forces" ``R_i``: a bar that is elongated reacts by pulling its ends inwards. The forces are "el-on-nod" (element on node).  `Muscade` uses the same convention for the *description of models*.  
+
+However, the *implementation of elements* in `Muscade` uses another convention.  This is because `Muscade` optimizes a Lagrangian, relative to a set of variables here collectively denoted as ``Z``.  A Newton step for seeking to make ``L(Z)`` stationary is naturaly written as
+
+```math
+G = \frac{\partial L}{\partial Z}\\
+H = \frac{\partial G}{\partial Z}\\
+H \cdot ΔZ = G\\
+Z \leftarrow Z - ΔX
+```
+
+Note the minus sign on the last lign.  As a consequence of this minus sign, in `Muscade`, an element returns ``R_e=-R_i``, and ``K`` is computed (by automatic differentiation, invisible to the element developer) as
+
+```math
+K = \frac{\partial R_e}{\partial X}
+```
+
+This implies that ``R_e`` are the "external forces": to elongate a bar one must pull its ends outwards. The forces are "nod-on-el" (node on element).  This has one implication that may be surprising: an element that for example implements an (external) point load ``F`` must return ``R_e = -F``, note the minus sign, so that the user of the element will interpret ``F`` as a classic external load.  The same applies to elements that connect unknown external loads ``U`` to the equilibrium equations.  Such an element must return ``R_e = -U``.  Further, elements that return a Lagrangian (see below) must return ``L = Q + \Lambda \cdot R_e``, note the plus sign.
+
 ## API
 
 The implementation of a element requires 
 
-- A **[`struct`](@ref struct)** defining the element. 
+- A **[`DataType`](@ref struct)** defining the element. 
 - A **[constructor](@ref constructor)** which is called when the user adds an element to the model, and constructs the above `struct`.
 - **[`Muscade.doflist`](@ref)** specifies the degrees of freedom (dofs) of the element.
 - **[`Muscade.residual`](@ref)** (either this of [`Muscade.lagrangian`](@ref)) takes element dofs as input and returns the element's additive contribution to the residual of a non-linear system of equations,
@@ -32,16 +70,6 @@ The implementation of a element requires
 
 Each element must implement *either* [`Muscade.lagrangian`](@ref) *or* [`Muscade.residual`](@ref), depending on what is more natural: a beam element will implement [`Muscade.residual`](@ref) (element reaction forces as a function of nodal displacements), while an element representing a strain sensor will implement [`Muscade.lagrangian`](@ref) (log-of the probability density of the strain, given an uncertain measurement).
 
-## No internal variables
-
-In classical finite element formulations, plastic strain is implemented by letting an element have plastic strain and a hardening parameter as an *internal variable* at each quadrature point of the element.  Internal variables are not degrees of freedom.  Instead, they are a memory of the converged state of the element at the previous load or time step, used to affect the residual computed at the present time step.  Internal variables are also used when modeling friction and damage processes.
-
-`Muscade` does not allow elements to have internal variables. The reason is that problem involing dofs of class `U` are not causal: our estimation of the state of a system a step `i` also depends on on measurements taken at steps `j` with `j>i`.  Hence "sweep" procedures, that is, procedures that solve a problem one load or time step at a time, are not applicable to such problems.  Solvers must hence solve for all dofs and steps at once.  When using Newton-Raphson iterations to solve problems of this class, internal variables make the Hessian matrix full: the value of a stress at a given stress depends (through the internal variable) on the strain at all preceeding steps. Or more formaly: Internal variables transforms the problem from *differential* to *integral*. A full matrix quickly leads to impossibly heavy computations.
-
-To model phenomena usualy treated using internal variable, it is necessary in `Muscade` to make the "internal" variable into a degree of freedom, and describe the equation of evolution of this degree of freedom. See [`examples/DecayAnalysis.jl`](DecayAnalysis.md) for an example of implementation.  
-!!! warning
-    Because the equations of evolution involves first order time derivative, one can not use a static solver in combination with such elements.
-
 ## [DataType](@id struct)
 
 For a new element type `MyELement`, the datatype is defined as
@@ -300,9 +328,12 @@ See the page on [automatic differentiation](Adiff.md).
 
 Elements *can* implement a [`Muscade.draw`](@ref) method. If no method is implemented, the element will be invisible if the user requests a drawing of the element.
 
+None of `Muscade` built-in elements implement methods for `draw`: because `Muscade` has no inherent interpretation of the various `X` dofs, there is no graphical representation associated to them.  On the other hand, it might make sense for an app developer (giving an interpretation to various dofs) to create such methods.
+
+Because `Muscade` provides no implementation of `draw` (with the exception of some demo elements), `Muscade` does not prescribe the use of any specific graphic package.  See [`Makie.jl`](https://docs.makie.org/) and [`WriteVTK.jl`](https://juliavtk.github.io/WriteVTK.jl/stable/) for candidates.
+
 While the API may remind that of [`Muscade.lagrangian`](@ref), there is one significant difference: 
-because it is more efficient to create few graphical object (few calls to `lines!`, `scatter!`) etc., the element's method for `draw` 
-will be called once to draw several elements of the same type. Multiple lines can be drawn in one call to `lines!` by using `NaN`s to "lift the pen".
+because it is more efficient to create few graphical object (in `Makie`: few calls to `lines!`, `scatter!`) etc., the element's method for `draw` will be called once to draw several elements of the same type. In `Makie` multiple lines can be drawn in one call to `lines!` by using `NaN`s to "lift the pen".
 
 ### Keyword arguments
 

docs/src/Solvers.md

@@ -1,12 +1,12 @@
-# Solvers
+# [Solvers](@id solvers)
 
 ## `SweepX`
 
 `SweepX{O}` is a non-linear solver for differential equations of order `O` in time. This can be used for static and quasi static problems without hysterertic behaviour (plasticity, friction).
 
 `SweepX{1}` is an implicit-Euler solver for differential equations of order `1` in time. This *must* be used for viscous problems, and for static and quasi static problem with hysteretic behaviour. The reason for this is that Muscade does not allow elements to have element-internal "state" variables (plastic strain, shear-free position for dry friction). Hence, where elements implement such physics, this is done by introducing the "state" as a degree of freedom of the element, and a corresponding equation.  This equation is the equation of evolution of the "state" variable, which involves the first order derivative of the variable in question *even in a static problem*.
 
-`SweepX{2}` is  Newmark-β solver for differential equations of order `2` in time. A typical application is in structural dynamics. 
+`SweepX{2}` is a Newmark-β solver for differential equations of order `2` in time. A typical application is in structural dynamics. 
 
 `SweepX` solves forward FEM problems (not optimisation-FEM) (see [Theory](@ref)).  However, `SweepX` can be applied to models that have ``U``- and ``A``-dofs. This is handled as follows: One input to `SweepX` is a `State`, which can come from [`initialize!`](@ref) or from the output of another solver. `SweepX` will keep the ``U``- and ``A``-dofs to the value in the input `State`. `initialize!` sets all dofs to zero, so when `SweepX` is given a `State` produced by `initialize!` the analysis starts with ``X``-dofs equal to zero, and ``U``- and ``A``-dofs are kept zero throughout the analysis. 
 
@@ -15,6 +15,7 @@
 See the reference manual [`SweepX`](@ref).
 
 ## `DirectXUA`
+
 `DirectXUA` is a solver for non-linear, static (`OX=0`), first order (`OX=1`) or dynamic (`OX=2`), optimisation-FEM problems. The same remarks on "state" variables and the choice of `OX` as for `SweepX` apply here. 
 
 `DirectXUA` is designed for load and model parameter identification. Given a model with costs (and possibly constraints) on U- and ``A``-dofs, the solver will determine response (``X``-dofs) and unknown loads for each step (``U``-dofs). If (`IA=1`), the algorithm will also estimate model parameters for the whole history (``A``-dofs).
@@ -27,4 +28,6 @@ See the reference manual [`SweepX`](@ref).
 !!! info
     Currently, the solver does not handle inequality constraints.
 
-See the reference manual [`DirectXUA`](@ref).
+See the reference manual [`DirectXUA`](@ref).
+
+

docs/src/Theory.md

@@ -1,5 +1,9 @@
 # Theory
 
+## Introduction
+
+This section presents the theory of how FEM-optimization problems are *defined* or formulated.  This is of interest both for the user of a `Muscade`-based application, and for someone implementing new elements to create such an application.
+
 ## Classes of degrees of freedom
 
 `Muscade` introduces 3 *classes* of degrees of freedom (dofs). 
@@ -94,7 +98,7 @@ When doing optimization-FEM in `Muscade`, the elements provide additive contribu
 
 ## Physical and optimisation constraints
 
-*Physical constraints* (including contact or Dirichlet/essential boundary conditions) are added to the model by using an element that adds adds a dof of class ``X`` to the model.  This new dof ``X_λ`` is a Lagrange multiplier for the constraint. If we note ``R``, ``X`` and ``Λ`` the list of dofs before adding the constraint element and ``R^*``, ``X^*`` and ``Λ^*`` after, then  
+**Physical constraints** (including contact or Dirichlet/essential boundary conditions) are added to the model by using an element that adds adds a dof of class ``X`` to the model.  This new dof ``X_λ`` is a Lagrange multiplier for the constraint. If we note ``R``, ``X`` and ``Λ`` the list of dofs before adding the constraint element and ``R^*``, ``X^*`` and ``Λ^*`` after, then  
 
 ```math
 L^*(\Lambda^*,X^*,U,A) = Q(X,U,A)  + R^*(X^*,U,A) \cdot Λ^* 
@@ -110,7 +114,7 @@ R^*(X^*,U,A) &= \left[R(X,U,A) - ∇g_x(X,U,A) \cdot X_λ \; , \; g_x(X,U,A)\rig
 \end{aligned}
 ```
 
-*Optimisation constraints* allow to define that some situations are impossible, or inacceptable. For example, in a design optimisation analysis, excessive stresses would lead to failure, so would be constrained to remain under a given threshold. In a tracking or optimal control problem, an actuator force may not exceed some limit for the actuator's capacity.
+**Optimisation constraints** allow to define that some situations are impossible, or inacceptable. For example, in a design optimisation analysis, excessive stresses would lead to failure, so would be constrained to remain under a given threshold. In a tracking or optimal control problem, an actuator force may not exceed some limit for the actuator's capacity.
 
 Optimisation constraints that have to be verified at every step (for example stresses that must remain below a critical level, at any time) require a Lagrange multiplier that changes over time, and that is thus of class ``U``. Optimisation constraints that act only on ``A``-dofs (for example, there is a limit to the strength of steel we can order)  require a Lagrange multiplier of class ``A``.  With optimisation constraints, the Lagrangian is of the form
 
@@ -130,7 +134,7 @@ Q^*(X,U^*,A^*) &= Q(X,U,A) +  g_u(X,U,A) \cdot U_\lambda + g_a(A) \cdot A_λ \\
 
 ## Literature
 
-For more details on the theory and example of applications, see [maincon04a](@cite), [maincon04b](@cite), [maree04](@cite), [barnardo04](@cite), [hauser06](@cite), [wu07](@cite), [maincon08](@cite), [hauser08](@cite), [wu08](@cite), [wu09](@cite), [maincon13](@cite). 
+For more details on the theory and examples of applications, see [maincon04a](@cite), [maincon04b](@cite), [maree04](@cite), [barnardo04](@cite), [hauser06](@cite), [wu07](@cite), [maincon08](@cite), [hauser08](@cite), [wu08](@cite), [wu09](@cite), [maincon13](@cite).