ML-based identification of the interface regions for coupling local and nonlocal models
Noujoud Nader Patrick Diehl Marta D'Elia Christian Glusa Serge Prudhomme

Outage Probability Analysis of  Wireless Paths with Faulty Reconfigurable Intelligent Surfaces
    Mounir Bensalem and Admela Jukan
Technische Universität Braunschweig, Germany;
{mounir.bensalem,  a.jukan}@tu-bs.de







    Version of April 23, 2024
==============================================================================================================================


Local-nonlocal coupling approaches provide a means to combine the computational efficiency of local models and the accuracy of nonlocal models. However, the coupling process can be challenging, requiring expertise to identify the interface between local and nonlocal regions.
This study introduces a machine learning-based approach to automatically detect the regions in which the local and nonlocal models should be used in a coupling approach. 
This identification process takes as input the loading functions evaluated at the grid points and provides as output the selected model at those points. Training of the networks is based on datasets provided by classes of loading functions for which reference coupling configurations are computed using accurate coupled solutions, where accuracy is measured in terms of the relative error between the solution to the coupling approach and the solution to the nonlocal model. 
We study two approaches that differ from one another in terms of the data structure. The first approach, referred to as the full-domain input data approach, inputs the full load vector and outputs a full label vector. In this case, the classification process is carried out globally. The second approach consists of a window-based approach, where loads are preprocessed and partitioned into windows and the problem is formulated as a node-wise classification approach in which the central point of each window is treated individually.
The classification problems are solved via deep learning algorithms based on convolutional neural networks. The performance of these approaches is studied on one-dimensional numerical examples using F1-scores and accuracy metrics. In particular, it is shown that the windowing approach provides promising results, achieving an accuracy of 0.96 and an F1-score of 0.97. These results underscore the potential of the approach to automate coupling processes, leading to more accurate and computationally efficient solutions for material science applications.











§ INTRODUCTION


Recent advances in engineering mechanics have seen an increasing emphasis on coupling nonlocal models with classical local models, an approach driven by the need to address the excessive computational costs of high fidelity models, such as nonlocal equations, for modeling and simulation of complex behaviors in materials science. Nonlocal models, such as nonlocal diffusion and peridynamics, offer a more comprehensive representation of phenomena that are not adequately captured by classical Partial Differential Equations (PDEs). These phenomena include fracture and cracks in mechanics <cit.>. The scope of nonlocal models extends beyond diffusion and mechanics, finding applications in diverse areas like subsurface transport <cit.>, phase transitions <cit.> and image processing <cit.>.

However, the implementation and simulation of nonlocal models present several challenges. These include the complexity of handling nonlocal boundary conditions <cit.> and the high computational costs associated with their numerical solutions. To address these challenges, coupling between local and nonlocal models approach has been proposed. Coupling methods aim to merge the computational efficiency of PDEs with the accuracy of nonlocal models, particularly in cases where nonlocal effects are confined to specific domain areas, and the rest of the system can be effectively described by a PDE.

For a general overview of local-nonlocal coupling methods, classified as either force-based or energy-based formulations, we refer to the recent survey <cit.>.  Our focus is on the coupling of the bond-based peridynamic model with the classical linear elasticity continuum model, using a force-based coupling formulation. In this context, various coupling approaches have been proposed; these are based on matching the displacements (the Matching Displacement Coupling Method) or the stresses (the Matching Stress Coupling Method) either over an overlap region or a sharp interface obtained by shrinking the size of the nonlocal horizon when transitioning from a nonlocal to a local region. The latter approach is known as the variable horizon coupling method (VHCM); the interested reader can find an extended description in <cit.>.


A key aspect of the coupling process is that it requires domain expertise to effectively choose the interface between the local and nonlocal regions. In this work, we propose to circumvent this challenge by letting machine learning detect such an interface region using a supervised approach. Specifically, the use of machine learning involves training neural network architectures on datasets containing local-nonlocal coupling regions that yield a certain accuracy of the coupled solution. Once trained, the networks can automatically identify the regions of the computational domain where nonlocal effects are likely to be significant, solely based on external loading information. This approach represents a significant advance, as it automates the identification of regions by removing the expert-in-the-loop.

Although to the best of our knowledge the use of machine learning for local-nonlocal coupling has never been pursued before, machine learning for nonlocal simulations and/or computational mechanics is not a new concept. In what follows, we briefly review the most relevant research. We first focus on machine learning for peridynamics. Haghighat et al. <cit.> presented a nonlocal physics-informed framework with the peridyanmic differential operator <cit.> using deep learning. Kim et al. <cit.> presented Peri-net to analyze crack patterns using Deep Neural Networks. Nguyen et al. <cit.> presented peridynamic-based machine learning approach for up to two-dimensional structures. Xu et al. <cit.>  presented a machine-learning-based framework for peridynamic material models with physical constraints. Concerning nonlocal models, the following literature is available. Tao et al. <cit.> investigated nonlocal neural networks, nonlocal diffusion, and nonlocal modeling. Feng et al. <cit.> analyzed stochastic nonlocal damage using machine learning. Zhou et al. <cit.> learned nonlocal constitutive models using neural networks. You et al. <cit.> used data-driven learning of nonlocal physics from high-fidelity synthetic data. In <cit.>, de Moraes et al. used machine learning to predict the evolution of nonlocal micro-structural defects in crystalline materials.
Lal et al. <cit.> predicted nonlocal elasticity parameters using machine learning and molecular dynamics simulations. 

Concerning local models, we restrict ourselves to classical continuum mechanics. For a review of machine learning and data-miming approaches for continuum material mechanics, we refer the reader to <cit.>; whereas for computational solid mechanics, we refer to <cit.>. A machine-learning approach for physics-based material models for multiscale solid mechanics was presented in <cit.>. 

