{'fraction_non_alphanumeric': 0.02995371061669404, 'fraction_numerical': 0.011617142003882336, 'mean_word_length': 5.004303388918774, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A Neural Network-enhanced Reproducing Kernel Particle Method for Modeling Strain Localization

Modeling the localized intensive deformation in a damaged solid requires highly refined discretization for accurate prediction, which significantly increases the computational cost. Although adaptive model refinement can be employed for enhanced effectiveness, it is cumbersome for the traditional mesh-based methods to perform while modeling the evolving localizations. In this work, neural network-enhanced reproducing kernel particle method (NN-RKPM) is proposed, where the location, orientation, and shape of the solution transition near a localization is automatically captured by the NN approximation via a block-level neural network optimization. The weights and biases in the blocked parametrization network control the location and orientation of the localization. The designed basic four-kernel NN block is capable of capturing a triple junction or a quadruple junction topological pattern, while more complicated localization topological patters are captured by the superposition of multiple four-kernel NN blocks. The standard RK approximation is then utilized to approximate the smooth part of the solution, which permits a much coarser discretization than the high-resolution discretization needed to capture sharp solution transitions with the conventional methods. A regularization of the neural network approximation is additionally introduced for discretization-independent material responses. The effectiveness of the proposed NN-RKPM is verified by a series of numerical verifications.


Introduction
Modeling many complicated and evolving localizations, such as, fracture and shear bands, is a challenging task. In fracture modeling, significant developments have been made over the past several decades, such as discrete crack approaches and diffuse crack approaches. The discrete crack approaches include meshfree methods with visibility/diffraction criterion 1,2 , local partition of unity-based enrichment 3,4 , and extended and generalized finite element methods [5][6][7] .
This class of methods can yield mesh-independent solutions but requires surface tracking, which becomes ineffective when the crack path is highly complicated. The diffuse crack approaches offer the possibility of modeling complicated crack patterns without surface tracking. However, the solution of this method-type is prone to being discretization-sensitive. If the discretization sensitivity is not properly addressed, numerical results including energy, force, and damage patterns can be strongly affected by the choice of discretization due to the loss of ellipticity in the mathematical formulation 8 . For the diffuse crack approaches, nonlocal models with intrinsic length scales, such as quantity averaging 9 , high order gradient models [10][11][12] , and phase field methods [13][14][15][16] have successfully remedied the discretization sensitivity. However, nonlocal methods often require very fine discretization, which significantly increases the computational cost. While adaptive model refinement can be employed for enhanced effectiveness, it is cumbersome for the traditional mesh-based methods to perform.
The deep artificial neural network has demonstrated its ability in information processing, such as image recognition, and is rapidly extending its regime to various fields, including mechanistic machine learning. Its successful widespread adoption is mainly attributed to the adaptive nature of manipulating the function space based on data via minimization of the loss function. Such flexibility allows a computational model to capture certain information that is hard to identify by conventional approaches. For example,   17 introduced the deep autoencoder to discover the latent low-dimensional data pattern embedded in noisy material datasets, which significantly enhanced the effectiveness of the physics-constrained data-driven approach developed by He and Chen (2020) 18 19 . Lee and Carlberg (2020) 20  Recently, the utilization of neural networks (NNs) for solving partial differential equations (PDE) in physics and mechanics has gained increased attention. This class of studies explores the NN's ability to be used as an approximation, considering its flexible function space construction. Raissi et al. (2019) 21 proposed a physics-informed neural network (PINN), which serves as an approximation of the solution of a PDE in the framework of a collocation method, and showed the effectiveness of the method on the one-dimensional Burgers equation using densely-connected, deep neural networks with the hyperbolic tangent activation function.
Haghighat and Juanes (2021) 22 introduced a Python package SciANN for scientific computing using PINN and showed that the method captured strain and stress localization produced in a perfectly plastic material. However, they used 4 densely-connected layers with 100 neurons per layer, which involves 100 million unknowns, to obtain the results. Samaniego et al. (2020) 23 showed that the potential-based loss functions lead to better results with significantly less unknowns than the collocation-based loss function that is widely used in PINN. However, the aforementioned studies employed multiple fully-connected hidden layers, which may not be an optimal choice to the solution of PDEs. Zhang et al. (2020) 24 devised the deep neural network to represent standard approximations such as finite elements and reproducing kernel particles, and treated the nodal positions as unknown network parameters. This element-or particle-wise network is a sparse network, indicating a potential of designing a more efficient network.
Although NN for solving PDEs in physics and mechanics has drawn increased attention, the research on development of NN for modeling discontinuities and localizations remains very limited and deserves further investigation.
While a neural network-based approximation can be highly effective in capturing localizations due to the adaptive nature of the function space, it may be inefficient to employ the NN-based approximation across the entire domain, in which the vast majority exhibits a smooth solution that can be accurately approximated by traditional methods, such as meshfree methods 25,26 and isogeometric analysis 27 . Particularly, in the reproducing kernel (RK) approximation 25,26 , the order of continuity of the approximation and the order of basis are independently controlled, and the RK approximation has been proven to effectively and accurately capture smooth solutions with a coarse domain discretization. In this work, the RK approximation will be utilized to approximate the smooth part of the solutions while a NN approximation will be NN designed to effectively capture sharp solution transitions with proper regularization for dissipation energy objectivity.
The paper is organized as follows. In Section 2, the basic equations will be provided, including the minimization problems under consideration and the reproducing kernel approximation. In Section 3, a neural network-enhanced reproducing kernel approximation is presented. In Section 4, the implementation details including the neural network structure and solution procedure are provided. This is followed by numerical examples in Section 5 and concluding remarks in Section 6.

Minimization problems under consideration
Let Ω ∈ ℝ be a domain with the space dimension and Ω = ∂Ω ∪ ∂Ω h be the boundary of Ω with the Dirichlet boundary ∂Ω and the Neumann boundary ∂Ω h . In this work, the following minimization problem is considered: where is the displacement vector, ( ) is the strain energy density functional, and ( ) is the external energy functional. The Dirichlet boundary condition = ̅ on Ω g is imposed by the penalty method with the penalty parameter . For linear elasticity, ( ) = ( ) and ( ) are taken to have the following form: where , , , and are the strain tensor, elasticity tensor, body force vector, and boundary traction vector, respectively. For tension-compression decomposed elasticity-damage problems, ( ) = ( ) are defined as follows 13 : where the damage parameter, , degradation function, ( ), and dissipation energy density functional, ̅ , are introduced, and 0 + and 0 − are tensile and compressive strain energies, respectively, denoted as with the principal strain tensor, ̅, Lamé's first and second parameters, and , and Macaulay brackets, 〈⋅⟩ = max(⋅ ,0).
In this work, the following damage model is employed: where , 0 , and are irreversible internal state variable, limit elastic strain, and critical equivalent strain, respectively, with satisfing the following Kuhn-Tucker conditions: where ̅ is defined as Note that, in one-dimensional uniaxial tension, ̅ = 11 . For this considered damage model, ̅ in (3) is defined as follows: where and are as follows: The stress defined as The dissipation functional defined in Eq. (10) is derived such that, for the adopted damage law in Eq. (5)- (7), the strain energy density functional satisfies = ∫ : , which leads to Note that with the dissipation functional (10), the minimization problem in (1) is consistent with the standard variational formulation shown below: In this work, Π( ) in Eq. (1) is taken as the loss function for the numerical procedures in the neural network framework. The dissipation functional ̅ introduces additional energy to the elastic energy, ( ( ( ))) 0 + ( ) + 0 − ( ), such that the stationary condition of the potential energy functional Π( ) = 0 yields the variational equation in (13) with the desired stress degradation in (12).

Reproducing kernel approximation
In this work, the standard reproducing kernel (RK) approximation 25,26 for solving the minimization problem will be enriched with a neural network (NN) approximation (see Section 3) to better capture the strain localization. In the proposed approach, the RK approximation is used to approximate the smooth part of the solution while the NN approximation is used to approximate the sharp solution transition near localization. In this subsection, the RK approximation is introduced as the basis for the smooth part of the solution; the NN enhancement will be discussed in the next section.
Let us consider a domain Ω discretized by nodes, as shown in Figure 1, with nodal coordinate with 1 ≤ ≤ . The RK approximation, ℎ ( ), of a function ( ) is where Ψ ( ) is RK shape function of node and is the generalized nodal coefficient of node . The RK shape function, Ψ ( ), is a correction of a kernel function, Φ ( − ), defined on the compact support of node with a support size of : where ( − ) is a basis function, = ( 1 , 2 , … , ) with dimension , and | | ≡ ∑ =1 .
is defined as The coefficients, ( ) for | | ≤ , are the solutions of the following set of reproducing conditions: which leads to the explicit form of Ψ ( ) as follows.
where ( ) is the moment matrix and ( − ) is the basis vector defined as The order of continuity of the RK approximation is determined by the kernel function, Φ ( − ), while the polynomial completeness of the approximation is determined by the basis vector, ( − ). Thus, it is straightforward to introduce high order continuity into the approximation space, which makes the RK approximation more appealing for approximating the smooth part of solution than the C 0 interpolation-type approximations used in finite element methods.

Neural Network-enhanced Reproducing Kernel Approximation
In this chapter, a neural network-enhanced reproducing kernel (NN-RK) approximation is proposed. Neural network (NN) architectures are developed with interpretable weights and biases to automatically detect the locations and orientations of strain localization, and to construct functions with sharp transitions near the localizations. This section presents the details of the network structure that incorporates the RK approximation.
where ( ) and ( ) are the RK approximation and NN approximation, respectively. In the above approximation, the NN approximation is constructed to detect localization and to introduce a fine-scale feature to the solution in a region close to a localization, while the RK approximation is employed to capture the overall smooth behavior of the solution. The RK approximation takes its traditional form discussed in Section 2: where Ψ ( ) and are RK shape function and RK generalized coefficient associated with RK node , respectively.

Considerations on the construction of neural network approximation
In the construction of the NN approximation, we focus on the following considerations: 1) the sharp solution transition associated with highly localized strain with weak discontinuity be approximated by NN approximation, 2) the position and orientation of complex localization paths be automatically captured, 3) the NN approximation influences only local regions close to localizations.
In Consideration 1, the approximation of weak discontinuities in the displacement field on the interface between damaged and undamaged materials is considered, as illustrated in Figure 2. To achieve Consideration 2, unknown parameters representing the localization position and orientation will be treated as unknowns and be included in the NN approximation. Also, a blocklevel NN approximation will be constructed to capture complex localization topology by a superposition of multiple block-level NN approximations, each of which are constructed to approximate relatively simple localization topology. Consideration 3 is related to the computational efficiency. As the smooth solution away from the localization can be efficiently obtained by the RK approximation, the domain of influence of the NN approximation is controlled to be in the regions near localizations. The NN approximation will be designed such that the domain of influence is controlled by adjustable NN parameters that are automatically determined by the optimization procedures. For example, when multiple sparsely-distributed localization clusters are populated in the domain, as shown in Figure 3, the domain of influence in the NN approximation is sparsely distributed near the localization clusters. In the construction of the NN approximation space that satisfies Considerations 1-3, parameters that control the location, orientation, and shape of a localization are included, and these parameters are determined by the potential energy minimization. In the next subsection, a design of a blocklevel NN approximation satisfying Considerations 1-3 is proposed.

