DEVELOPMENT OF A METHOD FOR TRIANGULATION OF INHOMOGENEOUS REGIONS REPRESENTED BY FUNCTIONS

елементів, інцидентних вузлах, координати яких було змінено, перевіряєть-ся виконання умови Делоне і за необхідності виконується операція зміни діагоналі «flip». Після видалення


Introduction
Designing new equipment with rational characteristics often employs the structures consisting of heterogeneous domains (for example, multilayer structures, composite materials, etc.). The process of designing such structures is associated with the need to numerically analyze their performance. Consequently, that necessitates the construction of appropriate discrete models (meshes).
It is possible to account for differences in the characteristics of domains in discrete models by using a macro-or micro-approach. Under a macro-approach, a discrete model must adequately represent the geometry of the structure, while accounting for the heterogeneities employs different models for averaging (homogenization) the characteristics of domains, which typically leads to compromised accuracy. A micro-approach is based on the explicit modelling of heterogeneities at domains via the elements of the discrete model, making it possible to improve modelling accuracy; however, it complicates the process of discretization. The first step in building a discrete model of an object is to describe its shape. Functional representation is one of the most universal ones to represent the shape of objects. Such a representation makes it possible to describe the shape of an object of arbitrary complexity using implicit functions.
An implicit function representing a structure's shape can be derived constructively using the R-functions theory [1,2]. Such a function makes it possible to check if an arbitrary point in space belongs to the structure.
Objects of constant thickness can be regarded as two-dimensional. The construction of discrete models for the shapes of two-dimensional objects most frequently employs triangular elements. Their choice is predetermined by the availability of simple effective procedures that make it possible to build high-quality meshes. At the same time, existing methods [3][4][5][6][7][8] do not make it possible to triangulate the heterogeneous regions represented functionally. Therefore, it is a relevant task to construct methods for triangulation considering the geometric patterns and boundaries of domains, as well the boundaries of heterogeneous subdomains.

Literature review and problem statement
The Delaunay triangulation is one of the most common procedures for constructing discrete models. Such a triangulation warrants that the circle that enclosures an arbitrary triangle does not capture those nodes that do not belong there. One of its properties is the possibility to derive models without «bad» triangles (with angles that are too sharp). There are effective methods for constructing the Delaunay triangulation with constraints for the regions assigned by polygons [3]. The parallel Delaunay triangulation methods are being actively developed that have constraints for domains set by planar graphs [4][5][6]. The parallel Delaunay triangulation with constraints implies «ear clipping» a polygon [4]. To increase the rate of triangulation, there are protocols for a simultaneous insertion of nodes and edges into initial triangulation [5]. There are the algorithms developed for a graphic processor to enable the faster execution of common operations [6]. However, representation of domains of complex shape in the form of polygons or planar graphs is not always convenient. At the same time, transition from the functional to the boundary representation in the form of a planar graph is a separate task.
Current methods for building discrete models of the domains assigned functionally [7,8] make it possible to construct models that rather accurately approximate boundaries. It is possible, in order to triangulate an arbitrary region, set functionally, to apply a background mesh. One eliminates all outer nodes and elements in a background mesh, and then fills the near-boundary layer of space [7]. For parallel triangulation, a background mesh can be divided into pairs of disjoint parts [8]. In methods based on the background mesh, positions of inner nodes are defined by the coordinates of nodes at the background mesh. As a result, when simulating heterogeneous regions, the geometric features of subdomains might not be approximated by the internal nodes and edges.
Paper [9] constructed a method to triangulate surfaces that are assigned by implicit functions. The authors suggested minimizing the error in the surface representation to search for geometric patterns. However, a Laplace smoothing, used in that work, weighted for the areas of adjacent elements, does not necessarily eliminate triangles with sharp angles.
Studies [10,11] suggested using the hierarchical structures (quadtree) to split the two-dimensional regions, set by using a boundary representation (brep), into quadrangles, taking into account the geometrical features of subdomains. Under a given method, the background cells that overlap the subdomains, are recursively divided into parts. Next, one performs the nodes' coordinate correction so that there are no sections of boundaries that are close to them. At the final step, one finds the points of intersection between edges and boundaries, and recovers the topological mesh correctness. The method was generalized for a three-dimensional case [12] to build meshes of hexagons. Its main disadvantage is the large number of nodes and elements near the boundaries (a consequence of the application of hierarchical structures), as well as a possibility to obtain blunt angles in elements (larger than 160°).
Paper [13] devised a method for constructing discrete models taking into account the directions of composites' orthotropy using micro-tomography. The disadvantage of this method is the impossibility to use it for heterogeneous domains, which include disparate subdomains.
Thus, our analysis of publications suggests that there are no identified effective methods to triangulate the functionally assigned domains considering the heterogeneity of the materials.