Finally, the use of machine learning for coupling methods (not involving nonlocal or peridynamics models) is also not new; in fact, Raymond <cit.> used machine learning to model coupled solid and fluid mechanics with mesh-free methods.

The paper is organized as follows. Section <ref> briefly introduces the models, coupling approaches, and discretization. In Section <ref> the methodology and the ML model which is the convolution neural network (CNN) are described. Section <ref> covers the data generation. In Section <ref>, we present numerical results for two cases using multiple node classification in Subsection <ref> and the main case of this work using node-wise classification in Subsection <ref>. Finally, Section <ref> concludes the paper.




§ DESCRIPTION OF THE MODELS AND COUPLING APPROACH



The model problem that we consider in this work deals with the quasi-static simulation of a one-dimensional bar with mixed boundary conditions at the two extremities. The deformations within the bar are assumed to be small enough to be described by linear elasticity. This section briefly describes the local and nonlocal models: linear elasticity and  linearized bond-based peridynamics <cit.>. We present the continuous formulations first and then provide their discretizations using finite differences for the local model and a collocation approach for peridynamics. 





 §.§ Classical linear elasticity model


We assume a one-dimensional bar of length ℓ occupying the open domain Ω=(0,ℓ) ⊂ℝ. Denoting the closure of Ω as Ω̅=[0,ℓ], we seek the displacement = (x), ∀ x ∈Ω̅, such that:

    - EA ”(x) = f_b(x),      ∀ x ∈Ω, 
    (x) = 0,      at x=0,
    
    EA'(x) = g,      at x=ℓ,

where E is the modulus of elasticity, A denotes the cross-sectional area of the bar, f_b=f_b(x) is the external body force density (per unit length), and g∈ℝ is a traction force applied at x=ℓ. We assume that f_b is chosen smoothly enough so that the solution  can be differentiated as many times as needed. For the sake of simplicity in the presentation, but without loss of generality, we assume a constant Young's modulus E and a constant cross-sectional area A for the bar. The model is presented here with a homogeneous Dirichlet boundary condition at x=0 and a Neumann condition at x=ℓ, but we shall also consider examples with homogeneous Dirichlet conditions at both ends in the numerical experiments.





 §.§ Linearized bond-based peridynamic model


The static peridynamic equation <cit.> reads in one dimension as

    - ∫_H_δ(x)κ(y) - (x)/|y-x| dy = f_b(x),

where the so-called horizon δ>0 determines the neighborhood of point x

    H_δ(x) = { y ∈ℝ; | y - x | < δ}.

In addition, κ=κ(x) is the peridynamic stiffness parameter evaluated at point x and u is the displacement for the deformed configuration. Following <cit.>, the peridynamic stiffness parameter κ is chosen such that the solution to the linearized peridynamic bond-based model is compatible with that to the continuum local model, i.e.,

    κδ^2/2 = EA,

or, equivalently,

    κ = 2EA/δ^2.

As far as the boundary conditions are concerned, we consider the same local boundary conditions (<ref>) and (<ref>), which are imposed using the variable horizon method described below. The solution to this problem will serve as the reference solution. For more details about boundary conditions for the peridynamic model, we refer the reader to, e.g., <cit.>.




 §.§ Variable horizon coupling approach




We consider in this work the Variable Horizon Coupling Method (VHCM) as introduced in <cit.>. A similar coupling approach was proposed in <cit.>. For a broader review of coupling methods, including optimization-based and energy-based coupling methods, we refer to <cit.>. As shown in Figure <ref>, we define the region in which we shall use the nonlocal model (<ref>) as Ω_δ = (a,b) and the region for the local model (<ref>) as Ω_e = Ω\Ω̅_δ. We suppose here that Ω_e is decomposed into the two subdomains Ω_1 = (0,a) and Ω_2=(b,ℓ) so that the end points of the bar are both adjacent to a domain governed by the local model. The main idea of the method is to make the horizon δ decrease to zero as one approaches the interfaces x=a and x=b so as to avoid the need to introduce an overlapping region around the interfaces. In this section, we describe the main ingredients of the VHCM method.