Block-level approximation
The proposed NN approximation is constructed by the block-level NN approximations as follows.
where is -th block-level approximation with unknown weights and biases lumped into , and is the number of blocks. Each NN block is constructed to capture a certain portion of the localizations distributed throughout the domain and to enrich the solution near localizations.
Here, the introduction of NN blocks is intended to reduce the number of unknown weights and biases in representing a localization network. Also, the superposition of multiple block-level approximations is performed to capture complex localization topology.
In this work, the block-level approximation takes the following form: where is a monomial function and ̂ is the -th normalized NN kernel function of theth NN block, both are to be constructed automatically with the neural networks. The role of the kernel function ̂ is to capture the location, orientation, and the transition shape of the localization. Here, ̂ involves two sets of unknown parameters: the localization location and orientation parameter set, ≡ { } =1 , and the kernel shape parameter set, ≡ The other NN parameter sets, and , for enhancement of the RK approximation for capturing localizations will be introduced next.

NN kernel functions by
Let us consider the following one-dimensional kernel function for the NN approximation that is constructed by the multiplication of two regularized step functions ̅ 1 and ̅ 2 (see Figure 4): for -th NN kernel function of block , where is the 1-D coordinate and is a set of shape control parameters to be determined by the NN optimization. For the remaining section, the subscript will denote that the corresponding quantity is associated with -th regularized step function, ̅ .
where and ( ; ) are the normalized coordinate and Softplus function, respectively. The Softplus function ( ; ) is defined as where controls the sharpness in the transition of derivative. Given values of , ( ; ) and ̅ ( ; ) are plotted in Figure 5 where larger values lead to sharper derivative transition.
where ̂( ) is the normalized NN kernel and is a small number introduced to prevent sparsely distributed NN kernel functions from having too large domains of influence (see Figure   9), so that the NN approximation only influences regions near localizations. must be small enough to approximate the partition of unity property for the region where the material point is covered by an enough number of NN kernel functions. In this work, = 10 −3 is used.

