# RFEM / set' 12th / Las Flores

- Dejar armado para spring-mass:
- [x] modificar ejemplo springmass de ONSAS para sacar las matrices necesarias
- [x] "Formato" para intercambiar modelos entre ONSAS y JuliaReach
- [ ] Dejar integrado el ejemplo: matrices ensambladas se leen desde ReachabilityAnalysis.jl
- [ ] Agregar F sinusoidal al spring-mass.

- Dejar armado para heat:
- [ ] implementar heat element en ONSAS

## Formato

$$
Mu''(t) + Cu'(t) + f_{int}(u) = f_{ext}(t).
$$

en el caso lineal:
$$
f_{int} (u) = Ku
$$
entonces queda:
$$
Mu''(t) + Cu'(t) + Ku(t) = f_{ext}(t).
$$

### Del lado de RA (ReachabilityAnalysis.jl)

- Se llama la [función/ejemplo](https://github.com/ONSAS/ONSAS/blob/master/examples/onsasExample_springMass.m) de ONSAS `onsasExample_springMass( dirOnsas, scalarParams )` con dos argumentos:
- `dirOnsas`: string con la ruta relativa del ONSAS
- `scalarParams`: vector con los siguientes parámetros reales:
1. `k`: constante elástica de resorte
1. `c`: amortiguamiento de resorte
1. `m`: masa
1. `omegaBar`: frecuencia angular del forzamiento
1. `u0`: estiramiento inicial del resorte
1. `p0`: amplitud del forzamiento (ver función en línea 62)
1. `dt`: delta de tiempo

- Lee un `.MAT` y carga los datos:
- `M`: (nN x nN) matriz de masa (mass matrix) donde nN es el número de grados de libertad de Neumann
- `C`: (nN x nN) matriz de viscosidad (viscosity matrix)
- `K`: (nN x nN) matriz de rigidez (stiffness)
- `fext`: (1 x nT, opcional) vector de valores de factores de carga, donde nT es el número de pasos de tiempo
- `us`: (nN x nT, opcional) historia de desplazamientos
- `udots`: (nN x nT, opcional) historia de velocidades
- `neumdofs`: (nN x 1, opcional) vector de enteros con grados de libertad de las incógnitas del problema

octave-cli --eval "addpath(genpath('$ONSAS_PATH')), onsasExample_springMass('$ONSAS_PATH')"

octave-cli --eval "prueba('$ONSAS_PATH')"

using MAT

function read_fem(filename)
    M = matopen(filename, "M")
    C = matopen(filename, "C")
    K = matopen(filename, "K")
    return SecondOrderLinearContinuousSystem(M, C, K)
end

---

function input(t, times, evaluations)
    m = length(times)

    i = 1
    while i <= m && t > times[i]
        i += 1
    end
     
    if i <= m
        return Singleton(evaluations[i])
    else
        throw(ArgumentError("time t = $t is not contained in the " *
                "given time span, [$(times[1]), $(times[end])]"))
    end
end

---

por ahora lo hacemos asi, calculando la matriz inversa directamente:

$$
x'(t) = M^{-1} Ax(t) + M^{-1} Bu(t), x \in X, u \in U
$$

----

- Agregar a `MathematicalSystems.jl` un nuevo tipo de sistema que reprsenta

$$
Mx'(t) = Ax(t) + Bu(t), x \in X, u \in U
$$

$$
x'(t) = M^{-1}Ax(t) + M^{-1}Bu(t), x \in X, u \in U
$$

$$
\Phi = exp(M^{-1}A \delta)
$$

$$
\rho(d, \Phi X) = \rho(d, exp(M^{-1}A \delta) X) = \rho([exp(M^{-1}A \delta)]^T d, X) = \rho(exp(A^T (M^{-1})^T \delta) d, X)
$$

----

$$
\begin{bmatrix} \mathbf{A} & \mathbf{B} \\ \mathbf{C} & \mathbf{D} \end{bmatrix}^{-1} = \begin{bmatrix} \mathbf{A}^{-1}+\mathbf{A}^{-1}\mathbf{B}(\mathbf{D}-\mathbf{CA}^{-1}\mathbf{B})^{-1}\mathbf{CA}^{-1} & -\mathbf{A}^{-1}\mathbf{B}(\mathbf{D}-\mathbf{CA}^{-1}\mathbf{B})^{-1} \\ -(\mathbf{D}-\mathbf{CA}^{-1}\mathbf{B})^{-1}\mathbf{CA}^{-1} & (\mathbf{D}-\mathbf{CA}^{-1}\mathbf{B})^{-1} \end{bmatrix},
$$

$$
\begin{bmatrix} \mathbf{M}^{-1} & -\mathbf{M}^{-1} \mathbf{C} \\ 0 & \mathbf{I} \end{bmatrix}
$$

---------

# RFEM / set' 14th / Virtual

## Spring mass example

### Stacktrace (Marcelo)

$ pwd
/home/mforets/Projects/escritoire/2020/Week38

$ ./spring_mass.sh 
hofsfssssss =  1
dtCrit =  0.31834
error: run: file SCRIPT must exist and be a valid Octave scriptfile

$ cat spring_mass.sh 
octave-cli --eval "addpath(genpath('$ONSAS_PATH')), onsasExample_springMass('$ONSAS_PATH')"

$ echo $ONSAS_PATH
/home/mforets/Projects/ONSAS

Corriendolo desde el directorio de ejemplos:

$ pwd
/home/mforets/Projects/ONSAS/examples

$ octave-cli
GNU Octave, version 5.1.0
...

octave:1> onsasExample_springMass
hofsfssssss =  1
dtCrit =  0.31834
...
| Solving problem:  springMass
|  - input variables verification ...  done.
|  - Cleaning output directory ... done. 
ans = /home/mforets/Projects/ONSAS
|-------------------------------------------------|
| TimeSteps progress: 1|                   |  |
                        =====================
| end time-step 202 - (max,avg) iters: (  2, 1.02) | 
 ----------------------------------------------- 
  Analytical solution verification ... PASSED.
  Numerical solution error:   6.7654e-03 <   5.00e-02 
----------------------------------------------- 
  - writing report file ...  done. 
|-------------------------------------------------|
|  ONSAS finished in: 2.8e+00 seconds / 0.05 mins |
|=================================================|


octave:2> exit()

---

$ ./spring_mass.sh

...
| Solving problem:  springMass
|  - input variables verification ...  done.
  - Creating directory ./output/ ... done. 
|  - Creating output directory ...outputDir = ./output/springMass/
 done. 
ans = /home/mforets/Projects/escritoire/2020/Week38
|-------------------------------------------------|
| TimeSteps progress: 1|                   |  |
                        =====================
| end time-step 202 - (max,avg) iters: (  2, 1.02) | 
 ----------------------------------------------- 
  Analytical solution verification ... PASSED.
  Numerical solution error:   6.7654e-03 <   5.00e-02 
----------------------------------------------- 
  - writing report file ...  done. 
|-------------------------------------------------|
|  ONSAS finished in: 3.3e+00 seconds / 0.05 mins |
|=================================================|

---

$ ls
Manifest.toml  output  outputMatrices.mat  Project.toml  spring_mass2.sh  spring_mass.sh  Untitled.ipynb

[mforets@localhost Week38]$ julia
               _
   _       _ _(_)_     |  Documentation: https://docs.julialang.org
  (_)     | (_) (_)    |
   _ _   _| |_  __ _   |  Type "?" for help, "]?" for Pkg help.
  | | | | | | |/ _` |  |
  | | |_| | | | (_| |  |  Version 1.5.0-beta1.0 (2020-05-28)
 _/ |\__'_|_|_|\__'_|  |  Official https://julialang.org/ release
|__/                   |

julia> using MAT

julia> vars = matread("outputMatrices.mat")
Dict{String,Any} with 6 entries:
  "M"     => …
  "C"     => …
  "us"    => [0.0 0.000248295 … -0.00109504 -0.00245687]
  "K"     => …
  "fext"  => [0.0 9.94656 … -0.214363 9.73878]
  "udots" => [0.0 0.000248295 … -0.00109504 -0.00245687]

julia> size(vars["M"])
(1, 1)

julia> size(vars["M"]) == size(vars["C"]) == size(vars["K"])
true