The aim and objectives of the study
The aim of this study is to devise a method for the triangulation of functionally assigned heterogeneous domains, whose application would result in that that boundaries of inner subdomains are approximated by edges of the inner elements. That would improve accuracy of modeling the performance of inhomogeneous domains by taking into account the shape of subdomains representing different materials.
To accomplish the aim, the following tasks have been set: -to construct an algorithm for the initial triangulation of a functionally-assigned domain; -to construct an algorithm for the adaptation of initial triangulation to the boundaries of subdomains representing different materials; -to verify the constructed algorithms.

Functional approach to representing the shape of inhomogeneous structures
In a general case, one can assume that an arbitrary two-dimensional domain Ω is assigned functionally if there is a defined function F that is greater than zero only at the inner points of Ω. In this case, function F equals zero at the boundary of Ω and is less than zero at outer points of Ω. Accordingly, the following set is defined: Domains of a simple shape can be described using elementary functions. For complex regions, one can build the compositions of implicit functions using the operations of disjunction, conjunction, and negation. Such operations are proposed in the R-functions theory [1,2] in the form of systems, one of the most common of which is: where f 1 and f 2 are the values for implicit functions. For example, a rectangular domain can be represented as a conjunction of two bands: To represent the geometric structure of a heterogeneous domain that includes heterogeneous subdomains, it is possible to apply the system of sets in form (1): where n is the number of subdomains; w 0 is the implicit function that represents the initial domain Ω 0 ; w 1 , …, w n are the implicit functions representing the subdomains Ω 0 , …, Ω n , in this case, Ω 0 ∩ Ω i ≠ ∅, if 1 ≤ i ≤ n, and Ω i ∩ Ω j = ∅, if 1 ≤ i ≤ n, 1 ≤ j ≤ n and i ≠ j. That is, the first function describes the domain, and the rest -sub-domains.
Hereafter, we assume that a set is assigned if the appropriate implicit function is defined. Accordingly, a heterogeneous domain with subdomains will be defined by the system of functions. For example, system: describes a rectangular region of width w, height h, centered at point (x 0 , y 0 ), which includes a subdomain ( Fig. 1) located above the diagonal with top (x 0 , y 0 ). Thus, a system of implicit functions can be applied to describe the shape of an arbitrary heterogeneous structure and its subdomains.