Remarks 3.1:
1. Note that is not a trainable parameter; the value of does not evolve during the potential energy minimization.
2. The localization width and intensity, and the transition shape of the NN kernel function are controlled by the aforementioned shape control parameter that evolves in the potential energy minimization process: the localization width is controlled by ̅ 2 − ̅ 1 as shown in Figure   8, the localization intensity is controlled by as shown in Figure 6, and the localization transition shape is controlled by as shown in Figure 7. The extension of the NN kernel to two-dimension with coordinate = ( 1 , 2 ) is straightforward: for the -th kernel in the -th block, } is the kernel shape parameter set associated for the -th direction and with -th regularized step function. Similar to Eq. (26), the regularized step function, ̅ , takes the following form where the superscript, , is dropped for brevity: where ̅ is defined in Eq. (26), and = (−1) ( − ̅ )/ is similarly introduced.

Parametrization by
In multi-dimensional problems, highly complex localization topological patterns can be present.
In this approach, the complicated localization patterns will be projected onto a parametric space in each NN block "A" with parametric coordinates, Α : → Α , where ∈ ℝ and Α ∈ ℝ are physical coordinates and parametric coordinates of block Α, respectively. Figure

Approximation ability of a single neural network block
As discussed in Section 3.2.2, the superposition of block-level NN parametrizations is introduced to capture the topological localization patterns. For two-dimensional problems, each NN block contains four NN kernels defined in (29)-(30) along with the parameterization network defined in (31)- (33). Figure 11 shows the topological geometries that can be captured by a single four-kernel NN block. The approximation ℎ is obtained by the standard 2 minimization to obtain NN parameters = { L , S , P }.
The domain is discretized using 10x10 uniformly spaced RK nodes and 72x72 Gauss cells with 2x2 Gauss integration for domain integration in (34). For all convergence studies discussed in this section, the 2 error norm is integrated using 8x8 Gauss integration and the same 72x72 Gauss cells. The convergence rates are calculated with respect to the square root of the varied NN variables.
For the approximation of TF1, is varied as 8, 16, and 32. The observed 2 error norm decreases as the number of neurons is increased at an average rate of convergence of 2.282, as seen in Figure 13. convergence rate for the first four cases is 0.430, but the error reduces significantly as is increased to 128, leading to an average convergence rate of 1.620 with all 5 cases considered.
TF2 presents a more challenging function to approximate than TF1, as the neural network must learn the transition zone location as well as the varying magnitude of the jump. As such, a more refined neural network is needed in the parametrization layer. Once a sufficiently large number of neurons is incorporated, the error reduces significantly. In the approximation of TF3, is varied across 1, 2, 4, and 8. The average rate of convergence is shown in Figure 15 as 3.578. The 2 error norm monotonically decreases as is increased, as shown in Figure 15. The 2 error decreases most significantly between the 2-and 4-block cases. Due to the existence of three triple-junctions in TF3, increasing beyond 4 appears to be less effective.

