## Materials and Methods
### Macroscopic Scale Model
\right\} \tag{12}\]  
Similarly, the elemental decomposition of \(\chi\) would be  
\[\chi_{f}=\sum\chi_{f}^{e}\text{ and }\chi_{s}=\sum\chi_{s}^{e} \tag{13}\]  
considering \(\chi^{e}\) as the contribution of element \(e\) to the total integration. Interestingly, the relation between the quadratic function of the nodal values for each element would be linear:  
\[\left\{ \begin{array}{l}\frac{\partial\chi_{f}}{\partial\phi_{i}}\\ \frac{\partial\chi_{f}}{\partial\phi_{i}}\\ \frac{\partial\chi_{f}}{\partial\phi_{i}}\\ \frac{\partial\chi_{f}}{\partial\phi_{m}}\end{array} \right\}=\frac{\partial \chi_{f}^{e}}{\partial\{\phi\}^{e}}=S^{e}\{\phi\}^{e} \tag{14}\]  
where the matrix \(S\) depended on the shape function \(N\), which in fact was associated with the nodal elements as \(N=[N_{i},N_{j},N_{l},N_{m}]\). The matrix \(S^{e}\) was symmetric and could be represented as  
\[S^{e}_{rs}=\left(\int\frac{\partial N_{r}}{\partial x}\frac{\partial N_{s}}{ \partial x}d\Omega+\int\frac{\partial N_{r}}{\partial y}\frac{\partial N_{s}} {\partial y}d\Omega+\int\frac{\partial N_{r}}{\partial z}\frac{\partial N_{s} }{\partial z}d\Omega\right)D \tag{15}\]  
The term \(\chi_{s}\) was defined by  
\[\chi_{s}=-\int\frac{1}{2}\mu_{a}\phi^{2}d\Omega \tag{16}\]  
Equation (16) could be redefined for other nodes of any element as  
\[\frac{\partial\chi_{s}^{e}}{\partial\phi_{i}}=-\int\mu_{a}N_{i}\big{[}N_{i}, N_{j},\ N_{m},\ N_{e}\big{]}d\Omega\{\phi^{e}\} \tag{17}\]  
Using the matrix form, it would be possible to state the total contribution for a 4-noded element as  
\[\left\{ \begin{array}{l}\frac{\partial\chi_{s}}{\partial\phi_{i}}\\ \frac{\partial\chi_{s}}{\partial\phi_{j}}\\ \frac{\partial\chi_{s}}{\partial\phi_{i}}\\ \frac{\partial\chi_{s}}{\partial\phi_{i}}\end{array} \right\}^{e}=-\int\mu_{a} \left[ \begin{array}{cccc}N_{i}^{2}&N_{i}N_{j}&N_{i}N_{m}&N_{i}N_{e}\\ N_{j}N_{i}&N_{j}^{2}&N_{j}N_{m}&N_{j}N_{e}\\ N_{m}N_{i}&N_{m}N_{j}&N_{m}^{2}&N_{m}N_{e}\\ N_{e}N_{i}&N_{e}N_{j}&N_{e}N_{m}&N_{e}^{2}\end{array}   
From Equations (12), (14), (19),  
\[S\{\phi\}=\mu_{a}M\{\phi\} \tag{21}\]  
The solution of Equation (21) gave the eigenvalues of \(\mu_{a}\). This could be further extended for any number of associated nodes. Here, 4-noded tetrahedron elements were used. The mathematical description for a 4-noded tetrahedron element is given in Appendix A. The seamless mesh containing the tetrahedron elements is shown in Figure 3.  
The employed mesh consisted of 468,345 domain elements including 34,056 boundary elements and 820 edge elements. The term \(D\) could be expressed in terms of absorption (\(\mu_{a}\)) and reduced scattering (\(\mu_{s}\)) coefficient as  
\[D=\frac{1}{3(\mu_{a}+\mu_{s})} \tag{22}\]  
The present model was implemented and solved using the MUMPS solver, which allowed better exploitation of computer resources and reduced numerical fluctuations associated with the FEM technique. The solver configurations included: (1) a memory allocation factor of 1.2; (b) pivoting with a pivot threshold of 0.1; and (c) hyper-threading enabled for the multiprocessor cluster. Due to the employment of a fine mesh, the numerical simulations were executed on supercomputer with Intel(r) Xeon E5-2670 2.60 GHz using 32 physical cores and hyper-threaded to 64.