The assembling procedure from a local element matrix
to a global matrix
has the same routines for two- and three-dimensional structures. Therefore, the assembling procedure is demonstrated on a simpler two-dimensional example as shown in Fig. 5.3. The dimension of the local matrix
, where is the number of grid nodes on the finite element ( for triangles and for tetrahedrons) and is the number of unknown variables on a grid node. The dimension of the global matrix
, where is the total number of grid nodes in the discretized domain.
If it is assumed that there is only one sought variable
The local matrix uses the local node indexes which are 1, 2, and 3 for every finite element. The global indexes for these grid nodes are different. There must exists a transformation which projects the local indexes and of the components to the global indexes and . For example, the element 66 in Fig. 5.3 with its local nodes 1, 2, and 3 has the global nodes 42, 30, and 25 and the index transformation is
Through the Dirichlet boundary conditions the values on the surface grid nodes are already fixed with the so-called Dirichlet value. Therefore, it is not necessary and even not allowed to recalculate the values on these grid nodes from the global equation system , because it is impossible to obtain the same Dirichlet values by solving the equation system. These surface grid nodes must be treated differently with the Dirichlet value. If on the global node there is a Dirichlet boundary value , the global equation system must be changed to
In practice it is more comfortable to multiplicate with a transformation matrix
The simulated structures generally consist of several segments. Normally a segment is a continuous region of one material and so there are everywhere the same material characteristics. Therefore, the electrical and mechanical behavior can be described with the same parameters within a segment.
The used meshing module makes a separated grid for every segment. This means that for a global grid point at the interface between two different segments there exist two different indices, because each segment has its own global index, as shown in Fig. 5.4.
The calculated displacements would be wrong, because the entries from segment A for this interface points in the global stiffness matrix do not take into account that there is also a stiffness from the other segment B due to the separate assembling of the segments. The same but reverse explanation is valid for the entries from segment B. Therefore, the stiffness in matrix and also the force in must be corrected in the way that the entries with second index are added to the first index of the grid point and then all entries with second index are set to 0.
In the following the correcting procedure in is demonstrated for the representative interface point with Index 35 in Segment A and Index 85 in Segment B (see Fig. 5.4).
After solving the linear system for displacements (5.102), has the correct value and . In order to also have the correct displacement for the second Index 85 (segment B), must be set to .
The oxidation problem induces 5 unknown variables: the oxidation concentration , the normalized silicon concentration , and the three components , , and (in x-, y-, and z-direction) of the displacement vector . If the dopant redistribution is taken into account, than the unknown variables for the concentrations of the different species have also to be considered. In case the five-stream diffusion model these are the 5 variables for , , , , and (see Section 4.2).
The following explanation is focused on the pure oxidation problem. Therefore, the discretization of the simulation domain with grid nodes leads to 5 variables on each grid node. After assembling of the whole equation system there are in total 5 unknown variables and, as already described in Section 5.3.1, the dimension of the global matrix is 55.
As shown in Fig. 5.5 the global matrix consists of the two sub-matrices and . contains the coupled entries from the oxidant diffusion (5.65) and the -dynamics (5.74). Because of the unknown variables and the size of is 22. is the global stiffness matrix for the mechanics with the three displacement components , , and , and so its dimension is 33.
The structure of the local system matrix for one finite element is the same as for the global one (see Fig. 5.5). The only difference is that because of the 4 nodes of a tetrahedron. So the dimension of the local matrix is 2020 with the 88 and 1212 sub-matrices.
The first part of the equation system which describes the oxidation diffusion and the -dynamics, is non-linear because of the coupling between and in the form (see (5.65) and (5.74)). The variables in the non-linear system can not be calculated directly, only their increments and . The second part of the system which is responsible for the mechanical problem is linear, as marked in Fig. 5.5.
For a quasi-stationary time step there is no coupling between the --system and the mechanical system, because the equations for and are not functions of displacements. Although the normalized additional volume after a time step is calculated with and , the mechanical equations for the displacements also are not functions of and . Because of the missing coupling the off-diagonal sub-matrices in the global equation system are zero (see Fig. 5.5).
On the right-hand side of the non-linear subsystem are the residuals calculated with the results from the last Newton iteration as explained in the next section. The right-hand side of the linear subsystem contains the (internal) forces on the grid nodes.
In contrast to the linear mechanical subsystem, where the displacements can be calculated directly with one of the various standard methods like Gaussian elimination, the solving of the non-linear part demands other routines like the used Newton Method .
The global non-linear subsystem can be written in the form
With the Newton formula the solution vector for the actual time step becomes
Transforming the Newton formula to the form
Because the Newton method is only a first order approximation of the solution, an iteration never provides the exact results. Therefore the right-hand side of the non-linear system (5.103) can not be 0, instead there always exist residuals .
The quality of the approximation is increased with each iteration, but the number of iterations must be limited with termination conditions. With these conditions also the accuracy of the approximation can be controlled. It makes sense to use the following termination conditions 
For one finite tetrahedral element with its 4 nodes () the non-linear system (5.103) is