Algorithm for constructing the initial triangulation of a domain
Let some two-dimensional inhomogeneous domain be assigned in form (2). Triangulation of a heterogeneous region, taking into account the subdomains, can be performed using a background mesh of triangles. A background mesh of triangles, for example, might be structured, but it should fully incorporate the domain.
The basic idea of the algorithm implies the iterative displacement of nodes, closest to the boundaries of domain and subdomains, by validating the Delaunay condition and by optimizing coordinates for neighboring nodes. The method's algorithm is as follows.
1. Generate for a background mesh a set of edges that are intersected by the domain's boundary Ω 0 (the roots of function Ω 0 are available): for which the distance from one of its vertices to the point of intersection of this edge by boundary Ω 0 is minimal.
2. 2. Displace the closest node to the point of intersection between edge e j and boundary Ω 0 .
2. 3. For each node v e , which is adjacent to any node at the vertices of edge e j while w 0 (v e ) ≠ 0, the coordinates are computed as a result of minimizing the functional: where is the normalized vectors, which are defined along the sides of triangles, incident to node v e , τ is the parameter, which, to account for the size of the model, can be considered equal to 1 α (α is the mean area of a triangle).
2. 4. For the triangles incident to the nodes processed at step 2. 3, check meeting the Delaunay condition and, if necessary, perform the procedure «flip» to shift the diagonal.
2. 5. Renew the set of edges at discrete model E.
3. Remove those nodes and elements that are outer relative to Ω 0 .
4. For each boundary node v b (w 0 (v b ) = 0), perform the smoothing by minimizing the following functional: where κ is the method parameter; δA is the set of boundary nodes adjacent to v b ; projection(p, w 0 (x, y)) is the projection of some point p onto the boundary of the domain, assigned by function Ω 0 .
5. The end of the algorithm. When searching for the edge's vertices closest to the boundary at step 2. 1, in order to improve the speed of algorithm execution, one can use a linear approximation.
At step 2. 2, in order to find the intersection of boundary Ω 0 by the edge, one can apply a binary search (function Ω 0 at the edge's vertices accepts values with different signs).
For each node at step 2. 3, a set of triangles adjacent to it is considered. For each triangle, one determines the sum of exponents of the vector products of vectors that are defined at its vertices (Fig. 2). For each i-th vertex, value of s i = a i × b i (i = 1, 3) would be positive if the angle at this vertex varies from 0° to 90°, and negative -if the angle varies in the range from 90° to 180°. As shown in [14], the minimization of functional (4) makes it possible to eliminate degenerate elements and ensure an even distribution of angles in the elements.

Fig. 2. Scheme for constructing vectors at the sides of a triangle
At step 2. 4, one applies the operation «flip» to change a common edge between triangles [9]. If the circle that enclosures one triangle includes all vertices form the neighboring vone, their shared edge is changed (Fig. 3). The «flip» operation makes it possible to ensure that the Delaunay criterion is met locally.

Fig. 3. The «flip» operation that changes the shared edge
At the fourth step of the algorithm, smoothing the boundary nodes employs the minimization of the functional, which takes into account the distance to the boundary section approximated by the edge. In order to find a minimum of this functional, it is possible, for example, to use a conjugate gradient method.
The result of the algorithm execution is the built initial discrete model of a domain.

Algorithm for the adaptation of initial triangulation to the boundaries of subdomains
Adaptation of the initial triangulation to the boundaries of subdomains can be divided into two stages. The first step implies performing the adaptation of boundary edges to the projection of intersection points onto a domain's boundary. The second stage implies the adaptation of inner edges. Thus, on can use the following algorithm.
2. While B i ≠ ∅: 2. 1. For edges from set B i -follow steps 2. 1-2. 5, described for constructing the initial triangulation. In this case, at step 2. 2, find the projection of the intersection point between edge e j and boundary Ω i onto a domain's boundary Ω 0 , next, move the closest node.
5. Perform the smoothing of the nodes' coordinate at the subdomain's boundary.
For each boundary node v b from subdomain Ω i (w i (v b ) = 0), perform the smoothing by minimizing the following functional: where κ is the method parameter; δA i is the set of boundary nodes from subdomain Ω i ; adjacent to v b ; projection(p, w(x, y)) is the projection of some point p onto the domain's boundary set by function w. In this case, one searches for the projection of the point onto the boundary of intersection between domain Ω 0 and subdomain Ω i . II. The end of the algorithm.

Results of applying the method of triangulation of inhomogeneous domains
It is possible to use, for the triangulation of a non-uniform rectangular domain, assigned by formula (3), the structured background mesh (Fig. 4). For example, one can use a background mesh built for a square of 4 × 4, which contains 121 nodes and 200 elements (Fig. 4, a, color shows the distribution of values for function w 0 in formula (3)). The initial triangulation would involve 95 nodes and 156 elements (Fig. 4, b). The result of the adaptation to the boundary of the sub-domain takes the form shown in Fig. 4, c.
where R is the radius of the circle described around a hexagon; r is the radius of the circular subdomain; function: , , , , cos , sin If one uses the unstructured grid that contains 210 nodes and 333 elements as the initial one (Fig. 5, a, color shows the distribution of values for function w 0 ), the result of adaptation would take the form shown in Fig. 5, b. The friction disk of a transmission assembly is a complex domain with a hole and ten round subdomains (Fig. 6). Its shape can be assigned by an implicit function describing the point symmetry of the cyclic type: where S, R, r, h, c and F are the dimensions, indicated in the drawing (Fig. 6, a); n is the number of teeth; clutch(x, y, S, R, r, h, n) is the function of form: clutch , , , , , , is the partial sum of a Fourier series [2]; θ = arctg 2 (y, x) is the arc tangent of ratio y x. The result of triangulating the domain assigned by function (6), when using a structured mesh, is shown in Fig. 6, b.  Table 1 gives results from studying the influence of the number of nodes and elements in the initial (structured) mesh on the value for a minimum angle in the element. One can notice that the value for a minimum angle less than 18° was not obtained. This is due to the local optimization of angles and to the application of operation to change a shared edge.
Thus, the constructed method makes it possible to build discrete models of heterogeneous domains that are functionally assigned. In the resulting discrete models, the boundaries of a domain and the sub-domains are approximated by the elements' edges. Table 1 Influence of the number of nodes and elements in the initial mesh on the value for a minimum angle in element at triangulation

Discussion of results of applying the method of triangulation of inhomogeneous domains
As the result of the application of the devised method, in the derived discrete models the boundaries of subdomains in inhomogeneous domains (in contrast to results reported in [7,8]) are explicitly represented by the elements' edged. As opposed to papers [10][11][12], the devised method does not require the thickening of discrete models at the boundaries of subdomains.
The result of the application of the method for triangulating heterogeneous domains is the derived discrete models whose elements include angles that are less than 20°. Applying the Delaunay triangulation methods with constraints [3][4][5][6] makes it possible to obtain models whose elements' angles are not less than 26°. However, the models of boundaries of a non-uniform domain and its subdomains, built by using the devised method, could be used as source data (in the form of planar graphs), for the triangulation methods described in [4][5][6].
The accuracy of approximation and characteristics for the resulting discrete model depend on an initial grid. It is obvious that the size of an element in the initial grid must be less than the geometric patterns of the simulated domain.
The disadvantage of the proposed method might be its ite rative nature. Consequently, its implementation for parallel computer systems is problematic.
When smoothing the coordinates of nodes at the boundaries of a f domain and the subdomains, one can apply the method proposed by authors in [15].
The method reported in the current paper could be used as a preprocessor to the finite element analysis of inhomogeneous regions. The prospect of further advancement is to generalize the method for cases of three-dimensional bodies or shell structures.

Conclusions
1. We have constructed an algorithm for the initial triangulation of a functionally-assigned domain. A given algorithm employs a consistent correction of coordinates for the nodes at initial triangulation. Initial triangulation can be arbitrary, but it must fully capture the structure. At each step of the algorithm, the closest node is displaced to the boundary of the structure or the subdomain. Following the displacement of each node, coordinates for the adjacent nodes are determined by minimizing the functional of exponents of doubled areas considering the sign. In this case, for the elements, incident at nodes with modified coordinates, meeting the Delaunay condition is checked and, if necessary, the operation «flip» is executed to change the diagonal. The result of such an algorithm is the discrete model of a domain, in which the triangles satisfy the Delaunay condition.
2. An algorithm for adapting the initial triangulation to the boundaries of subdomains representing different materials has been constructed. Underlying the algorithm is the consistent displacement of nodes to the boundaries of subdomains. The first phase implies the displacement of boundary nodes from the initial triangulation. In this case, one searches for the point of intersection between a domain and the subdomains. The second phase implies the displacement of inner nodes. Both phases involve the local smoothing of nodes adjacent to the displaced one and checking if the Delaunay condition is satisfied. The result is that the subdomains are represented by elements of a discrete model, which is needed to improve the accuracy of modeling heterogeneous domains.
3. Verification of the constructed algorithms has shown that the resulting discrete models lack «bad» triangles (with very sharp angles). Uniform and non-uniform discrete models could be used as the background models.
The devised method for the triangulation of inhomogeneous functionally-assigned structures makes it possible to build discrete models, in which the subdomains' boundaries are represented by the elements' edges. From an applied point of view, such a representation could improve the accuracy of a finite element analysis of inhomogeneous structures. A limitation of the devised method is the requirement that in the background triangulation a mesh size should be less than the size of the smallest geometrical feature of the structure.
It should be noted that the functional approach makes it possible to simulate heterogeneous domains of arbitrary complexity. In this case, the shape of domains and subdomains can be changed through the variation of model parameters. The result is the simplified investigation and optimization of the shapes of domains in the design process.