# Three Dimensional CS-FEM Phase-Field Modeling Technique for Brittle Fracture in Elastic Solids
## 2 Methods
### Three-Dimensional Cell Based Smoothed-Finite Elements  
In S-FEM, the idea of strain smoothing is combined with the standard underlying mesh of FEA by dividing each element into a number of sub-cells [1]. We use the eight-noded hexahedral elements with eight smoothing cells (Figure 2). They have already been proven to have a better accuracy than the standard eight-noded hexahedral element used in FEA with eight Gaussian points [38] using the same set of elements. Classical FEM involves isoparametric transformation for every element, making the problem very much mesh dependent because it involves inverting the Jacobian matrix. The solutions tend to vary with minor mesh distortions and for better modeling of certain topological features one needs mesh refinement or higher order elements, which makes the problem computationally expensive. The S-FEM can solve majority of computational problems using a base of linear triangular/tetrahedral mesh, which can be generated automatically. Since the shape functions are created on the basis of radial point interpolation, creating special elements to treat special topological features [18] is of minor hassle and does not need any transition or patch elements. In S-FEM, a smoothing operation is performed on the gradient of the displacement field for each smoothing cell in an element. Subsequently, the interior integration on each of the smoothing cell is transferred to the boundary surface area using the divergence theorem. Herein lies one of the biggest advantages of the S-FEM. This feature makes the solution pretty impervious to mesh distortion since we can avoid the isoparametric transformation and also contributes in saving computational cost. In this case, the gradient smoothing is the spatial average of the strain, over a smoothing hexahedral cell. A smoothing operation is performed to the gradient of displacement for each smoothing cell in an element:  
\[\
abla u^{h}(\mathbf{x}_{C})=\int_{\Omega}\
abla u^{h}(\mathbf{x})\Phi(\mathbf{x}-\mathbf{x}_{ C})d\mathbf{\Omega} \tag{18}\]  
where \(\mathbf{x}_{C}\) and \(\mathbf{\Omega}\) represent the center and volume of the smoothing domain, respectively, and \(\Phi\) is a smoothing function. Integration by parts on Equation (18) yields:  
\[\
abla u^{h}(\mathbf{x}_{C})=\int_{\Gamma}u^{h}(\mathbf{x})\mathbf{n}(\mathbf{x})\Phi(\mathbf{x}- \mathbf{x}_{C})d\mathbf{\Gamma}-\int_{\Omega}u^{h}(\mathbf{x})\
abla\Phi(\mathbf{x}-\mathbf{x }_{C})d\mathbf{\Omega} \tag{19}\]  
For simplicity, a piecewise constant smoothing function is applied, which is assumed to be constant within \(\Omega_{C}\) and vanish everywhere else,  
\[\Phi(\mathbf{x}-\mathbf{x}_{C})=\left\{ \begin{array}{ccc}\frac{1}{V_{ C}}&\text{for}&\mathbf{x}\in\mathbf{\Omega}_{c}\\ 0&\text{for}&\mathbf{x}\
otin\mathbf{\Omega}_{c}\end{array} \right. \tag{20}\]  
where \(V_{C}=\int_{\mathbf{\Omega}_{c}}^{d}\mathbf{\Omega}\) and \(\mathbf{\Omega}_{c}\) is the smoothing cell. Substituting Equation (20) into Equation (18), we can get the smoothed gradient or strain field  
\[\
abla u^{h}(\mathbf{x}_{C})=\int_{\Gamma_{C}}u^{h}(\mathbf{x})\mathbf{n}(\mathbf{x})\Phi(\bm {x}-\mathbf{x}_{C})d\mathbf{\Gamma}=\frac{1}{V_{c}}\int_{\Gamma_{c}}u^{h}(\mathbf{x}) \mathbf{n}(\mathbf{x})d\mathbf{\Gamma} \tag{21}\]  
with \(d\mathbf{\Gamma}_{c}\) denoting the boundary of the smoothing cell. It is noteworthy that the choice of the smoothing function makes the second term of Equation (19) vanish and the area integration is thus converted to line integration along the boundary of the smoothing cell. The smoothed strain can be obtained as \[\tilde{e}^{h}(\mathbf{x}_{C})=\sum_{I=1}^{n}\tilde{B}_{I}(\mathbf{x}_{C})\mathbf{d_{I}} \tag{22}\]  
where \(\tilde{B}_{I}\) is the smoothed strain gradient matrix given by  
\[\tilde{B}_{I}=\left[