The quality of mesh generation has long been considered a vital aspect in providing engineers with reliable simulation results throughout the history of the Finite Element Method (FEM). The element extraction method, which is currently the most robust method, is used in business soft- ware. However, in order to speed up extraction, the approach is done by ﬁnding the next element that optimizes a target function, which can result in local mesh of bad quality after many time steps. We provide TreeMesh , a method that uses this method in conjunction with reinforcement learning (also possible with supervised learning) and a novel Monte-Carlo tree search (MCTS) ([Coulom(2006) ], [Kocsis and Szepesvári(2006) ], [Browne et al.(2012) ]). The algorithm is based on a previously proposed approach ([ Pan et al.(2021) ]). After making many improvements on DRL (algorithm, state-action-reward setting) and adding a MCTS, it outperforms the former work on the same boundary. Furthermore, using tree search, our program reveals much preponderance on seed-density-changing boundaries, which is common on thin-ﬁlm materials. 1 Introduction 1.1 FEM mesh generation In Figure 1, there are two techniques to mesh a given boundary: structured ([ Thompson et al.(1998) ]) and unstructured meshes ([ Owen(1998) ]). All elements in a structured mesh are ﬁrst arranged in a regular way before being mapped to the real-world boundaries. As a result of the complex boundaries, the mesh quality may be poor. Unstructured mesh can be categorized into indirect and direct methods. Indirect methods require an intermediate triangular mesh, which is prone to instability and has a restricted number of vertices ([ Garimella et al.(2004) ]), whereas direct methods generate quadrilateral pieces directly, avoiding some potential problems ([ Zeng and Cheng(1993) ]). Our research is based on the element extraction method, which is the most popular of the direct methods ([ Docampo-Sanchez and Haimes(2019) ]). 1.2 Previous works MCTS with deep neural networks (NN) can defeat world champions in the game of Go, according to previous research in [Silver et al.(2016) ] and [ Silver et al.(2017) ]. NN may train with self-play from a random policy to a top player using Deep Reinforcement Learning (DRL). [ Yao et al.(2005) ] and [ Pan et al.(2021) ] also reported on the notion of utilizing the Artiﬁcial Neural Network (ANN) and DRL on FEM mesh generation. The Advantage Actor-Critic (A2C) approach, as shown in Figure 2, trains a policy and a value NN to provide solutions to diverse boundaries. After that, they perform behavior cloning by extracting experience (also known as data-set) and training another NN. With a self-evolving procedure, the supervised-learning-based NN FreeMesh-S can address new boundary challenges. Undergraduate Student, Theoretical and Applied Mechanics, Dynamics F1, 443 Huangshan Rd, Shushan District, Hefei, Anhui, China Generate Plane Quad Mesh with Neural Networks and Tree Search A P REPRINT de-mappingmapping mesh (a) Structured mesh removedistribu... mesh (b) Unstructured mesh Figure 1: Two ways to mesh a given boundary. Indirect unstructured meshes have a transformation process while direct ones generate quadrilateral elements directly. A2C meshing resultExperience extractionFNN training FNN modelFNN meshing result FNN meshingStep 1 Step 2 Figure 2: A previously proposed method to supervise an FNN with experience extracted from itself and another DRL program. 1.3 Our approach According to [ Hornik(1991) ], standard multi-layer feed-forward networks with activation function may arbitrarily well approximate any continuous function if there are enough hidden units available. Therefore, once enough experience has been extracted from Reinforcement Learning (RL), we can construct an MLP or ResNet ([ Shan et al.(2016) ]) architecture to solve new boundaries. To improve the quality of data-set in FreeMesh-S ,TreeMesh , our semi-open- source program makes many improvements on FreeMesh-S . We replace the A2C DRL algorithm with the more powerful DPPO algorithm and add a brute force MCTS to ﬁnd nearly-optimal meshes. We use a new state tensor, an action tensor with a smaller range, and a sophisticated reward function to speed up convergence. On the whole, TreeMesh aims to use DRL and MCTS to discover an optimal, or nearly-optimal, solution to some problematic boundaries in a large-scale boundary (e.g. sharp angles, bottleneck region, unevenly distributed segments, and holes) with a pre-trained model or from scratch. Further technical improvements in DRL and MCTS algorithms and solving procedures are described in Section 2, mesh results are exhibited in Section 3. We have released a single-threaded version on GitHub. 2 Methodology 2.1 Reinforcement learning in TreeMesh Our method continues to use the actor-critic framework, including a policy network (St;)and a value network V(St;w). Both NN use the same input tensor, which includes standardized node coordinates and mesh progression. The value network evaluates the mean return value, whereas the policy network outputs the coordinates for a new element. On one hand, we use the back propagation of the value loss function to assist it in providing increasingly precise evaluations after receiving trajectories. On the other hand, we calculate the time difference error and the probability of choosing each action to back propagate the policy loss function. We adopt linear layers and leaky ReLUs in the NN architecture. 2 Generate Plane Quad Mesh with Neural Networks and Tree Search A P REPRINT a. Sampling S0 S1 A0~π0 V0S2 A1~π1 V1···At~πt VtSt Rt ···R0 R1 b. Training Figure 3: PPO workﬂow. a.The program is given a state and is asked to perform an action and predict cumulative rewards. A move is sampled using a multivariate normal distribution, with the action prediction serves as the mean value, a decay variance for exploration, and the raw reward serving as a criterion for the move’s legality. b.Adaptive Laplacian Smoothing ([ Xi et al.(2021) ]) is used to optimize coordinates and recalculate each rewards. The loss functions are calculated using new rewards and the probability of sampling each moves to update parameters in NN. NN in FreeMesh-S are trained with an actor-critic framework, which was introduced in [ Konda and Tsitsiklis(2000) ]. We use a distributed version of proximal policy optimization (DPPO), an steady off-policy algorithm based on actor- critic. According to [ Schulman et al.(2017) ], PPO has outperformed A2C, A2C+Trust Region, and CEM on seven MuJoCo environments, including HalfCheetah-v1, Hopper-v1, InvertedDoublePendulum-v1, InvertedPendulum-v1, Reacher-v1, Swimmer-v1, and Walker2d-v1. Each of Nparallel actors collects Titime steps of data in each episode. Then we calculate each trajectory, shufﬂe them, and back propagate mini-batches with Adam ([ Kingma and Ba(2014) ]), forKepochs. In Figure 3, a workﬂow of Algorithm 1 is shown. Lines 3 to 9 are referred to as the sampling section a, while lines 11 to 20 are referred to as the training section b. Algorithm 1 Distributed Proximal Policy Optimization with clipping 1:forepisodei= 1;2;;Mdo 2: foragentj= 1;2;;Ndo 3: Get an initial state Sij 1(iandjare omitted hereafter); 4: fortime stept= 1;2;;Tdo 5: Collect value network prediction V(St;w) 6: Sample action At(jSt;); 7: Get the raw reward Rr tand next state St+1; 8: If done, break the loop; 9: end for 10: end for 11: Optimize coordinates with Adaptive Laplacian Smoothing ([ Xi et al.(2021) ]), and calculate optimized rewards Ro t; 12: Calculate the TD error t=Tt1P i=1i[Ro t+V(St+i+1;w)V(St+i;w)] 13: Calculate cumulative rewards Ut=TP i=1Ro t 14: formini-batch number k= 1;2;;Kdo(komitted hereafter) 15: forupdate epoch l= 1;2;;Ldo 16: Calculate policy loss J=P tminn new l(AtjSt;) old(AtjSt;)t;cliph new l(AtjSt;) old(AtjSt;);1;1 +i to 17: Calculate value loss JV= smooth _L1_loss [V(St;w);Ut] 18: Back propagate w; withJ+JV 19: end for 20: end for 21:end for In compared to business software meshes and the mesh in [ Pan et al.(2021) ], the parameters w; we acquire after DPPO are strong enough to give a state-of-the-art mesh. Many investigations have revealed, however, that higher-quality meshes have higher variance, because the goal of reinforcement learning is for the agent to learn an optimal, or 3 Generate Plane Quad Mesh with Neural Networks and Tree Search A P REPRINT Extraction1. Selection Q0+U0 Q1+U1 max2. Expansion & Simulation Q0+U0 Q1+U1 max Q0+U0 Q1+U1 1st Pπ3. Backpropagation Q max Q 1st QQ RepeatTraining Policy Network Prediction DPPOActor- Critic Network Predictiona c b Figure 4: MCTS workﬂow 1. a.Behavior cloning is a quick and easy technique to create effective policies given limited data. b.We train an actor-critic network with the DPPO algorithm and store the policy model, either from scratch or from a pre-trained model. c.In Selection, each simulation traverses the tree by selecting the edge with the highest sum of an expected quality Qplus an Upper Conﬁdence Bound (UCB) Ubased on a prior probability Pand visit count N for that edge. In Expansion and Simulation, the leaf node is expanded and evaluated by the quality of one meshing to the end. In Back Propagation, quality Qis utilized to update all evaluations Qin the sub-tree under the current episode. Once the search is complete, search times and expected qualities are returned, and the user can choose the next move. nearly-optimal, policy that maximizes the expected total reward, as described in [ Russell and Norvig(2010) ], the policy may never ﬁnd the ideal solution with a high variance. Thus, we combine NN with an MCTS to search for coordinates surrounding projected action, resulting in more dependable ﬁndings. 2.2 Monte-Carlo tree search in TreeMesh We use MCTS to perform a brute-force search to check if there are any better meshes. The prediction point will be the center of a 55dot array with a unit side length of 0:1. As shown in Figure 4 and Figure 5, the MCTS uses the policy network(jSt;)to guide its simulations. Each edge in the search tree holds a prior probability 0(AtjSt)with the same mean value as (AtjSt;). Each episode starts at the root and ﬁnds the child node that maximizes an upper conﬁdence bound before it reaches a leaf node: at= arg max At[Q(St;At) +U(St;At)] = arg max At U0+cPUCT0(AtjSt)pNt 1 +Nt+1(1) cPUCT is a hyper-parameter that determines the proportion of UCB to expected quality. U0is the cumulative reward from time step 0. If the leaf node has been visited, the algorithm will extend all possible moves and choose the ﬁrst child node. Finally, the algorithm will run all the way to the ﬁnish, updating nodes along the way. The whole pseudo code is in Algorithm 2. 2.3 State and action setting As shown in Figure 6 (a), the state tensor has 18 dimensions, including: 1. xP0xmin xmaxxmin;yP0ymin ymaxymin 2. xPr1xmin xmaxxmin;yPr1ymin ymaxymin 3.(x0 Pr2;y0 Pr2) 4.(x0 Pl1;y0 Pl1) 4 Generate Plane Quad Mesh with Neural Networks and Tree Search A P REPRINT Algorithm 2 Monte-Carlo tree search 1:function UCBS EARCH (s0) 2: Create a root node n0withN(n0) = 0;Q(n0) = 0;n’s parent NULL 3: while Within computational budget and ( 500-Episodeior conﬁrm continue) do 4:n SELECTION (n0) 5:n;Q =ROLLOUT(n) 6: UPDATE (n;Q) 7: end while 8: returnn0 9:end function 10:function SELECTION (n) 11: whilenis non-terminal ( N(n)1)do 12:n arg max n02children of n[Q(St;At) +U(St;At)] 13: end while 14: ifN(n) = 1 then 15: return EXPAND (n) 16: else 17:N(n) = 1 18: returnn 19: end if 20:end function 21:function EXPAND (n) 22: Addn02children ofnton’s children list 23: returnn’s 1st child 24:end function 25:function ROLLOUT(n) 26: while not done do 27:n uniform(n0s child) 28: end while 29: return noden, qualityQ 30:end function 31:function UPDATE (n;Q) 32: whilenis not NULL do 33:N(n) =N(n) + 1 34:Q(n) =Q(n)[N(n)1]+Q N(n) 35:n n’s parent 36: end while 37:end function 5.(x0 Pl2;y0 Pl2) 6.(x0 P1;y0 P1) 7.(x0 P2;y0 P2) 8.(x0 P3;y0 P3) 9.m N 10.n N As demonstrated in Figure 6 (b), from 1 to 6 are coordinates of Pr2;Pl1;Pl2;P1;P2;P3, standardized by coordinates of P0andPr1(by dragging P0toO, andPr1to(1;0)). The standardization procedure is as follows: "x0 y0 1# ="cossin0 sincos0 0 0 1#2 641p (xP0xr1)2+(yP0yr1)20 0 01p (xP0xr1)2+(yP0yr1)20 0 0 13 75"1 0xP0 0 1yP0 0 0 1#"x y 1# =E"x y 1#(2) 5 Generate Plane Quad Mesh with Neural Networks and Tree Search A P REPRINT Begin at S 0 at leaf go to the child node with max a tyes novisited noyesexpand and choose the first child node rollout update Figure 5: MCTS workﬂow 2. Notice that the ﬁrst "no" corresponds to step 1 in Figure 4, from the ﬁrst "yes" to "roll out" is step 2, and "update" relates to step 3. P0P1 Pr1Pr2 Pl1Pl2Pθ2 Pθ3 Pθ1 (a) State setting OXY θ P0Pr1 Pr2Pl1Pl2 (b) Coordinate transformation Figure 6: A state tensor with 18 dimensions and a sketch showing the coordinate transformation. The tensor includes 12 standardized coordinates, 4 re-scaled coordinates, and 2 progress indicators to aid in the evaluation of the value network. 3 steps are required to standardize coordinates: drag P0toO, rotateP0Pr1to x-axis, and re-scale Pr1to(1;0). There are a few notations that need to be clariﬁed. From 1 to 6, PliandPrj indicate the ith point on the left hand (clockwise) and the jth point on the right hand (anti-clockwise) of P0, respectively. These points are particularly relevant to the selection of point P1. All points on the current boundary except P0;Pl1;Pr1are traversed to ﬁnd P1;2;3. They are each in one of the three trisections of \Pr1P0Pl1and are close enough to point P0. (a composite function is used to sort out the three points). To speed up convergence, the three points are clipped to a circle with a centerP0and a radius of 3(in the transformed coordinate system). 7 and 8 are coordinates of P0andPr1, re-scaled by coordinates of the top right and bottom left corner points on the original boundary. 9 and 10 indicates the progress of meshing. mis the number of sides of the remaining polygon, nis the number of meshed quadrilaterals, and Nis the maximum step permitted in meshing. The actor-critic network outputs a value and an action. The action is a 2-dimensional tensor of the next point (x0;y0), and the value V(St;w) =E(UtjSt)is an expectation of cumulative rewards. Because the minimum angle for all boundaries is smaller than , both coordinates are constrained in [0;1:5]with a hyperbolic tangent function, the new point is most likely in the ﬁrst quadrant. An inverse matrix E1will be used to transform the new point back to the previous coordinate system. 6 Generate Plane Quad Mesh with Neural Networks and Tree Search A P REPRINT P0P1 P0P1P1 P0 (a) Three possible moves θ aA θ aA (b) Variables considered in judge- ment Figure 7: At any given time step, there are three possible moves. All of the moves generate a new quadrilateral, with the exception of the ﬁrst, which has no effect on the number of edges, and the other two, which reduce the number of edges by two. As a result, the starting boundary should have an even number of nodes. 2.4 Reward setting The reward function is a critical signal for the NN since it allows the NN to determine whether a movement is positive or harmful. We propose an improved reward function upon [ Pan et al.(2021) ]. This option delivers superior convergence performance, according to empirical analysis. Noticing that on Algorithm 1 line 7 and line 11, we calculate a raw rewardRr tand an optimized reward Ro t. The raw reward, with a range f5g[[0;1], is purely for judging a movement in Figure 7 (a), or give a punishment reward (e.g. in Figure 7 (a), the third move is illegal, the other two moves are legal): Rr t=8 >>< >>:5;illegal 4vuut4Y i=1" lipAtsgn(pAtli)# 1jaj 2j  2 ;legal(3) At time step t, the three raw rewards will be multiplied by their respective factor2Q i=1 min(k; 3)  3 , and the highest-valued move will be chosen. It’s worth noticing that we can also allow the NN to decide which action to take, but this will make convergence more difﬁcult. On line 11, we have collected a trajectory Rr 0;Rr 1;;Rr T. Assume there are a total of T0potential steps (with the remainingTT0rewards being -5). Then, using Adaptive Laplacian Smoothing ([ Xi et al.(2021) ]), we optimize all internal coordinates and update every effective reward to Ro t=Rr0 t T0.is a hyper-parameter with a value of 500that controls the cumulative rewards in the [100;1000] range. The penalty for illegal rewards remains at 5. 3 Result 3.1 Some DPPO results Here, we exhibit our work on two boundaries: a " 11!5" seed-density-changing rectangle, and a pentagonal plate that was studied in [Pan et al.(2021)]. Table 1 compares TreeMesh (DPPO algorithm used only) with ABAQUS on a " 11!5" seed-density-changing boundary. Equation 3 is used to compute the mean quality of the optimized mesh. We outperform ABAQUS in terms of quality by 13:8573% . Table 2 compares TreeMesh (DPPO algorithm used only) with ABAQUS and FreeMesh-S on a pentagonal boundary. We outperform ABAQUS in terms of quality by 5:3204% , and FreeMesh-S by1:4822% in quality, and by around 90% in time. 7 Generate Plane Quad Mesh with Neural Networks and Tree Search A P REPRINT Table 1: A " 11!5" seed-density-changing boundary with 4 MLP layers. Program Raw Mesh Optimized Mesh Quality Quantity Time ABAQUS 0:7036 22 1s TreeMesh 0:8011 34 103s 3.2 A DPPO+MCTS time line A more detailed experiment is conducted on the ﬁrst boundary (Table 1), with a 3-residual-block architecture. In 1500s, we output the raw mesh created by the DPPO algorithm. Weight ﬁles at three time points are rendered to execute a brute force MCTS. At this boundary, a better mesh is discovered and shown in Figure 8’s upper right corner. It has been appraised as a state-of-the-art mesh by some engineers. 3.3 Empirical Analysis of DPPO training reward in Table 1 and Table 2 We continued the training for approximately 4 hours. Figure 9 shows the smoothed training reward on both boundaries. The co-variance matrix of the multivariate normal distribution starts at 0:05Iand decays at a constant rate of 0:999per 100episodes during the training. This slows convergence slightly but allows the policy NN to explore more policies. The policy NN will eventually stop at a nearly-optimal policy if the initial learning rate is properly set ( 106) and decays at the same rate. Here, however, we remain the learning rate constant ( 106) to accelerate convergence.