In the VHCM, the horizon δ_v(x) depends on the position x ∈Ω; it is usually constant in the nonlocal region far from the local-nonlocal interface and shrinks to zero when approaching the interfaces x=a and x=b. Several choices for the transition of δ_v(x) are possible. In one dimension, it was suggested in <cit.> to consider a piece-wise linear function, that is, given δ∈ℝ^+ constant,

    δ_v(x) = {[             x-a,     a < x ≤ a+δ,;               δ, a+δ < x ≤ b - δ,;           b - x,     b-δ < x < b, ].

as shown in Figure <ref>. It follows that Equation (<ref>) for the stiffness constant κ(x) now becomes:

    κ̅(x) = 2EA/δ_v(x)^2.




Using the variable horizon method, the coupling problem consists then in finding the displacements ∈Ω̅_e and u ∈Ω̅_δ such that

    - EA ”(x)    = f_b(x), 
      ∀ x ∈Ω_e, 
    
    - ∫_x-δ_v(x)^x+δ_v(x)κ̅(x) (y) - (x)/|y-x| dy    = f_b(x), 
      ∀ x ∈, 
    (x)    = 0, 
      at x = 0, 
    
    EA '(x)    = g, 
      at x = ℓ, 
    (x) - (x)    = 0, 
      at x = a, 
    (x) - (x)    = 0, 
      at x = b, 
    σ^+(u)(x) - EA'(x)    = 0, 
      at x=a, 
    σ^-(u)(x) - EA'(x)    = 0, 
      at x=b.

In other words, the displacements and stresses from the two models are matched at the interfaces, i.e. (x) = u(x) and σ^±(x) = E'(x), where the peridynamic stresses σ^±(x) <cit.> are defined as

    σ^±(u)(x)
    = δ/2∫_x-δ^x∫_x^z±δκu(y) - u(z)/|y-z| dydz - κδ^4/48 u”'(x).






 §.§ Discretization



Several approaches for discretizing the coupled problem (<ref>) are available <cit.>. Here we choose to discretize the classical linear elasticity model by a second-order finite differences scheme and the peridynamic model by a collocation approach, as introduced by Silling and Askari <cit.>. 

We now describe the discrete model associated with (<ref>) on the configuration shown in Figure <ref>. As customarily done in the literature <cit.>, we fix the value of the horizon δ and introduce a uniform grid spacing h chosen such that δ be a multiple of h, i.e. δ/h = m, with m a positive integer. Moreover, we choose the grid size h to be the same in each of the subregions of the computational domain. In other words, the numbers of intervals n_δ, n_1, and n_2 in Ω_δ, Ω_1, and Ω_2, respectively, are taken such that

    h = b-a/n_δ = a/n_1 = ℓ-b/n_2.

Since we assume that the grid points from the two models coincide at the interfaces x=a and x=b and that the grid is uniform in the three regions Ω_δ, Ω_1, and Ω_2, the grid points are simply given by:

    x_k = k h,   k=0,1,…,n,

where n = n_1 + n_δ + n_2. The numbering of the nodes is shown in Figure <ref>. In particular, we note that x_n_1 = a and x_n_1+n_δ = b. We consider the same numbering for the degrees of freedom associated with the displacements, that is, the nodal values will be denoted by:

    u_k ≈{   (x_k),          k = 0,…, n_1-1
       (x_k)= u(x_k),          k = n_1
       u(x_k),          k = n_1+1,…, n_1+n_δ-1 
       (x_k)= u(x_k),          k = n_1+n_δ
       (x_k),          k=n_1+n_δ+1,…,n
    .

where we have used the continuity of the displacement at the two interfaces, that is, (x) = u(x) at x=a and x=b. The total number of degrees of freedom is then equal to n+1, being each degree of freedom associated with one grid point.




The discretization of the VHC problem (<ref>) leads to the following system of equations, in which we assume m=2:


  * Dirichlet BC at x=0:

    u_0 = 0.



  * In Ω_1: 

    - EA u_k-1-2u_k+u_k+1/h^2 = f_b(x_k),    k=1,…,n_1-1.



  * At interface x=a:

    EA 
    - 2u_k-3 
    + 9u_k-2 
    -18u_k-1 
    +11u_k/6h 
    - σ_h^+(x_k) = 0,    k=n_1.



  * In Ω_δ: 

    - κδ^2/2u_k-1-2u_k+u_k+1/h_δ^2 = f_b(x_k),
             k=n_1+1,
       - κδ^2/2u_k-2+4u_k-1-10u_k+4u_k+1+u_k+2/8h_δ^2 = f_b(x_k),
             k=n_1+2,…,n_1+n_δ-2,
       - κδ^2/2u_k-1-2u_k+u_k+1/h^2 = f_b(x_k),
             k=n_1+n_δ-1.



  * At interface x=b: 

    EA 
    -11u_k
    +18u_k+1
    - 9u_k+2
    + 2u_k+3/6h 
    - σ_h^-(x_k) = 0,    k=n_1+n_δ.



  * In Ω_2: 

    - EA u_k-1-2u_k+u_k+1/h^2 = f_b(x_k),    k=n_1+n_δ +1,…,n-1.



  * Neumann BC at x=ℓ:

    EA 
    - 2u_k-3
    + 9u_k-2
    -18u_k-1
    +11u_k/6h = g,    k=n.



We recall that if E and A are constant, then the quantity κ̅(x)δ_v(x)^2 = κδ^2 remains also constant. Moreover, the stresses σ_h^± in Equations (<ref>) and (<ref>) are obtained as approximations of σ_± (<ref>). More specifically, assuming that the solution u(x) is sufficiently regular in Ω_δ, the stresses are approximated using one-sided third-order finite differences stencils of the first derivative as

    σ_h^+(x_k)
    = κδ^2/2
    -11u_k
    +18u_k+1
    - 9u_k+2
    + 2u_k+3/6h,    k=n_1,
       σ_h^-(x_k)
    = κδ^2/2
    - 2u_k-3
    + 9u_k-2
    -18u_k-1
    +11u_k/6h,    k=n_1+n_δ.

This choice ensures that the discretization method has a degree of precision of three, in other words, that any polynomial solution of degree up to three will be computed exactly.




§ METHODOLOGY



We propose to employ deep learning techniques to identify the regions in which one should use the local model or the nonlocal model, by solely using load information (i.e. external forces), as illustrated in Figure <ref>. More specifically, we shall use CNNs as our deep learning tool. We describe below CNNs and the specific architecture that we considered to classify the nature of the grid points.




 §.§ Overview


Our main objective is to develop an automated approach to identify the nonlocal region Ω_δ and the local regions Ω_1 and Ω_2, shown in Figure <ref>, to be used in a coupling method. From now on, the nodal values of the solution located in the regions Ω_1 and Ω_2 will be indexed by LM (Local Model), whereas the nodal values located in Ω_δ will be indexed by NLM (Nonlocal Model). Although this is a simplified setting that does not fully reflect the numerical tests presented below, we choose to describe the method using one single nonlocal region for the sake of clarity. The extension to more than one nonlocal domains is straightforward.




An overview of the proposed automated identification process is shown in Figure <ref>. The input data consists of the external load vector f_b(x_k), k=1,…,n-1, evaluated at all interior grid points x_k, see Figure <ref>a. The output of the process is a label associated with each node that takes the value 0, if the node is classified as belonging to the local region, or 1 otherwise (i.e. the node belongs to the nonlocal region). It is essential to note that the nature of our approach is supervised learning, that is, each load input during training corresponds to a known classification. As this requires training data, the load input (see Figure <ref>a) is also used to generate the ground-truth data that will serve as the reference configuration during training. The generation, as shown in Figure <ref>b, is conducted using the following procedure: a coupling configuration (local-nonlocal domain assignment) is added to the training set when it provides a coupled solution whose error, with respect to the fully nonlocal solution, is below a given tolerance. More specifically, we first compute the reference displacement u_NLM using the nonlocal model in which all nodal points are NLM nodes ∘. Then, given a coupling configuration (in which the nodes are either LM nodes ∙ or NLM nodes ∘), we compute the displacement u_VHCM using the method described in the previous section, see Figure <ref>b.

As mentioned above, the reference configuration is included in the training dataset if the relative error in the ℓ_2 norm is below a given tolerance ε, namely


    ℰ(u_NLM,u_VHCM) =‖ u_NLM - u_VHCM‖_ℓ_2/‖ u_NLM‖_ℓ_2 < ε.




For more details about the data generation, the reader is referred to Section <ref>. The loads and the reference configurations of the local-nonlocal splitting provide the training data for the CNN in a supervised manner. Once trained, the CNN model takes as input an (unseen) external load f_b(x), see Figure <ref>c, and outputs the predicted configuration (i.e. it classifies the nodes as either LM ∙ or NLM ∘ at each grid point) in the domain, see Figure <ref>d. More details of the CNN architecture are given in the following section.






 §.§ Convolution neural networks



In recent years, deep learning has yielded impressive results across a spectrum of challenges, such as in visual recognition, speech understanding, natural language processing, and even inference of physical systems responses. Among different deep neural network architectures, CNNs have been investigated and analyzed in depth <cit.>. CNNs are particularly effective for tasks involving data with a grid-like topology, such as images, audio spectrograms, and sequential data arising for instance in natural language processing,  with applications ranging from image classification and object tracking to speech and natural language processing <cit.>.

The CNN structure consists of various fundamental components, including convolutional, pooling, and fully connected layers. The standard structure involves a series of convolution layers followed by a pooling layer. The layers are then succeeded by one or more fully connected layers. The phase in which input data undergoes transformation into output through the layers is the so-called forward propagation process <cit.>. In the next two sections, we give a brief description of CNNs' features for the non-expert reader.




  * Convolutional Layers: 
Convolutional layers are the cornerstone of CNNs. These layers use learnable kernels that convolve over the input data, enabling the network to detect local patterns and features. More precisely, a kernel is typically a small array, usually of size 3 × 3, which is systematically passed through an input array, also referred to as a tensor. At each position within the tensor, an element-wise multiplication is performed between the elements of the kernel and the corresponding elements in the input tensor. A bias is then added to the sum of these products in order to yield the output value, which occupies the corresponding position in the output tensor—termed “feature map”. This process is reiterated by employing multiple kernels to generate diverse feature maps. These feature maps encapsulate various attributes of the input tensors, such as distinct kernel functions serving as distinct feature extractors. The operation of a convolutional layer can be expressed as:

    A_j = σ_ReLU(∑_i=1^N C_i K_i,j + B_j ),

where C_i denotes an input matrix, K_i,j a kernel matrix,  N is the number of feature maps from the previous layer, and B_j a bias value. Finally, the rectified linear unit activation function σ_ReLU <cit.>,

    σ_ReLU(x)=max(0,x),

is applied element-wise to the matrix, resulting in the output matrix A_j.
The primary task of the learning process involves identifying sets of suitable kernel matrices capable of extracting distinctive and discriminative features for the output tasks. To achieve this, the back-propagation algorithm, commonly employed to optimize the connection weights in neural networks, can be adapted. In this context, it serves to train both the kernel matrices and biases, functioning as shared neuron connection weights.


  * Pooling Layers:
Pooling layers downsample the dimensions of feature maps while retaining the salient information; thus, their main effect is an intelligent feature-dimension reduction. To effectively reduce the number of output neurons, pooling algorithms usually combine neighboring elements in the convolution output matrices. 
Notable among these algorithms are max-pooling and average-pooling <cit.>.  In this study, we implement a max-pooling layer with a  2 × 1 size, which selects the highest value from the two neighboring elements of the input feature map, thereby generating a single element for the output feature map. 




  * Fully Connected Layers:
Fully connected layers process the flattened feature vectors to make the final predictions. A fully connected layer operation with input z and output y_out is represented as

    y_out = σ(Gz + v)

where y_out is the output, z the input vector, G the weight matrix, v the bias vector, and σ a given activation function <cit.>.








  
The selected CNN architecture: We propose to use the following CNN architecture for the identification of local and nonlocal regions. First, as a preprocessing step, the input of the network is normalized to a load vector with unit variance and zero mean. The architecture is composed of convolutional layers with ReLU activation functions for feature extraction, followed by pooling layers for spatial reduction, and fully connected layers for the final output. This is shown in Fig. <ref>. The first convolution layer uses 32 filters with 3 × 1 kernels. The second convolution layer has 64 filters with 3 kernels, and the third convolution layer has 128 filters with 3× 1 kernels. Each convolution layer is followed by a max pooling layer with pool size of two. The output tensor from these three conv-pooling blocks has a dimension of 30 × 128. The following two layers consist of one flattened layer, followed by a fully connected layer with 64 neurons. The last layer of our CNN model is a fully connected (dense) layer, where the number of neurons matches the output size. In our study, we explore two cases with distinct output sizes. In the first case, the CNN receives as input the full vector of the second derivative of the load functions at the grid points. The resulting output is a vector of labels, having the same size of the input load vector. Herein, the output size is 257, leading to a fully connected layer with 257 neurons. 
Meanwhile, in the second case, the CNN receives a window as input. The resulting output is a label for the central node. Herein, the output size is then one, corresponding to a single neuron in the final layer. Further details on these cases are discussed in Section <ref>. The sigmoid activation function is used in the last layer across both cases. It is important to note that our experimental investigations explore a variety of network architectures, using different numbers of layers and neurons. Out of all the architectures evaluated, the one that produces the best results is the one that is presented herein.



  
Computational resources and model parameters:
The CNN model is trained by using the Adam optimizer <cit.> and a batch size of 32. The training proceeds for a maximum of 200 epochs; however, an early stopping callback is implemented to avoid overfitting. When early stopping is used, training will end as soon as the validation dataset's loss does not improve <cit.>. The learning rate is set to the commonly used value of 0.001. The proposed model is implemented in TensorFlow 2.10, with
Keras in Python 3.10.11. 
































 §.§ Classification evaluation metrics


















In the context of binary classification problems such as ours, the classification capability of the trained network is usually evaluated with a confusion matrix.

The confusion matrix displays the correctly or incorrectly predicted results. It is based on four key elements: true positives (TP, correctly predicted LM), false positives (FP, wrongly predicted LM), true negatives (TN, accurately predicted NLM), and false negatives (FN, wrongly predicted NLM). To evaluate the confusion matrix, we use two induced metrics, namely the accuracy metric and the F1-score, which are the most used metrics in classification <cit.>:



  * The accuracy metric measures the
ratio of correct predictions over the total
number of instances evaluated. Higher values for accuracy are indicative of better performance. Formally,

    Accuracy = TP + TN/TP+FP+TN+FN.















  * The F1-score is the harmonic mean of “precision” and “recall”. It ranges from 0 to 1, where a higher value indicates better classification performance. Formally,

    F1-score = 2 ×Precision×Recall/Precision + Recall,

where 

    Precision = TP/TP + FP, 
       Recall = TP/TP + FN.



We note that in case 1, see Section <ref>, we report the average over all samples for all metrics (i.e. confusion matrix, accuracy, and F1-score).




§ DATA GENERATION



The CNN used for classification is trained on input-output pairs of external loads and corresponding coupling configurations (inducing a coupled solution with a certain accuracy); this pipeline is illustrated in Figure <ref>. In this section, we introduce the definition of the load functions and present the steps involved in the generation of the coupling regions.





 §.§ Families of loading functions



In this section, we consider several classes of loading functions, including singular functions that induce discontinuous solutions with a finite jump to mimic fracture phenomena in higher dimensions. One advantage of peridynamics is that the model does not feature any spatial derivatives thus allowing for solutions with jumps. 



  * 






            


The first family of load functions, f_1, is based on functions <cit.> featuring a singularity at x = 0.5 and defined as:

    f_1(x)=
    
    0,     x∈ [0, 0.5  - δ),
    1/2δ^2 -δ +3/8 + (2δ - 3/2  -lnδ)x 
         + (3/2 +lnδ) x^2 -(x^2 -x)ln((1/2-x)),     x∈ [0.5-δ, 0.5),
    1/2δ^2 -δ +3/8 + (2δ + 3/2  +lnδ)x 
         - (3/2 +lnδ) x^2 +(x^2-x) ln(x-1/2),     x∈ [0.5, 0.5+δ),
    
    1,      x∈ [0.5+δ, 1.0].



The function f_1 and the corresponding fully nonlocal solution u_NLM are shown in Figure <ref>(a-b) for δ = 0.03 and h=δ/8. 

In cases where the singularity is located at a point x_jump different from 0.5, we will consider the loading function f_1(x+(0.5 -x_jump)).








            




  * The second family of load functions, f_2, featuring a singularity at x = 0.5 <cit.>, is defined by:

    f_2(x)=
    0, 
        x∈ [0, 0.5 - δ),
    1/2δ^2(lnδ-ln(1/2-x) ), 
        x∈ [0.5-δ, 0.5),
    1/2δ^2(ln(x-1/2)-lnδ )),
        x∈ [0.5, 0.5+δ),
    
    0,  
        x∈ [0.5+δ, 1.0].




The function f_2 and the corresponding full nonlocal solution u_NLM are plotted in Figure <ref>(c-d) for δ = 0.03 and h=δ/8. 

In cases where the singularity is located at a point x_jump different from 0.5, we will evaluate the loading function as f_2(x+(0.5 -x_jump)).







       






















 





























       







    





  * The study also includes the family of loads characterized by polynomial solutions of degree three and lower. The loading functions f_3(x) are in this case linear of the form:

    f_3(x)=c_1x + c_2,

an example of which is shown in Fig. <ref>(e-f) with c_1=18 and c_2=4.
 





  * We consider as well the following family of load functions:

    f_4(x) = tanh1/t(x - 1/2),

where the parameter t controls the variation in the solution around the point x=0.5. An example of function f_4 is shown in Fig. <ref>(a-b) for t=0.0005.

  * Finally, we also include a family of smooth load functions that induce continuous solutions with local behavior: 

    f_5(x) = e^-400(x-c)^2.

Such a function is a smooth representation of a point-wise load applied at the point x=c. The load and corresponding solution is shown in Fig. <ref>(c-d) for c=0.6.








 §.§ Generation of the coupling regions



We present the steps used to generate the output data for the CNN training, i.e. we describe how the reference configurations of the coupling regions are constructed based on the loading functions.

For the sake of simplicity, all experiments focus on a one-dimensional bar with a material property of E=1 and cross-sectional area A=1, resulting in κ = 2/δ^2. We consider the configuration of Figure <ref> with ℓ= 1 and interface locations a and b. We choose the values of a and b such that grid points coincide the two interfaces. The domain is partitioned into sub-intervals of size h = 1/n, with n given. The horizon is taken as a multiple of h, i.e. δ=mh, with m=8. 

We acknowledge that implementing non-uniform grids in our approach poses significant challenges, particularly in integrating them with the standard CNN architecture and our coupling method. Therefore, and for the purposes of this proof-of-concept, we will employ a fixed discretization grid to evaluate the innovative aspects of our proposed approach. 

To avoid any issues with applying the local boundary conditions, we assume that Ω_δ = (a, b) is such that a > 2δ and b < ℓ-2δ. We will nevertheless vary the positions a and b of the interfaces in order to consider solutions with a discontinuity located in the interval (2δ,ℓ-2δ). 



We now describe the process of generating the reference coupling regions used in the training for the special cases f_1 and f_2.
Since these loads are smooth away from their respective singularities, it seems sensible that the region where one would use the nonlocal model in the local-nonlocal coupling approach should contain the singular point.
We thus compute the fully nonlocal reference solution u_NLM and choose the smallest interval around the singularity, i.e. Ω_δ=(x_jump-α,x_jump+α) where α>0 is a multiple of the mesh size h, such that the relative error ℰ(u_NLM,u_VHCM) (<ref>) between the fully nonlocal and coupled solutions is below a given tolerance ϵ. In our experiments, we fix ϵ to 0.01.
It is clear that the proposed process is directly influenced by the specific features of the chosen loading functions. 

In Figure <ref>, we show an example of an optimal nonlocal region configuration using f_1 and f_2 with singularity at x = 0.5. The displacements u_NLM and u_VHCM using f_1 are shown in Figure <ref>(a). The value of ℰ(u_NLM,u_VHCM) using f_1 is 7.026 × 10^-3.  Figure <ref>(b) presents the displacement u_NLM and u_VHCM using the loading type f_2. The value of ℰ(u_NLM,u_VHCM) in this case is 7.027 × 10^-3. In both cases, the error in the coupled solution is below the selected threshold. We also show the pointwise error between the nonlocal solution (NLM) and the coupled solution (VHCM) in Figure <ref>(c-d), to illustrate that the local error remains relatively small at all points.




As mentioned in Section <ref>, the solutions to the local and the nonlocal models coincide for polynomials of degree up to three; furthermore, it is known that the fourth derivative of the solution is an indicator of the extent of the nonlocality of the solution <cit.>. As a pre-processing step, we leverage this information by taking the second derivative of the loading functions rather than the original functions themselves. Here, the second derivatives are simply evaluated using finite differences approximations of the loading functions. This step will help in accurately capturing the underlying behavior for the load and serve as a nonlocality indicator. The second derivatives of these loads are then used as inputs for our ML model.




We now proceed to explain the process when using load f_3. The coefficients c_1 and c_2 in Eq. (<ref>) are random values within the range [-30,30] and [-10,10], respectively.

In this case, the reference configuration is  characterized by a fully local region as the local model and the nonlocal model lead to the same solution. This  choice is made to anchor the training of our model on an instance where the behavior is fully local. This is pivotal for enriching the learning process of the model, enabling it to discern and accurately predict the nuances of fully local configurations alongside the more complex local-nonlocal configurations. Notably, the second derivative of f_3 is always zero and is used as input for our ML model.





 §.§ Approaches for input-output selection

We recall that we consider two strategies. In the first one, we utilize the complete dataset over the full computational domain as input. In this scenario, the nature of the problem is akin to the simultaneous classification of all nodes. In the second approach, we transform the datasets into windowed inputs, and solve a point-wise classification problem. 



  §.§.§ Full-domain input data












The first strategy consists of considering the second derivative of the load functions at all grid points simultaneously following the generation of data as explained in Section <ref>. In this context, we increase the model's data volume and enhance its ability to generalize by introducing specific transformations to the functions f_1 (<ref>) and f_2 (<ref>). Specifically, we incorporate the negation operation, which complements our existing data by providing contrasting cases.

The second derivatives of f_1, f_2, and f_3 are estimated and used here as input vectors.
In this approach, the input-output pairs are given by a full second derivative of the load vector-full labels vector. The CNN receives the full second derivative  vector as input. The resulting output is a vector of labels, having the same size of the input vector. Each output node has either the LM or NLM label, according to the region (local or nonlocal) they belong to. 
After prediction, we implement a post-processing step where nonlocal regions containing fewer than eight nodes are treated as local (LM) since we are using the variable horizon coupling method (VHCM) with m=8.




























                       




























                       




   
















  §.§.§ Window-based input data



The second (and more effective) strategy proposed in this paper is the windowing approach. This approach proves to be a highly effective strategy for handling data with several singularities. 

In this approach, each node, denoted as x_i, is treated individually.
In the windowing approach, we reformulate the classification problem as point-wise classification; an illustration of the window selection is presented in Figure <ref>. By utilizing window configurations around each data point, we are effectively segmenting the load's second derivative, enabling a focused and precise classification strategy at each point. This method allows for a more focused classification approach, providing a detailed insight into the load characteristics within specific windowed segments. The induced input-output pairs are then given by window-label. The CNN receives a window of load's second derivative as input. The resulting output is a label for the central node.  The label is either the LM or NLM, according to the region (local or nonlocal) the central node belong to. In a post-processing step, we apply the same process detailed in Section <ref>.
This approach aims at enhancing the model's capability to detect patterns and variations in the data, making our analysis more precise and adaptable to different load scenarios.

We present in Figure <ref> an illustration of the window selection for the window-based approach. Considering two points and their corresponding windows, the window spans two horizons on the left and on the right of the central point, thereby covering a total of four horizons. Figure <ref> shows two central nodes x_i and x_j, each with its associated windows W(x_i) and W(x_j). These windows are used as input data for the CNN enabling, classification labels based on their region. The output labels will be then LM for x_i since it is located in the local region (represented by ∙) and NLM for x_j since it is located in the nonlocal region (represented by ∘). 




We note that this approach can be also applied and used in higher dimensions. For instance, in two or three dimensions, the windows might transition from one-dimensional segments to two-dimensional areas or even three-dimensional volumes. By varying these windows, we aim to capture and understand how nonlocal effects act within specific regions or volumes across the broader dimensional space.




 §.§ Case studies


In this section, we describe the different scenarios and case studies considered in our tests. A summary is reported in Table <ref>.




  §.§.§ Case 1: Full-domain input data



In this case, we generate the data as explained in Section <ref>. These simulations result in a total of 3,084 data samples, as shown in Table <ref>(A). The data is first split into training, testing, and validation sets. The training set is used to train the CNN model, where the model learns the relationship between inputs and outputs, adjusting its parameters (like weights and biases) accordingly. The validation set is used during the training process; it helps in evaluating the performance of the model on a dataset it has not be subjected to before. The testing set is used to check how well the model can generalize and helps prevent overfitting. In other words, the testing set allows one to evaluate the performance of the model on new, unseen data <cit.>.  We split the dataset by assigning 75% of the samples to the training set, 10% to the validation set and 15% to the testing
set. The size of each data sample is 257, resulting in a training set of size 589× 257 and a validation set of size 78× 257, while the testing set has size 119× 257. The generation and partitioning of the dataset into training, validation, and testing sets necessitates a total of 138.77 seconds (Table <ref>(B)).

In the full-domain input data approach, the CNN receives as input the full vector of the second derivatives of the load functions at the grid points. The resulting output is a vector of labels, having the same size of the input load vector. Each output node has either the LM or NLM label, according to the region (local or nonlocal) they belong to.
The CNN model is structured as multiple node classification model; the results are detailed in Section <ref>. A summary of this test case is shown in Table <ref>.







 








  §.§.§ Case 2: Window-based input data



In this case, we use the window-based approach, in which the dataset is generated and preprocessed using the methodology explained in Section <ref>. The initial dataset, consisting of 786 complete input-output vector pairs, is partitioned into distinct sets for training, validation, and testing. The training set constitutes 75% of the total data, the validation set allocates 10%, and the remaining 15% is designated for the testing set. We then implement the windowing approach, segmenting each vector into windows of data. Subsequently, we perform a deduplication process, removing duplicated input windows-output label  in each set and between sets (train-test, train-validation, and validation-test), ensuring unique samples across all sets.
The resultant windowed dataset comprises a total of 30776 data samples (Table <ref>(A)). The training set, the validation set, and the testing set contain 22802, 3476, and 4498  data samples, respectively. We note that the size of input samples depends on the window size.
The generation and data splitting into windows require a total of 182.37 seconds (Table <ref>(B)). We point out that we introduced the deduplication step to ensure fairness in the testing process. While straightforward in  1D experiments with uniform grids, in practice and especially in higher dimensions, deduplication is nontrivial. However, it is often the case that solutions have common features in certain regions of a computational domain so that having duplicates is very natural and actually helps the learning process. Thus, our results show that even when we stress test the model by removing favorable data, the model still performs well (Section <ref>).

The CNN inputs are the windows and the corresponding output is the label (local or nonlocal) of the central node. Thus, the model has been restructured as a node-wise classification model. We recall that the labels are either LM or NLM depending on the classification (local or nonlocal). An overview of this case is summarized in Table <ref>.




§ NUMERICAL RESULTS






 §.§ Case 1: Full-domain input data
 












In the process of training our CNN model, both training and validation loss were computed to monitor the model's learning progress and its ability to generalize. The training proceed for a maximum of 200 epochs; however, an early stopping callback is implemented to avoid overfitting. The training time is notably brief, amounting to only 30.19 seconds. However, testing the model requires a duration of 0.6 seconds (Table <ref>(B)). Figure <ref>(a) presents the training and validation losses over the training epochs.
The performance of the model is assessed using the evaluation metrics in Section <ref> and is shown in Table <ref>. The model's performance evaluation produces an average accuracy of 0.99 and an F1-score value of 0.94 (Table <ref>).
We present the corresponding average confusion matrix in Figure <ref>(b). Within the confusion matrix, the model predicts correctly 99.73% of the LM instances and 96.34% of the NLM ones. These percentage-based values illustrate the model's high precision and recall scores, and show its robustness in accurately classifying nodes within local and nonlocal regions for the training dataset.




To evaluate the results after regions detection, we estimate the error in Eq (<ref>) between the fully nonlocal displacement u_NLM and the coupled displacement u_VHCM, whose nonlocal region is the one predicted by the CNN. An example of the results for both f_1 and f_2 is reported in  Figure <ref>, where, in both cases, the discontinuity is at x=0.71. The value of ℰ(u_NLM,u_VHCM) using f_1 is 5.5282 × 10^-5; while using f_2 is 5.5285 × 10^-5.

For the test set, the predicted nonlocal region consistently induces an error in the coupled solution much lower than the tolerance used during training . Across all cases, the estimated errors fall within the range of 2.18 × 10^-6 to 5.32 × 10^-3 using f_1 and 2.17 × 10^-6 to 1.36 × 10^-4 using f_2 (Table <ref>).









 §.§ Case 2: Window-based input data




The results in the previous Section <ref> showed promise that machine learning could support the selection of local-nonlocal coupling regions; however, the full-load strategy does not have satisfactory generalization properties when testing using unseen data. In this section, we illustrate how the windowing approach circumvents these limitations.

As before, the number of epochs is set to of 200; however, an early stopping callback is used to avoid overfitting.  Training the model with windowed data take 137.17 seconds while testing requires a duration of 10.9 seconds (Table <ref>(B)). Figure <ref>(a) presents the training and validation losses over the training epochs. While both losses decrease overall, the close tracking between training and validation loss suggests that the model is not overfitting to the training data. The plot shows that the callback effectively prevented overfitting by stopping training at epoch 50 when the validation loss showed no further improvement.
The performance of the model was assessed using the evaluation metrics are shown in Table <ref>. The performance metrics show an accuracy of 0.96 and an F1-score of 0.97;



the corresponding confusion matrix, reported in Figure <ref>(c),  shows that the model predicted correctly 94.03% of the LM instances and 97.18% of of the NLM ones. Also in this case, the model has high precision and recall scores.

For the test set, the predicted nonlocal region induces an error in the coupled solution much lower than the tolerance ϵ used during training all cases using both f_1 and f_2. Across all cases, the estimated errors fall within the range of 3.50 × 10^-6 to 2.34 × 10^-3 using f_1 and 3.49 × 10^-6 to 1.97 × 10^-5 using f_2 (Table <ref>). The lowest errors are obtained when using the original version of loads f_1 and f_2, which means without transformations. 





To evaluate the windowing approach's ability to identify multiple discontinuities, we test our pre-trained model, which has not been trained on loads with more than one discontinuity, using a load with two discontinuities. The results are presented in Figure <ref>. The predicted configurations of local and nonlocal regions are shown in Figure <ref>(a), while the displacements u_NLM (represented by ▪) and u_VHCM (represented by ∙) are shown in Figure <ref>(b). The corresponding error, denoted as ℰ(u_NLM, u_VHCM), is 3.299 × 10^-3.



In addition, we test our model using load f_4 (Eq. <ref>) and f_5 (Eq. <ref>); these loading functions do not belong to the training set. Thus, with these examples, we test mild generalization properties of the algorithm. 





For f_4, we consider two scenarios; one with t=0.05 and another with t=0.0005. In all cases, the second derivative is estimated and subsequently segmented into windows, serving as the input for the pre-trained model. Once again, to evaluate the precision of the region detection, we estimate the error between the fully nonlocal displacement (u_NLM) and the coupled displacement (u_VHCM) using Eq. (<ref>). Figure <ref>(Left) shows the predicted configurations of local and nonlocal regions. Figure <ref>(Right) shows u_NLM, represented by ▪ and u_VHCM, represented by ∙.

For f_4 with t=0.05, the prediction labels and configurations of local and nonlocal regions are showcased in Figure <ref>(a-b). The corresponding error value, denoted as ℰ(u_NLM, u_VHCM), is 5.30 × 10^-5. When testing f_4 with t=0.0005, the findings are similarly reported in Figures <ref>(c-d), with an error value of 2.55 × 10^-5.

For f_5, as defined in Eq. <ref>, we follow the same process. The results are presented in Figure <ref>(e-f), and the corresponding error ℰ(u_NLM, u_VHCM) is 1.30 × 10^-4.












































§ CONCLUSION


This paper shows a proof of concept for a ML-based identification of local-nonlocal coupling regions that is independent of the local and nonlocal models utilized. Specifically, the load used as identification input is independent of the two models and the method used to generate local and nonlocal solutions is irrelevant. Furthermore, while we only considered a non-overlapping approach, a coupling approach with overlapping regions could also be used (examples can be found in  <cit.>). 

Among the two approaches we presented, the full-domain approach yields promising results for the classification of the local and nonlocal regions as long as the solutions do not feature multiple discontinuities. To overcome this limitation, we introduced the windowing approach based on a node-wise classification. The latter shows promising results even in the presence of several discontinuities without retraining (see Section <ref>). This model demonstrated robust performance, as reflected by high evaluation metrics, with an accuracy of 0.98 and an F1-score of 0.97. The results were further validated by the corresponding confusion matrix (see Figure <ref>(c)), which provides a detailed breakdown of the model's classification accuracy.
Nevertheless, a minor limitation in our model's predictive precision is noticed. Specifically, our analysis identified a misclassification rate of 5.21% of instances as LM (nodes in Local region) and 1.14% as NLM (nodes in nonlocal region). These relatively low values of the total predictions highlight areas of potential improvement.



The proof of concept approach introduced in this work opens several avenues of improvement and extension. In the results presented in this study, the horizon δ=m· h is constant; for more realistic tests, varying values should be investigated (the work in <cit.> can be used as a guide for the choice of h and δ). As for the boundary conditions, we use homogeneous Dirichlet boundary conditions, while in real-world applications mixed conditions are usually prescribed. More importantly, an extension to two- and three-dimensional problems is mandatory and it is the subject of our current work. Furthermore, stress testing the approach by using more complex forcing functions from a larger subspace should also be explored (e.g. randomly generated functions with prescribed regularity properties).

Another relevant aspect not directly addressed in this work is the presence of cracks and fractures. These features were mimicked in the current study by the presence of discontinuities. In two- and three-dimensional settings, we plan to incorporate these features into the machine learning model and validate our results against real experiments.

In terms of machine learning algorithms, the windowing approach showed promising results and is the candidate strategy to proceed with in higher dimensions. However, the open question is the choice of the shape of the window. Possible candidates in two dimensions could be rectangles or ellipses. More generally, the shape of the window should be chosen such that it matches the support of the nonlocal kernel function. Alternatively, one could parameterize a family of two-dimensional shapes, in which case the parameters would become learnable parameters in the training process. Another open question is the architecture of the CNN for higher dimensions as the input layer and output layer need to be adapted to the higher dimensions.

Additionally, due to the fact that CNNs fail in the presence of nonuniform unstructured grids, our future exploration includes the use of graph neural networks (GNNs) <cit.>. Given the graph-like structure of our problem, leveraging GNNs could offer a compelling approach to model and predict complex interactions within nodes in local and nonlocal regions. The adaptability of GNNs to capture relationships among various regions' nodes makes them a promising approach especially for higher dimensions.



§ FUNDING STATEMENT


Sandia National Laboratories is a multimission laboratory managed and operated by National Technology & Engineering Solutions of Sandia, LLC, a wholly owned subsidiary of Honeywell International Inc., for the U.S. Department of Energy’s National Nuclear Security Administration under contract DE-NA0003525.

This paper describes objective technical results and analysis. Any subjective views or opinions that might be expressed in the paper do not necessarily represent the views of the U.S. Department of Energy or the United States Government.

abbrv
  



























