fixes

nafems4.md

-$\renewcommand\vec[1]{\mathbf{#1}}$
+dnl $$\renewcommand\mathbf[1]{\mathbf{#1}}$$
 
 # Background and introduction
 
 Finite elements are like magic to me. I mean, I can follow the whole derivation of the equations, from the strong, weak and variational formulations of the equilibrium equations for the mechanical problem (or the energy conservation for heat transfer) down to the algebraic multigrid preconditioner for the inversion of the stiffness matrix passing through Sobolev spaces and the grid generation. Then I can sit down and program all these steps into a computer, including the shape functions and its derivatives, the assembly of the discretised stiffness matrix assembly, the numerical solution of the system of equations and the computation of the gradient of the solution. Yet, the fact that all these a-priori unconnected steps once gets a pretty picture that resembles reality is still astonishing to me.
 
 
-Again, take all this information as coming from a fellow that has already taken such a journey from college’s pencil and paper to real engineering cases involving complex numerical calculations. And developing, in the meantime, both an actual working finite-element [back-end](https://www.seamplex.com/fino) and [front-end](https://www.caeplex.com) from scratch (the dean of the engineering school I attended used to say “It is not the same to read than to write manuals, andwe should aim at writing.”).
+Again, take all this information as coming from a fellow that has already taken such a journey from college’s pencil and paper to real engineering cases involving complex numerical calculations. And developing, in the meantime, both an actual working finite-element [back-end](https://www.seamplex.com/fino) and [front-end](https://www.caeplex.com) from scratch (the dean of the engineering school I attended used to say “It is not the same to read than to write manuals, and we should aim at writing.”).
 
 
 ## Tips and tricks
 
 In any case, what we should understand (and imagine) is that external forces lead to internal stresses. And in any three-dimensional body subject to such external loads, the best way to represent internal stresses is through a $3 \times 3$ _stress tensor_. This is the first point in which we should not fear mathematics. Trust me, it will pay back later on.
 
-Does the term _tensor_ scare you? It should not. A tensor is just a slightly more complex vector, and I assume you are not afraid of vectors. If you recall, a vector somehow generalises the idea of a scalar in the following sense: a given vector $\vec{v}$ can be projected into any direction $\vec{n}$ to obtain a scalar $p$. We call this scalar $p$ the “projection” of the vector $\vec{v}$ in the direction $\vec{n}$. Well, a tensor can be also projected into any direction $\vec{n}$. The difference is that instead of a scalar, a vector is now obtained.
+Does the term _tensor_ scare you? It should not. A tensor is just a slightly more complex vector, and I assume you are not afraid of vectors. If you recall, a vector somehow generalises the idea of a scalar in the following sense: a given vector $\mathbf{v}$ can be projected into any direction $\mathbf{n}$ to obtain a scalar $p$. We call this scalar $p$ the “projection” of the vector $\mathbf{v}$ in the direction $\mathbf{n}$. Well, a tensor can be also projected into any direction $\mathbf{n}$. The difference is that instead of a scalar, a vector is now obtained.
 
 Let me introduce then the three-dimensional stress tensor:
 
  * There exists a particular coordinate system in which the stress tensor is diagonal, i.e. all the shear stresses are zero. In this case, the three diagonal elements are called the principal stresses, which happen to be the three [eigenvalues](https://en.wikipedia.org/wiki/Eigenvalues_and_eigenvectors) of the stress tensor. The basis of the coordinate system that renders the tensor diagonal are the eigenvectors.
 
 
-What does this all have to do with mechanical engineering? Well, once we know what the stress tensor is for every point of a solid, in order to obtain the internal forces per unit area acting in a plane passing through that point and with a normal given by the direction $\vec{n}$, all we have to do is “project” the stress tensor through $\vec{n}$. In plain simple words:
+What does this all have to do with mechanical engineering? Well, once we know what the stress tensor is for every point of a solid, in order to obtain the internal forces per unit area acting in a plane passing through that point and with a normal given by the direction $\mathbf{n}$, all we have to do is “project” the stress tensor through $\mathbf{n}$. In plain simple words:
 
 > If you can compute the stress tensor at each point of our geometry, then... Congratulations! You have solved the solid mechanics problem.
 
 So for each time of the operational transient, the pipes are subject to
 
  a. an uniform internal pressure\ $p_i(t)$ that depends on time,
- b. a uniform internal temperature $T_i(t)$ that gives rise to a non-trivial time-dependent temperature distribution\ $T(\vec{x},t)$ in the bulk of the pipes, and
- c. internal distributed forces\ $\vec{f}=\rho \cdot \vec{a}$ at those times where the design earthquake is assumed to act.
+ b. a uniform internal temperature $T_i(t)$ that gives rise to a non-trivial time-dependent temperature distribution\ $T(\mathbf{x},t)$ in the bulk of the pipes, and
+ c. internal distributed forces\ $\mathbf{f}=\rho \cdot \mathbf{a}$ at those times where the design earthquake is assumed to act.
 
 ## Thermal transient {#sec:thermal}
 
 
 
 
-As an example, let us consider the system depicted in\ [@fig:valve-cad1;@fig:valve-scls1] where there is a stainless-carbon steel interface at the discharge of the valves. Instead of solving the transient heat-conduction problem with the internal temperature of the pipes equal to the temperature of the water in the reference transient condition of the power plant and an external condition of natural convection to the ambient temperature in the whole mesh of\ [@fig:valve-mesh1], a reduced model consisting of half of one of the two valves and a small length of the pipes at both the valve inlet and outlet is used. Once the temperature distribution\ $\hat{T}(\vec{x},t)$ for each time is obtained in the reduced mesh ([@fig:valve-temp], which has the origin at the centre of the valve), the actual temperature distribution\ $T(\vec{x},t)$ is computed by an algebraic generalisation of $\hat{T}(\vec{x},t)$ in the full coordinate system (where the origin is shown in\ [@fig:valve-cad1]). As stated above, those locations which are not covered by the reduced model are generalised with a time-dependent uniform temperature which is the average of the inner and outer temperatures at the inlet and outlet of the reduced mesh. The result is illustrated in figure\ [@fig:valve-gen].
+As an example, let us consider the system depicted in\ [@fig:valve-cad1;@fig:valve-scls1] where there is a stainless-carbon steel interface at the discharge of the valves. Instead of solving the transient heat-conduction problem with the internal temperature of the pipes equal to the temperature of the water in the reference transient condition of the power plant and an external condition of natural convection to the ambient temperature in the whole mesh of\ [@fig:valve-mesh1], a reduced model consisting of half of one of the two valves and a small length of the pipes at both the valve inlet and outlet is used. Once the temperature distribution\ $\hat{T}(\mathbf{x},t)$ for each time is obtained in the reduced mesh ([@fig:valve-temp], which has the origin at the centre of the valve), the actual temperature distribution\ $T(\mathbf{x},t)$ is computed by an algebraic generalisation of $\hat{T}(\mathbf{x},t)$ in the full coordinate system (where the origin is shown in\ [@fig:valve-cad1]). As stated above, those locations which are not covered by the reduced model are generalised with a time-dependent uniform temperature which is the average of the inner and outer temperatures at the inlet and outlet of the reduced mesh. The result is illustrated in figure\ [@fig:valve-gen].
 
 ::::: {#fig:valve-temp-gen}
 ![Reduced mesh around the valve refined around the interface where the transient heat conduction problem is solved.](valve-temp.png){#fig:valve-temp width=80%}
    5.767255
 ```
 
-Did I convince you? More or less, right? The third eigenvalue seems to fit. Let us not throw all of our beloved linearity away and dig in further into the subject. There are still two important issues to discuss which can be easily addressed using fresh-year linear algebra (remember, do not fear maths!). First of all, even though principal stresses are not linear with respect to the sum they are linear with respect to pure multiplication. Once more, think what happens to the the eigenvalues and eigenvectors of a single stress tensor as all its elements are scaled up or down by a real scalar. They are the same! So, for example, the [Von\ Mises stress](https://en.wikipedia.org/wiki/Von_Mises_yield_criterion) (which is a combination of the principal stresses) of a beam loaded with a force\ $\alpha \cdot \vec{F}$ is\ $\alpha$ times the stress of the beam loaded with a force\ $\vec{F}$. Please test this hypothesis by playing with your favorite FEM solver an play. Or even better, take a look at the stress invariants $I_1$, $I_2$ and $I_3$ (you can search online or peek into the source code of [Fino](https://www.seamplex.com/fino) and grep for the routine called `fino_compute_principal_stress()`) and see (using paper and pencil!) how they scale up if the individual elements of the stress tensor are scaled by a real factor\ $\alpha$.
+Did I convince you? More or less, right? The third eigenvalue seems to fit. Let us not throw all of our beloved linearity away and dig in further into the subject. There are still two important issues to discuss which can be easily addressed using fresh-year linear algebra (remember, do not fear maths!). First of all, even though principal stresses are not linear with respect to the sum they are linear with respect to pure multiplication. Once more, think what happens to the the eigenvalues and eigenvectors of a single stress tensor as all its elements are scaled up or down by a real scalar. They are the same! So, for example, the [Von\ Mises stress](https://en.wikipedia.org/wiki/Von_Mises_yield_criterion) (which is a combination of the principal stresses) of a beam loaded with a force\ $\alpha \cdot \mathbf{F}$ is\ $\alpha$ times the stress of the beam loaded with a force\ $\mathbf{F}$. Please test this hypothesis by playing with your favorite FEM solver an play. Or even better, take a look at the stress invariants $I_1$, $I_2$ and $I_3$ (you can search online or peek into the source code of [Fino](https://www.seamplex.com/fino) and grep for the routine called `fino_compute_principal_stress()`) and see (using paper and pencil!) how they scale up if the individual elements of the stress tensor are scaled by a real factor\ $\alpha$.
 
 The other issue is that even though in general the eigenvalues of the sum of two matrices are not the same as the eigenvalues of the matrix sum, there are some cases when they are. In effect, if two matrices\ $A$ and\ $B$ commute, i.e. their product is commutative
 
 
 ### Solving
 
-The linear FEM problem leads of course of a system of\ $NG$ linear equations, cast in matrix form by the stiffness matrix\ $K$ and a right-hand size vector\ $\vec{b}$ containing the loads (both volumetric and the ones at the surfaces from the boundary conditions):
+The linear FEM problem leads of course of a system of\ $NG$ linear equations, cast in matrix form by the stiffness matrix\ $K$ and a right-hand size vector\ $\mathbf{b}$ containing the loads (both volumetric and the ones at the surfaces from the boundary conditions):
 
-$$K \cdot \vec{u} = \vec{b}$${#eq:kub}
+$$K \cdot \mathbf{u} = \mathbf{b}$${#eq:kub}
 
-The objective of the solver is to find the vector\ $u$ of nodal displacements that satisfy the momentum equilibrium. Luckily (well, not purely by chance but by design) the stiffness matrix is almost empty. It is called a [sparse matrix](https://en.wikipedia.org/wiki/Sparse_matrix) because most of its elements are zero. If it was fully filled, then a problem with just 100k nodes would need more than 700\ Gb of RAM just to store the matrix elements, rendering FEM as practically impossible. And even though the stiffness matrix is sparse, its inverse is not so we cannot solve the elastic problem by “inverting” the matrix. Particular methods to represent and more importantly to solve linear systems involving these kind of matrices have been developed, which are the methods used by finite-element (and the other finite-cousins) programs. These methods are [iterative](https://en.wikipedia.org/wiki/Iterative_method) in nature, meaning that at each step of the iteration the numerical answer improves, i.e. $\left|K \cdot \vec{u} - \vec{b}\right| \rightarrow 0$. Even though in particular for [Krylov-subspace](https://en.wikipedia.org/wiki/Krylov_subspace) methods, $K \cdot \vec{u} - \vec{b} = \vec{0}$ after\ $NG$ steps, a few dozen of iterations are usually enough to assume that the problem is effectively solve (i.e. the residual is less than a certain threshold).
+The objective of the solver is to find the vector\ $u$ of nodal displacements that satisfy the momentum equilibrium. Luckily (well, not purely by chance but by design) the stiffness matrix is almost empty. It is called a [sparse matrix](https://en.wikipedia.org/wiki/Sparse_matrix) because most of its elements are zero. If it was fully filled, then a problem with just 100k nodes would need more than 700\ Gb of RAM just to store the matrix elements, rendering FEM as practically impossible. And even though the stiffness matrix is sparse, its inverse is not so we cannot solve the elastic problem by “inverting” the matrix. Particular methods to represent and more importantly to solve linear systems involving these kind of matrices have been developed, which are the methods used by finite-element (and the other finite-cousins) programs. These methods are [iterative](https://en.wikipedia.org/wiki/Iterative_method) in nature, meaning that at each step of the iteration the numerical answer improves, i.e. $\left|K \cdot \mathbf{u} - \mathbf{b}\right| \rightarrow 0$. Even though in particular for [Krylov-subspace](https://en.wikipedia.org/wiki/Krylov_subspace) methods, $K \cdot \mathbf{u} - \mathbf{b} = \mathbf{0}$ after\ $NG$ steps, a few dozen of iterations are usually enough to assume that the problem is effectively solve (i.e. the residual is less than a certain threshold).
 
 So the question is... how hard is to solve a sparse linear problem? Well, the number of iterations needed to attain convergence depends on the [spectral radius](https://en.wikipedia.org/wiki/Spectral_radius) of the stiffness matrix\ $K$, which in turns depend on the elements themselves, which depend mainly on the parameters of the elastic problem. But it also depends on how sparse\ $K$ is, which changes with the element topology and order. [@Fig:test] shows the structure of two stiffness matrices for the same linear elastic problem discretised using exactly the same number of nodes with first and second-order elements. The matrices have the same size (because the number of nodes is the same) but the former is more sparse than the latter. Hence, it would be harder (in computational terms of CPU and RAM) to solve the second-order case.
 
 
 ### Stress computation
 
-In the [displacement-based formulation](http://web.mit.edu/kjb/www/Books/FEP_2nd_Edition_4th_Printing.pdf), the solver just “solves” for the displacements\ $\vec{u}(\vec{x})$ which are the principal unknowns. But from\ [@sec:tensor] we know that we actually “solve” the problem when we have the stress tensors at every location\ $\vec{x}$, which are the secondary unknowns. So the FEM program has to compute the stresses from the displacements using the materials’ constitutive equations which involve the Young Modulus\ $E$, the Poisson ratio\ $\nu$ and the spatial derivatives of the displacements\ $\vec{u}$. This sounds easy, as we (well, the solver) knows what the shape functions are for each element and then it is a matter of computing nine derivatives (that of the three displacements in each of the three directions) and multiplying by something involving\ $E$ and\ $\nu$. Yes, but there is a catch. As the displacements\ $u$ are computed at the nodes, we would like to also have the stresses at the nodes. However,
+In the [displacement-based formulation](http://web.mit.edu/kjb/www/Books/FEP_2nd_Edition_4th_Printing.pdf), the solver just “solves” for the displacements\ $\mathbf{u}(\mathbf{x})$ which are the principal unknowns. But from\ [@sec:tensor] we know that we actually “solved” the problem when we have the stress tensors at every location\ $\mathbf{x}$, which are the secondary unknowns. So the FEM program has to compute the stresses from the displacements using the materials’ [strain-stress constitutive equations](https://en.wikipedia.org/wiki/Constitutive_equation#Stress_and_strain) which involve the Young Modulus\ $E$, the Poisson ratio\ $\nu$ and the spatial derivatives of the displacements\ $\mathbf{u}=[u,v,w]$. This sounds easy, as we (well, the solver) knows what the shape functions are for each element and then it is a matter of computing nine derivatives (that of the three displacements in each of the three directions) and multiplying by something involving\ $E$ and\ $\nu$. Yes, but there is a catch. As the displacements\ $u$, $v$ and\ $w$ are computed at the nodes, we would like to also have the stresses at the nodes. However,
 
- i. the displacements\ $\vec{u}(\vec{x})$ are not differentiable at the nodes, and
- ii. if the node belongs to a material interface, the\ $E$ and\ $\nu$ stuff is not defined.
+ i. the displacements\ $\mathbf{u}(\mathbf{x})$ are not differentiable at the nodes, and
+ ii. if the node belongs to a material interface, neither\ $E$ nor\ $\nu$ are defined.
 
-Let us focus on the first item and leave the second one for a separate discussion in\ [@sec:two-materials]. The finite-element method computes the principal unknowns at the nodes and then says “interpolate the nodal values inside each element using its shape functions.” It sounds (and it is) great, but a node belongs to more than one element (you can now imagine a 3D mesh composed of tetrahedra but you can also simplify your mind pictures by thinking in just one dimension: a node is shared by two segments). So the slope of the interpolation when we move from the node into one element might (and it never is) the same as when we move from the same node into another element. Neither in linear, nor quadratic nor higher-order elements. So what is the derivative of the displacement\ $v$ in the\ $y$ direction with respect to the $z$ coordinate? One answer would be the average of the derivatives computed in each of the elements that share a common node. But what if one of the elements is very small? Or it has a very bad quality (i.e. it is deformed in one direction) its derivatives cannot be trusted? Should the “average” be weighted? How?
+Let us focus on the first item and leave the second one for a separate discussion in\ [@sec:two-materials]. The finite-element method computes the principal unknowns at the nodes and then says “interpolate the nodal values inside each element using its shape functions.” It sounds (and it is) great, but a node belongs to more than one element (you can now imagine a 3D mesh composed of tetrahedra but you can also simplify your mind pictures by thinking in just one dimension: a node is shared by two segments). So the slope of the interpolation when we move from the node into one element might (and it never is) the same as when we move from the same node into another element. Stop and think. Now let us take a look at [@fig:derivatives]. It shows four finite-element solutions for a certain problem that has a cosine as the analytical solution. All the cases use eight elements: either uniformly distributed along the domain $x \in [-1:1]$ or only one element in the negative half and the remaining seven evenly distributed between zero and one, configuring a rather extreme non-uniform grid to better illustrate the point. Both first and second-order elements were used for each mesh, hence obtaining four combinations. We know the derivative of the analytical solution is zero at $x=0$. In the uniform cases of [@fig:slab-1-0;@fig:slab-2-0], if we took the average of the derivatives at each side of the origin we would obtain zero as expected. But in the non-uniform cases\ [@fig:slab-1-1;@fig:slab-2-1], plain averaging fails even in the quadratic case.
 
 ::::: {#fig:derivatives}
 ![](slab-1-0.svg){#fig:slab-1-0 width=48%}\  
 Solution of a problem using FEM using eight linear/quadratic uniform/non-uniform elements. The reference solution is a cosine. Plain averaging works for uniform grids but fails in the non-uniform cases.
 :::::
 
+Now proceed to picturing the general three-dimensional cases with unstructured tetrahedra. What is the derivative of the displacement\ $v$ in the\ $y$ direction with respect to the $z$ coordinate at a certain node shared by many tetrahedra? What if one of the elements is very small? Or it has a very bad quality (i.e. it is deformed in one direction) its derivatives cannot be trusted? Should we still average? Should this average be weighted? How?
+
 Detailed mathematics show that the location where the derivatives of the interpolated displacements are closer to the real (i.e. the analytical in problem that have it) solution are the elements’ [Gauss points](https://en.wikipedia.org/wiki/Gaussian_quadrature). Even better, the material properties at these points are continuous (they are usually uniform but they can depend on temperature for example) because, unless we are using weird elements, there are no material interfaces inside elements. But how to manage a set of stresses given at the Gauss points instead of at the nodes? Should we use one mesh for the input and another one for the output? What happens when we need to know the stresses on a surface and not just in the bulk of the solid? There are still no one-size-fits-all answers. There is a very interesting [blog post](http://tor-eng.com/2017/11/practical-tips-dealing-stress-singularities/) by Nick Stevens that addresses the issue of stresses computed at sharp corners. What does your favourite FEM program do with such a case?
 
 In any case, this step takes a non-negligible amount of time. The most-common approach, i.e. the node-averaging method is driven mainly by the number of nodes of course. So all-in-all, these are the reasons to use the number of nodes instead of the numbers of elements as a basic parameter to measure the complexity of a FEM problem.
 :::::
 
 
-To simplify the discussion that follows, let us replace for one moment all the full $3 \times 3$ tensor and the nine partial derivatives of the displacement by just one strain\ $\epsilon$ and let the strain-stress relationship to be the scalar expression
+To simplify the discussion that follows, let us replace for one moment all the full $3 \times 3$ tensor and the nine partial derivatives of the displacement by just one strain\ $\epsilon$ and let the [linear elastic strain-stress relationship](https://en.wikipedia.org/wiki/Elasticity_(physics)#Linear_elasticity) to be the scalar expression
 
 $$
 \sigma = E \cdot \epsilon
 
  This would be like forcing a non-differentiable function to behave smoothly, or
 
- 3. internally duplicate the nodes at the interface and compute one stress for each material. This would result in a stress distribution which is not a real function because the same location\ $\vec{x}$ will be associated to more than one stress value, or
+ 3. internally duplicate the nodes at the interface and compute one stress for each material. This would result in a stress distribution which is not a real function because the same location\ $\mathbf{x}$ will be associated to more than one stress value, or
 
  4. duplicate the surface elements at the interfaces and solve the problem using a contact formulation. This would render the problem non-linear add the complexity of having to find appropriate penalty coefficients.
 
 So here comes the five Feynmann-Ohno questions:
 
  1. Why do you want to perform a parametric computation?  
-     So we can study how the linearised stresses change with the inclusion of a branch of a certain diameter (and because [we can](https://www.youtube.com/watch?v=BVd-rYIqSy8), using the `PARAMETRIC` [keyword]((https://www.seamplex.com/wasora/reference.html#parametric)) in [Fino](https://www.seamplex.com/fino)).
+     So we can study how the linearised stresses change with the inclusion of a branch of a certain diameter (and because [we can](https://www.youtube.com/watch?v=BVd-rYIqSy8), using the `PARAMETRIC` [keyword](https://www.seamplex.com/wasora/reference.html#parametric) in [Fino](https://www.seamplex.com/fino)).
 
  2. Why do you want to know how the stresses change with the inclusion of a branch?  
      Because, if we have the time, it is worth to do something harder than originally asked (football joke: it is like training throw-ins with watermelons during the week so you can reach the area with a regular ball on the weekends).  
 
 So here we are again with the case study where we have to compute the linearised stresses at certain SCLs located near the interface between two different kinds of steels during operational and incidental transients of the plant. The first part is then to “bake” the pipes, modeling a thermal transient with time-dependent temperature or convection (it depends on the system) boundary conditions. This steps gives a time and space-dependent temperature\ $T(x,y,z,t)$.
 
-Then we proceed to “shake” the pipes, i.e. to compute the natural frequencies and associated oscillations modes. Employing the floor response spectra and combining the individual contributions with the SRSS method discussed in section\ [@sec:seismic], we obtain a distributed load vector\ $\vec{f}(x,y,z)$ which is statically equivalent to the design earthquake.
+Then we proceed to “shake” the pipes, i.e. to compute the natural frequencies and associated oscillations modes. Employing the floor response spectra and combining the individual contributions with the SRSS method discussed in section\ [@sec:seismic], we obtain a distributed load vector\ $\mathbf{f}(x,y,z)$ which is statically equivalent to the design earthquake.
 
 Finally we attempt to “break” the pipes successively solving many steady-state elastic problems for different times\ $t$ of the operational transient. Some remarks about this step:
 
  3. The temperature distribution\ $T(x,y,t,z)$ for bullets 1 & 2 is the generalisation of the reduced-model procedure explained in\ [@sec:thermal].
  4. The internal faces of the pipes are subject to an uniform pressure\ $p(t)$ given by the definition of the transient
  5. There are mechanical supports throughout the pipe system. Depending on the type of the support (i.e. vertical, lateral, axial, full, etc.) one or more degrees of freedom (i.e. $u$, $v$ and/or $w$) are fixed to zero. The ends of the CAD models were chosen always to have axially-null displacements.
- 6. The earthquake-equivalent volumetric force\ $\vec{f}(x,y,z)$ should only be applied at the time\ $t$ where the maximum stresses occur. Note that due to the discussion from\ [@sec:linearity] we cannot compute the stresses raised just by\ $\vec{f}(x,y,z)$ and then add them to the main stresses. The force has to be included into the “break” step. An educated guess of the time where the maximum stress occur is usually enough. Anyway, it might be necessary a trial and error scheme to find the sweet spot.
- 7. According to ASME\ III, the seismic load has to be applied during two seconds with the two possible signs. That is to say, apply $\vec{f}(x,y,z)$ during two seconds and then  $-\vec{f}(x,y,z)$ during two further seconds when the main stresses are maximums.
+ 6. The earthquake-equivalent volumetric force\ $\mathbf{f}(x,y,z)$ should only be applied at the time\ $t$ where the maximum stresses occur. Note that due to the discussion from\ [@sec:linearity] we cannot compute the stresses raised just by\ $\mathbf{f}(x,y,z)$ and then add them to the main stresses. The force has to be included into the “break” step. An educated guess of the time where the maximum stress occur is usually enough. Anyway, it might be necessary a trial and error scheme to find the sweet spot.
+ 7. According to ASME\ III, the seismic load has to be applied during two seconds with the two possible signs. That is to say, apply $\mathbf{f}(x,y,z)$ during two seconds and then  $-\mathbf{f}(x,y,z)$ during two further seconds when the main stresses are maximums.
  8. A number of stress classification lines have to be defined. The locations were previously accorded with the plant owner and/or the regulator.
  9. The stress linearisation has to be performed individually for each principal stress\ $\sigma_1$, $\sigma_2$ and $\sigma_3$ to fulfill the requirements ASME\ III\ NB-3126 (see [@sec:in-air] below).  
  10. This “break” step is linear.
 
 $$ \sigma = E(T) \cdot \epsilon $$
 
-What changes with temperature is the slope of\ $\sigma$ with respect to\ $\epsilon$ (think and imagine!). On the other hand, we have a given non-trivial temperature distribution\ $T(\vec{x}, t)$ within the pipes that is a snapshot of a transient heat conduction problem at a certain time\ $t$ (think and picture yourself taking photos of the temperature distribution changing in time). Let us now forget about the time, as after all we are solving a steady-state elastic problem. Now you can trust me or ask a FEM teacher, but the continuous displacement formulation can be loosely written as
+What changes with temperature is the slope of\ $\sigma$ with respect to\ $\epsilon$ (think and imagine!). On the other hand, we have a given non-trivial temperature distribution\ $T(\mathbf{x}, t)$ within the pipes that is a snapshot of a transient heat conduction problem at a certain time\ $t$ (think and picture yourself taking photos of the temperature distribution changing in time). Let us now forget about the time, as after all we are solving a steady-state elastic problem. Now you can trust me or ask a FEM teacher, but the continuous displacement formulation can be loosely written as
 
-$$ K\big[E\left(T(\vec{x})\right), \vec{x}\big] \cdot \vec{u}(\vec{x}) = \vec{b}(\vec{x})$$
+$$ K\big[E\left(T(\mathbf{x})\right), \mathbf{x}\big] \cdot \mathbf{u}(\mathbf{x}) = \mathbf{b}(\mathbf{x})$$
 
-If you notice, the complex dependence of the stiffnes matrix\ $K$ can be reduced to just the spatial vector\ $\vec{x}$:
+If you notice, the complex dependence of the stiffnes matrix\ $K$ can be reduced to just the spatial vector\ $\mathbf{x}$:
 
-$$ K(\vec{x}) \cdot \vec{u}(\vec{x}) = \vec{b}(\vec{x})$$
+$$ K(\mathbf{x}) \cdot \mathbf{u}(\mathbf{x}) = \mathbf{b}(\mathbf{x})$$
 
-And this last equation is linear in\ $\vec{u}$. In effect, the discretization step means to integrate over\ $\vec{x}$. As\ $K$, $\vec{u}$ and\ $\vec{b}$ depend only on\ $\vec{x}$, then after integration one gets just numbers with the matrix representation of\ [@eq:kub]. Again, you can either trust me, ask a teacher or go through with the maths.
+And this last equation is linear in\ $\mathbf{u}$. In effect, the discretization step means to integrate over\ $\mathbf{x}$. As\ $K$, $\mathbf{u}$ and\ $\mathbf{b}$ depend only on\ $\mathbf{x}$, then after integration one gets just numbers with the matrix representation of\ [@eq:kub]. Again, you can either trust me, ask a teacher or go through with the maths.
 
 
 ::::: {#fig:case-cad}
  * take into account that even within the finite element method, there is a wide variety of complexity in the problems that can be solved
  * follow the “five whys rule” before compute anything, probably you do not need to
  * use engineering judgment and make sure understand the [“wronger than wrong”](https://en.wikipedia.org/wiki/Wronger_than_wrong) concept
- * play with your favourite FEM solver (mine is [CAEplex](https://caeplex.com>)) solving simple cases, trying to predict the results and picturing the stress tensor and its eigenvectors in your imagination
+ * play with your favourite FEM solver (mine is [CAEplex](https://caeplex.com)) solving simple cases, trying to predict the results and picturing the stress tensor and its eigenvectors in your imagination
  * grab any stress distribution from any of your FEM projects, compute the iso-stress curves and the draw normal lines to them to get acquainted with SCLs
  * first search online for “stress linearisation” (or “linearization” if you want) and then get a copy of ASME\ VIII Div\ 2 Annex 5-A
  * try to show that even though the principal stresses are not linear with respect to summation, they are linear with respect to multiplication (with pencil and paper)
  * clone the [parametric tee repository](https://bitbucket.org/seamplex/tee), understand how the figures from\ [@sec:parametric] were built and expand them to cover “we might go on...” bullets
  * think thermal-mechanical plus earthquakes as “bake, break and shake” problems
  * the elastic problem is still linear after all
+
+
+# Referenced works
+
+## Books
+
+## Articles
+
+## Reports
+
+## Blog posts
+
+## Websites
+
+esyscmd([[lynx -dump -listonly nafems4.html | grep http | awk '{print $2}' | sort -u | uniq | grep -v wikipedia | awk '{printf(" * <%s>\n", $1)}']])
+
+## Wikipedia
+
+esyscmd([[lynx -dump -listonly nafems4.html | grep http | awk '{print $2}' | sort -u | uniq | grep wikipedia | awk '{printf(" * <%s>\n", $1)}']])

simple.html

 $for(include-before)$
 $include-before$
 $endfor$
+<!-- <span class="math display">\renewcommand\vec[1]{\mathbf{#1}}</span> -->
 <div class="container">
  <div class="row">
   <div class="col text-center">

simple.tex

 
 \include{syntax}
 
+\renewcommand\vec[1]{\mathbf{#1}}
+
 \begin{document}
 % \thispagestyle{empty}
 % \thispagestyle{plain}