Regularization
To ensure the proposed NN-RK approximation is discretization-insensitive when applied to damage mechanics, a regularization is introduced into the NN kernel function. Although the NN kernel contains the sharpness parameter in Eqs. (26) and (30), the localization width is also affected by the parametrization, more specifically, the gradient of the parametric coordinate with respect to the physical coordinate . In this subsection, the parameters affecting the sharpness of kernel function in the parameterization are to be properly scaled by modifying Eqs. (26) and (30), so that the localization width is solely determined by the parameter . This introduces a length scale to the physical parameters that can then be used to define an objective energy dissipation.
To analyze the influence of the parametrization on the transition sharpness, consider the first order Taylor expansion of ( ) as Let us define ̂ as where ‖⋅‖ is the Euclidean norm. Introducing (35) into the normalized coordinate in the regularized step function ̅ in (26) and (30), where ̅ = (̅). Note that ̂| Here, the sharpness in the NN kernel function ̅ is entirely controlled by . With this modification, imposing a numerical length scale ℓ as a lower bound of regularizes the solution in problems with strain localization.
To obtain an objective energy dissipation, the relation among the physical parameters introduced in (6) the map ℝ → ℝ with input ∈ ℝ and output ℎ ∈ ℝ. As the RK approximation is used in approximating the smooth part of the solution, a relatively coarse, predetermined discretization is sufficient. Since in this approach the RK shape functions are designed not to evolve throughout the simulation, the shape functions {Ψ ( )} =1 are precomputed and directly entered into the network as the input as shown in Figure 17. The RK network has the RK generalized coordinates as its weights. On the other hand, the NN approximation directly takes the coordinate as its input. Instead of employing a densely connected deep neural network, multiple network blocks are constructed in parallel, as discussed in Section 0, which introduces sparsity to the system of equations associated with the weights obtained by an optimization solver. Figure 17. Entire network structure that incorporates the RK approximation and the NN approximation. The unknown parameters associated with each part of the network are shown in red. The thick black arrows denote that the network connections entail unknown parameters and the thin grey arrows denote that the weight of the associated connection is one (1) and do not change throughout the computation.
The network structure of a block is described in Figure 18. All the unknown parameters , , and introduced in Section 0 are involved in each sub-block and solved by the minimization procedure. Here the parameterization network involves the localization location and orientation parameters set , with it outputs passing into the sub-block for NN kernel to construct the localization shape parameter set ,while the polynomial basis set sub-block is independently constructed. The details of the NN kernel sub-block are given in Figure 19 where the localization shape parameters is included. In order to regularize the NN approximation as discussed in Section 3.3, the built-in function of Tensorflow tf.keras.constraints.MinMaxNorm is utilized to impose the lower bound of . The NN operation details of polynomial sub-block to obtain the polynomial basis functions are presented in Figure 20.
The neuron-wise multiplication of the polynomial layer and kernel function layer is performed by employing tensorflow.keras.layers.multiply available in TensorFlow. Figure 21 shows how the ℎ derivatives are computed for use in solving the potential energy minimization problem.
To compute the spatial derivatives of , the automatic differentiation function provided by TensorFlow is utilized. To compute the spatial derivatives of , the input of the RK network, Figure 21 is replaced by the pre-computed shape function derivatives Ψ , ( ).

Solution procedures
The optimization problem defined in (1)

A. RK Solution Stage
In the first stage, the problem is solved only for the RK coefficients to find the initial guess of the RK coefficients ̅ as follows: In this stage, the smooth part of the solution by RK approximation is obtained. Note that this leads to a standard Galerkin-based RKPM formulation and the problem can be solved by using a standard matrix solver. Since the RK approximation obtained in the first stage can oscillate around the localization due to the Gibbs-type phenomenon, the filtered RK coefficients are used as the initial guess for the next step as follows: where ̃ is the filtered RK coefficient at node I and Ψ ̅ ( ) is an RK filter 29 . As the purpose of filtering the RK coefficients is to provide a non-oscillatory initial guess to the next solution stage, many appropriate filters are possible. In this work, the RK shape function is used as the RK filter, that is, Ψ ̅ ( ) = Ψ ( ).

B. NN-RK Solution Stage
In the second stage, Eq.

One-dimensional elasticity with pre-degraded material
The following one-dimensional elastic bar problem is considered to investigate how the individual NN block approximations play a role in capturing localizations.
where , , ℎ, and are Young's modulus, body force, RK nodal spacing, and Dirichlet boundary value, respectively. The Dirichlet boundary values 1 = 0 and 2 = 0.5 are used. As shown in Figure 22, the material is pre-degraded locally with a small Young's modulus, and the bar is subjected to a smooth body force:   Figure 24 shows the total solution, NN approximation, and RK approximation. As expected, the NN approximation captures the very sharp solution transition while the RK approximation represents the overall smooth solution. As shown in Figure 25, the pure RK solution achieves a similar resolution with 801 unknowns while the transitions are insufficiently sharp even with 201 and 401 nodes, while the NN-enhanced approximation agrees very well with the exact solution with 57 unknowns, which is a 93% reduction in the number of unknowns. The number of unknowns required to capture localizations in multi-dimensional problems will be much more pronounced, and the proposed approach is expected to be even more effective compared to the standard RK approximation.

An elastic-damage bar under tension
To demonstrate the energy-based regularization effect along with the introduction of the NN length scale, a one-dimensional bar with imperfection under tension is considered. As shown in  Chen (2018) 28 . Contrarily, Figure 28 shows un-regularized load-displacement responses obtained in three cases with varying NN length scale proportional to the RK nodal spacing with the same damage model constants 0 and used. The localized strain field is shown in Figure   27 (b) where the strain distribution is converged as the domain is refined.
As shown in Figure 27 (c), the sudden jumps in the damage parameter, discussed in Section 3, are captured by the adaptive NN kernel functions. When smooth NN kernel functions with a fixed smoothness parameter = 5.0 (see Section 3.2.1) is used without NN optimization, severely oscillating stress field around the localization point is obtained as shown in Figure 29 (a) while the NN kernel with adaptive via NN optimization yields a significant improvement as shown in Figure 29 (b). Note that Wei and Chen (2018) 28 bypassed this issue by using a gradient smoothing technique implemented with smooth RK kernel functions to obtain stress fields without oscillation. This shows potential of implementing a gradient smoothing technique to integrate the NN-enhanced formulation with smooth NN kernel functions, which will be a task of future work.

Tensile specimen with asymmetric imperfection
As shown in Figure 30, a tensile specimen with an asymmetric imperfection is considered. The material properties used in Section 5. As shown in Figure 31, the proposed method yields consistent damage patterns on the four discretizations, which agree well with the damage pattern obtained by the implicit gradient RK regularization in Chen et al. (2000) 31 . The force-displacement curves are plotted in Figure 32 where all three models result in nearly identical reaction force profiles, which demonstrates the regularization effect of the proposed approach in multi-dimensional problems. Note that in the three employed RK discretization models, the RK nodal spacings are larger than the NN length scales, thus the sharp displacement transition due to the strain localization is entirely captured by the NN approximation (see Figure 33).

Pre-notched specimen under simple shear
To demonstrate the capability of the proposed method to model evolving localization paths, a pre-notched specimen under simple shear, as shown in Figure 34, is considered. In this work, the following damage model is applied:  The total CPU time of 572 minutes was used by the proposed method with a single GPU, NVIDIA A100 with 40 GB memory. In literature 32 , it was reported that CPU times of 1304 min and 957 min were required by a traditional phase field approach with 250×250 standard finite elements and a phase field method with an adaptive refinement, respectively. This comparison is encouraging although different program languages and hardware were employed. Also, it is worth noting that, in the proposed method, the NN approximation space is constructed based on a highly sparse neural network with interpretable NN weights and biases. As the computational efficiency that can be gained by utilizing sparse matrices becomes more significant with a larger problems size compared to using full matrices, it is expected that the CPU time saving achieved by the proposed method will also be more significant in solving larger scale problems.

L-Shaped Panel Test
To simulate relatively complex crack propagations, an L-shaped panel test 33 Figure 37(c). The numerically predicted damage fields are shown in Figure 38 where localization is initiated in a diagonal direction and changes its direction to the horizontal as the vertical deflection increases. As shown in Figure 38(d), the predicted damage path lies within the experimentally observed crack path range 33 . Figure 39 shows the loaddisplacement curve predicted by the proposed method with a good agreement with the experimental data 33 .

Conclusion
A neural network-enhanced reproducing kernel particle method has been proposed for modeling strain localization. The approximation space is constructed by enriching a standard Reproducing Kernel (RK) approximation with a neural network (NN) approximation. In the proposed function space, the NN approximation is constructed to capture the sharp solution transition near the region of strain localization while the RK approximation is utilized to approximate the smooth part of the solution. The capability of the RK approximation to accurately represent smooth functions with a coarse discretization allows the use of a small number of fixed RK shape functions for computational efficiency. In the NN approximation, adaptive NN kernel functions defined in a parametric coordinate introduce localized solution gradients with various types of topological patterns. The localization shape of the NN kernel functions is controlled by unknown shape parameters that are determined by the NN optimization. The parametrization is performed by the blocked parametrization network, with the weights of the network controlling the location and orientation of the localization, and they are also determined by the optimization.
One four-kernel NN block is capable of capturing a triple junction or a quadruple junction topological pattern. More complex topology can be captured by the superposition of multiple block-level NN approximations, each of which is defined in its own parametric domain. A series of function approximation testing problems showed that 1) the numerical solution converges with an increased number of neurons in the parameterization network and 2) the number of fourkernel NN blocks needed to capture complex localization topology with multiple triple junctions is greater than or equal to the number of the triple junctions.
To achieve a discretization-insensitive solution in modeling localization, a regularization of the NN-enhanced RK approximation has been introduced such that the localization width of NN kernels is determined by the NN length scale parameter, which leads to the dissipation energy objectivity. A set of numerical examples have been analyzed to demonstrate that the proposed method is effective in modeling evolving localization paths with discretization-insensitive solutions.

Appendix A. Target Function Construction
The expression used to construct target function TF3 used for the convergence study in Section 3.2.4 is provided. As shown in Figure