1. Introduction

Widely used mesh-based numerical methods such as finite element, finite difference, and finite volume methods, introduce a finite number of nodes to specify boundary conditions and perform numerical computations, and use spatial grids to approximate the geometric shape of a model. However, in modelling problems with complex geometry, difficulties often appear in creating good spatial grid that conforms to the shape of the model. To overcome this obstacle, a new class of numerical methods has been developed called meshfree or meshless methods. These methods may still use spatial grids to construct the basis functions and perform numerical computations, but such grids do not necessarily have to conform to the geometric model. To date, many different meshfree methods have been developed. Their detailed review and comparison can be found in many references [1], [2], [3], [4], [5]. Many meshfree methods use radial basis functions to represent solutions of engineering problems [6]. Using meshfree methods significantly simplified the meshing process, but, at the same time, it made the treatment of boundary conditions as demanding task [7].

Here, a numerical method for solving engineering problems that enables exact treatment of all prescribed boundary conditions at all boundary points and does not require numerical integration is presented. It combines meshfree method known as solution structure method, atomic basis functions (ABFs) and a collocation technique.

In the solution structure method, a solution is sought in the form of formulae called solution structure. The original idea is due to Kantorovich [8]. He proposed that the homogeneous Dirichlet conditions may be satisfied exactly by representing the solution as the product of two functions: (1) an real-valued function that takes on zero values on the boundary points; and (2) an unknown function that allows to satisfy (exactly or approximately) the differential equation of the problem. Such a solution structure was used by Kantorovich and his students to solve boundary value problems on geometrically simple domains.

Rvachev [9] suggested using R-functions--the real valued functions that behave as continuous analogy of logical Boolean functions. R-functions allow construction of a set of functions vanishing on the boundary that possess desired differential properties and may be assembled into a solution structure. Based on the theory of R-functions, the RFunction Method (RFM), also known as solution structure method, is developed which can be applied to problems with arbitrarily complex domains and boundary conditions. Over the last several decades, the theory of R-functions have been applied to numerous scientific and engineering problems by Rvachev and his students [10], [11]. The RFM has been applied to problems of thermo conduction, elasticity, magneto- hydrodynamics, various problems in inhomogeneous media, and many other areas [12].

This paper presents the use of atomic basis functions to approximate unknown differential component of the solution structure. They are infinitely-differentiable functions with compact support [13], [14], [15]. Rvachev and Rvachev [13], in their pioneering work, called these basis functions "atomic" because they span the vector spaces of all three fundamental functions in mathematics: algebraic, exponential and trigonometric polynomials. In numerical modelling, we applied Fup basis functions that belong to the atomic functions of algebraic type [16], [17], [18]. All derivatives of atomic Fup basis functions required by differential operators in the solution structure can be used directly in the numerical procedure. This fact allows to use procedures based on strong formulation. To determine the coefficients of linear combination in the solution structure, a collocation technique is used.

2. Solution Structure Method: Basic principle

The original idea of the solution structure method [8] is to express the solution of two-dimensional boundary value problem with homogeneous Dirichlet boundary conditions

u[|.sub.[GAMMA]] = 0 (1)

by formula called solution structure in the form of the product of two functions:

u = [omega][PHI] (2)

where [omega]: [R.sup.n] [right arrow] R is a known function that takes on zero values on the boundary of the domain [GAMMA] and is positive in the interior of the domain [OMEGA], and [PHI] is some unknown function that allows to satisfy (exactly or approximately) the differential equation of the problem.

In most practical situations, unknown [PHI] is represented by a linear combination of basis functions

[PHI] = [n.summation over (i=1)][C.sub.i][F.sub.i] (3)

where [C.sub.i] are scalar coefficients and [F.sub.i] are some basis functions. The solution structure does not place any constraints on the choice of basis functions. Numerical values of the coefficients [C.sub.i] can be obtained by using different numerical methods.

For complex domains, Rvachev [9] set the theory of R-functions which was the basis for the development of the Rfunction method.

2.1. Theory of R-functions

The question is how to generate functions that will simply describe the given domain and satisfy different boundary conditions on any part of the boundary? There is an elegant way by using [omega] functions.

Generating [omega] functions over the complex area is proposed by the Ukrainian scientist V. L. Rvachev. Rvachev [9] came up with the idea that logic operations of Boolean algebra are applied to the functions. In this way he created so called semi-algebra. Basic R-operations are shown in Fig. 1.

R--conjunction (section) ([f.sub.1] [conjunction] [f.sub.2]):

Boolean function is logical "and" ([conjunction]). It is defined in the form: [mathematical expression not reproducible]

R--disjunction (union) ([f.sub.1] [disjunction] [f.sub.2]):

Boolean function is logical "or" ([disjunction]). It is defined in the form: [mathematical expression not reproducible]

R--negation ([logical not]f = -f). The logical negation of the function is the change of sign of that function: [F.sub.3](f) = [bar.f] [equivalent to] -f.

Using these operations, we can determine the function [omega] over the very complex domains. Then such functions are called R-functions.

3. General solution structures

Once we constructed the solution structure u = [omega][PHI] for the boundary value problems with homogeneous Dirichlet boundary conditions, it is easy to obtain the solution structure for nonhomogeneous conditions

u[|.sub.[GAMMA]]= [[phi].sub.0] (4)

Let [phi] be the extension of the [[phi].sub.0] inside the domain [OMEGA]. Then the solution structure

u = [omega][PHI] + [phi] (5)

satisfies the prescribed boundary conditions exactly. In practice, the function [[phi].sub.0] may be specified in a piecewise fashion, with a different value [[phi].sub.i] prescribed on each portion of the boundary [[GAMMA].sub.i]. Such individual boundary conditions may be combined into a single global function [phi] [9], [12]:

[mathematical expression not reproducible] (6)

Detailed and systematic derivations of solution structures for more general boundary conditions can be found for example, in references [9] and [12]. General solution structure for the second- order boundary value problem with mixed boundary conditions:

[mathematical expression not reproducible] (7)

can be written in the form that interpolates boundary conditions on [[GAMMA].sub.1] and [[GAMMA].sub.2] according to [9]:

[mathematical expression not reproducible] (8)

where h is a part of the Robin boundary condition. [[omega].sub.1] is part of the primary function [omega] that belongs to the part of the boundary [[GAMMA].sub.1] while [[omega].sub.2] is part of the primary function [omega] that belongs to the part of the boundary [[GAMMA].sub.2]. For example, for the domain in Fig. 2 can be written:

[mathematical expression not reproducible].

Differential operator of the first order [mathematical expression not reproducible] is:

[mathematical expression not reproducible] (9)

Differential operators in the solution structures transmit information about the dynamic boundary conditions from the boundary to the domain. To this extension be consistent, the primary function ra must be normalized function [11]. From the general solution structure, one can derive solutions for all cases of boundary conditions. For h=0, Robin boundary condition becomes Neumann boundary condition:

[mathematical expression not reproducible] (10)

For these mixed boundary conditions, the solution structure is:

[mathematical expression not reproducible] (11)

For only Neumann boundary conditions, with substitutions [[omega].sub.1] = 1, [[omega].sub.2] = [omega], the solution structure can be written in this form:

u = [PHI] - [omega][D.sup.[omega].sub.1]([PHI]) - [omega][psi] + [[omega].sup.2][PHI] (12)

where [D.sup.[omega].sub.1]([PHI]) = [nabla][omega] x [nabla]([PHI]) is a differential operator in the direction of the internal normal to the boundary [GAMMA]. The resulting function [psi] interpolates the individual values [[psi].sub.i]:

[mathematical expression not reproducible] (13)

For only Dirichlet boundary conditions u[|.sub.[GAMMA]] = [[phi].sub.0], with substitutions [[omega].sub.1] = [omega], [[omega].sub.2] = 0, the solution structure can be written in the previously mentioned form u = [omega][PHI] + [phi].

4. Combining solution structure method with atomic basis functions and collocation method

A physical field being modelled is represented by a solution structure that can be written in a general form:

u = u([omega], [PSI], [PHI]) (14)

The solution component [omega] exactly introduces all information on the domain geometry. So, [omega] is called the primary function of the solution. The second component [PSI] accurately introduces all information about boundary conditions. The third or differential component [PHI] can be determined as a linear combination of chosen basis functions in the form [PHI] = [n.summation over (i=1)][C.sub.i][F.sub.i]. Approximation properties of the solution structure to a large extent are determined by selection of basis functions [{[F.sub.1]}.sup.n.sub.i=1]. Here, we propose to use [Fup.sub.n](x) basis functions because they are well conditioned, they have compact support and good approximation properties [19]. Index n denotes the highest degree of the polynomial that can be expressed exactly in the form of linear combination of [Fup.sub.n](x) basis functions. The length of the Fup function support is determined according to expression [-(n + 2) [2.sup.-n-1]; (n + 2) [2.sup.-n-1]].

For the Fup basis functions, a criterion of choice of collocation points exists [20]. It is optimal to perform collocation in natural knots of basis functions, i.e. vertices of basis functions situated in a closed domain as shown in Fig. 3. This selection of collocation points provides the simplest numerical procedure, banded collocation matrix is obtained, which is diagonally dominant and thus well conditioned. This selection also implies uniformly distributed nodes set.

To basis functions set can be complete, we must keep all basis functions with vertices outside the domain that have values inside the domain different from zero, and we have to write one conditional equation for each of them.

In the case of 1D problems, if [Fup.sub.4](x) basis functions are selected, an unknown component of the solution is:

[mathematical expression not reproducible] (15)

Then the function [PHI](x) is defined on the whole real axis. This linear combination can accurately represent an arbitrary algebraic polynomial of the fourth degree [P.sub.4](x). The R-function method satisfies exactly all boundary conditions using the solution structure, so it is only necessary to satisfy the differential equation in collocation points inside the domain [a, b].

For the coefficients of basis functions with vertices outside the domain (see Fig. 4), the following recursive formulas that represent connections between external and internal coefficients are used:

[mathematical expression not reproducible] (16)

These recursive formulas are obtained from the condition that the fifth derivative in the middle of the collocation points is equal zero ([[PHI].sup.V](x) = 0 in x = a + [DELTA]x/2 and x = b - [DELTA]x/2, see Fig. 4). In this way, an arbitrary solution function outside of the domain is naturally extended by corresponding polynomial of the fourth degree [P.sub.4](x).

The basis function for numerical analyses of two-dimensional problems is obtained from the Cartesian product of two one-dimensional Fup functions defined for each direction:

[Fup.sub.n](x, y) = [Fup.sub.n](x) x [Fup.sub.n](y) (17)

Calculation of all required derivatives of function [Fup.sub.n](x,y) can be written in an analogue form [16].

5. Numerical examples

In this section, we will investigate numerical properties of the proposed approach by solving boundary value problems for Poisson equation with homogeneous Dirichlet boundary conditions and nonhomogeneous Neumann boundary conditions and compare the obtained results with analytic solutions. To perform this comparison, we will choose a benchmark problem with known analytic (exact) solutions.

5.1. Example of homogeneous Dirichlet problem: torsion problem

We now apply R-function method to the torsion problem for a bar with the cross section shown in Fig. 5a). The elastic torsion of a bar is a classical problem in the theory of elasticity [21], [22]. This problem may be reduced to the boundary value problem with Poisson equation and homogeneous Dirichlet boundary conditions:

[mathematical expression not reproducible] (18)

where u(x, y) is the stress function, G is the shear modulus, while v is the angle of twist per unit length of a bar. Shear stress components are determined according to the following expressions:

[[tau].sub.xz] = [partial derivative]u/[partial derivative]y; [[tau].sub.yz] = - [partial derivative]u/[partial derivative]x (19)

Exact solution for this domain has been derived by algebraic polynomials in [23] and will be used here for comparison. An analytic solution for the stress function and maximum shear stress are expressed in term of a parameter a:

[mathematical expression not reproducible]. (20)

For G = 1.0, I = 1.0 and a = 12.0 an exact values are [u.sub.A exact] = 10.666, [[tau].sub.max exact] = 6.0.

We will represent approximate solution in the form (2), with function ra defined by (21). For the undetermined function [PHI] we choose a linear combination of [Fup.sub.4](x,y) basis functions on a uniform 26x29 grid shown in Fig. 5b). Only points within the domain are collocation points; for points outside of the domain (black dots in Fig. 5b) recursive formulas are written according to (16) while the remaining points (white dots in Fig. 5b) are not included in the procedure. So, the number of equations to be solved is less than 26 x29. Numerical solution obtained with this grid is [u.sub.A] = 10.663, [[tau].sub.max] = 6.003.

[mathematical expression not reproducible] (21)

Practically good enough numerical solution is obtained on a grid that has a minimum number of collocation points within the domain that allows the implementation of the procedure of solution structure method with the selected basis functions. The minimum number of collocation points means that within a domain there is a core of 5 x 5 points which allows writing recursive equations (16). It defines uniform 14 x 15 grid with 38 collocation points and a total of 176 equations. This grid gives numerical solution [u.sub.A] = 10.646 which is 98.81% of the exact solution. Such accuracy is expected because the exact solution (20) is a cubic polynomial and linear combination of [Fup.sub.4](x,y) basis functions can expressed exactly an algebraic polynomial of the fourth degree.

Figures 6a) and 6b) show the computed shearing stresses [[tau].sub.xz] and [[tau].sub.yz] for geometric domain with a = 12.0.

5.2. Dirichlet and Neumann boundary value problems

Let consider a bar with square cross-section length of sides 2d = 10 cm, G = 1.0 kN/[cm.sup.2], v = 1, Fig. 7a).

For this cross-sectional shape, arbitrarily exact solution can be found using the development of the stress function into infinite trigonometric Fourier series. According to the analytical expressions for maximum value of the stress function and maximum shear stress [21], the following values are obtained:

[mathematical expression not reproducible] (22)

First, the problem (18) was analyzed by proposed method using solution structure (2) with a different density of uniform collocation points which make up a grid that covers the given two- dimensional domain but does not conform to it. Fig. 7c) presents stress function over the domain. The number of collocation points is varied from a minimum required number 25 (Fig. 7b) to maximum number 2401. Numerical solution converges to the exact solution with the increase in the number of points N x N for N=5,7,13,25,49 as shown in Fig. 8.

The elastic torsion of a bar can be described also in formulation of warping function. This problem may be reduced to the boundary value problem with Laplace equation and nonhomogeneous Neumann boundary conditions:

[mathematical expression not reproducible] (23)

where w is the warping function, n is the outer normal to the boundary [GAMMA], [n.sub.1] and [n.sub.2] denote components of the outer normal in x and y directions, respectively. Shear stress components are determined according to the following expressions:

[[tau].sub.xz] = Gv([partial derivative]w/[partial derivative]x - y); [[tau].sub.yz] = GI([partial derivative]w/[partial derivative]y + x). (24)

Now, the problem (23) is analyzed by proposed method using solution structure (12). The resulting function [psi] that interpolates the individual values [[psi].sub.i] according to (13) transmit information about the dynamic boundary conditions from the boundary to the domain. Fig. 9 presents warping function obtained by 31x31 grid with 625 collocation points.

Numerical value for maximum shear stress converges to the exact solution with the increase in the number of points NxN for N=7,9,13,19,41,49 as shown in Fig. 10. Thus, we can obtain the numerical solution of the torsion problem in two ways, i.e. by means of the formulation using the stress function and by means of the formulation using the warping function. It is very useful that the exact solution is securely between these two numerical values, as can be seen in Fig. 8b) and Fig. 10.

5.3. Example of complex cross section

To demonstrate the applicability of the proposed method combined of solution structure method, atomic basis functions and a collocation technique, we analyzed practical engineering torsion problem for a rod with the cross section shown in Fig. 11a). This is a textbook problem [21] with many good approximations already known. For example, an approximate analytic expression for torque in terms of parameters r, a, and b has been derived for the same domain in [24] and will be used here for comparison.

The primary function of the solution [omega] which exactly introduces all information on the domain geometry is constructed using R-operations described in the part 2.1. Figure 11.b) shows the contour of the function [omega] which represents a surface "inflated" over the domain. The numerical solution of the Poisson equation (18) was obtained by using uniform 55x55 grid (1217 collocation points).

By using approximate analytic expression for this cross-sectional shape [24], we verified our numerical model. For the fixed values r, a and b (Fig. 11a), the computed maximum value of shear stress component [[tau].sub.yz] is 11.819 while the value predicted by the closed-form expression in [24] is 11.956. Figure 12 shows the diagram of the computed shearing stresses [[tau].sub.yz] along the horizontal axis of symmetry.

6. Conclusion

This paper presents numerical method created to solve the modelling of problems with complex geometry in a way that spatial grids, that are needed to construct the basis functions and perform numerical computations, do not have to conform to the geometric shape of the model and at the same time boundary conditions can be well treated.

The basic principle of R-function method is described. We have combined RFM with atomic basis functions and get a method that connects the advantages of solution structure method, collocation technique and good approximation properties of atomic Fup basis functions.

Numerical examples considered here have shown that properly constructed solution structures are complete in the sense that they converge to the exact solution of a problem. Proposed method is unique amongst meshfree methods because RFM solution structures can be constructed to satisfy all boundary conditions exactly. This property allows that the solution procedure does not require mapping the geometry of the domain. Combination of atomic basis functions and solution structure method gives numerical solutions that have characteristics of analytical solutions because solution structure method enables exact treatment of boundary conditions while ABFs ensure numerical solutions with desired level of accuracy.

RFM has one significant disadvantage when compared to the usual spatial discretization and other meshfree methods: the solution structure is constructed, differentiated, and integrated at runtime. Thus, a fully automatic implementation of RFM requires automatic construction of implicit functions and structures, given a geometric model and prescribed boundary conditions, automatic differentiation of the constructed functions and structures at various points of the domain, constructing and solving the resulting linear system for values of the coefficients [C.sub.i] in the undetermined functional component [PHI]. In our examples, [PHI] is a uniform set of atomic Fup basis functions, so the resulting sparse banded system can be solved in linear time using standard techniques. All these tasks are feasible today with existing technology, so the proposed method is promising numerical method for solving technical problems.

Our future research will be directed towards the development of numerical models for analysis 2D heat transfer problems, groundwater flow modelling in karst aquifers, nonlinear problems, plate bending problems, and towards automation of the numerical process with development of adaptive techniques.

DOI: 10.2507/27th.daaam.proceedings.036

7. Acknowledgments

This paper is a result of research in the frame of national project "Groundwater flow modelling in karst aquifers" funded by croatian Science Foundation.

8. References

[1] Atluri, S. N. (2005). Methods of Computer Modeling in Engineering & the Sciences, Volume I, Tech Science Press, ISBN 0-9657001-9-4, University of California, Irvine

[2] Griebel, M. & Schweitzer, M. A. (Eds.) (2003). Meshfree Methods for Partial Differential Equations, Vol. 26, Springer-Verlag, ISBN 3-540-43891-2, Berlin

[3] Liu, G. R. (2003). Mesh Free Methods: Moving Beyond the Finite Element Method, CRC Press LLC, ISBN 0-8493-1238-8, Boca Raton

[4] Belytschko, T.; Krongauz, Y.; Organ, D.; Fleming, M. & Krysl, P. (1996). Meshless methods: An overview and recent developments, Comput. Methods Appl. Mech. Eng., Vol. 139, No. 1-2, 1996, pp. 3-47

[5] Fries, T. P. & Matthies, H. G. (2004). Classification and Overview of Meshfree Methods, Institute of Scientic Computing Technical university Braunschweig, Brunswick

[6] Chen, W.; Fu, Z. J. & Chen, C. S. (2014). Recent Advances in Radial Basis Function Collocation Methods, Springer Briefs in Applied Sciences and Technology, ISBN 978-3-642-39571-0, Springer

[7] Gunter, F. C. & Liu, W. K. (1998). Implementation of boundary conditions for meshless methods, Comput. Methods Appl. Mech. Eng., Vol. 163, 1998, pp. 205-230

[8] Kantorovich, L. V. & Krylov, V. I. (1958). Approximate Methods of Higher Analysis, Interscience Publishers, New York

[9] Rvachev, V. L. (1982). Theory of R-functions and Some Applications, Naukova Dumka, Kiev. In Russian

[10] Rvachev, V. L.; Sheiko, T. I.; Shapiro, V. & Tsukanov, I. (2000). On completeness of RFM solution structures, Comput Mech, Vol. 25, 2000, pp. 305-317

[11] Rvachev, V. L. & Sinekop, N. S. (1990). R-functions Method in Problems of the Elasticity and Plasticity Theory, Naukova dumka, Kiev. In Russian

[12] Rvachev, V. L. & Sheiko, T. I. (1995). R-functions in boundary value problems in mechanics, Applied Mechanics Reviews, Vol. 48, No. 4, 1995, pp. 151-188

[13] Rvachev, V. L. & Rvachev, V. A. (1971). On a finite function, Dokl. Akad. Nauk Ukrainian SSR, ser. A, 1971, pp. 705

[14] Rvachev, V. L. & Rvachev, V. A. (1979). Non-classical Methods for Approximate Solution of Boundary-Value Problems, Naukova dumka, Kiev. In Russian

[15] Kravchenko, V. F. (2003). Lectures on the Theory of Atomic Functions and Their Some Applications, Radiotechnika, Moscow

[16] Kozulic, V.; Gotovac, B. & Colak, I. (2006). Multilevel mesh free method for the torsion problem, In: DAAAM International Scientific Book 2006, Katalinic, B., (Ed.), Chapter 29, pp. 365- 384, DAAAM International Vienna, ISBN 3-901509-47-X, Vienna

[17] Kozulic, V.; Gotovac, B. & Sesartic, R. (2008). Mesh free modeling of the curved beam structures, In: DAAAM International Scientific Book 2008, Katalinic, B., (Ed.), Chapter 34, pp. 395- 408, DAAAM International Vienna, ISBN 978-3-901509-66-7, Vienna

[18] Kozulic, V. & Gotovac, B. (2015). Computational modeling of structural problems using atomic basis functions, In: Advanced Structured Materials, Vol. 70: Mechanical and Materials Engineering of Modern Structure and Component Design, Ochsner, A.; Altenbach, H., (Eds.), Chapter 17, pp. 207-230, Springer

[19] Gotovac, B. & Kozulic, V. (1999). On a selection of basis functions in numerical analyses of engineering problems, Int. J. Eng. Modelling, Vol. 12, No. 1-4, 1999, pp. 25-41

[20] Prenter, P. M. (1989). Splines and Variational Methods, John Wiley & Sons, ISBN 0-471-50402-5, New York

[21] Timoshenko, S. P. & Goodier, J. N. (1961). Theory of Elasticity, McGraw- Hill, New York

[22] Lurie, A. I. (1970). Theory of Elasticity, Nauka, Moskva

[23] Kostrencic, Z. (1982). Theory of Elasticity, School Book, Zagreb

[24] Pilkey, W. D. (2005). Formulas for Stress, Strain, and Structural Matrices, John Wiley & Sons

This Publication has to be referred as: Kozulic, V[edrana] & Gotovac, B[laz] (2016). Numerical Solution of Poisson's Equation in an Arbitrary Domain by Using Meshless R-Function Method, Proceedings of the 27th DAAAM International Symposium, pp.0245-0254, B. Katalinic (Ed.), Published by DAAAM International, ISBN 978-3-902734-08-2, ISSN 1726-9679, Vienna, Austria

Caption: Fig. 1. Basic R-operations

Caption: Fig. 2. Two-dimensional domain of arbitrary shape with mixed boundary conditions

Caption: Fig. 3. Layout of the Fup basis functions and positions of collocation points

Caption: Fig. 4. Layout of Fup4 basis functions in relation to the boundaries of the one-dimensional domain

Caption: Fig. 5. a) Triangular cross section; b) Uniform grid of points that covers the domain

Caption: Fig. 6. Stresses [[tau].sub.xz] (a) and [[tau].sub.yz] (b) computed by solution structure method

Caption: Fig. 7. a) Square cross section; b) Uniform grid with the minimum number of collocation points; c) Stress function computed by 55 x 55 grid

Caption: Fig. 8. Convergence diagram for: a) a value of the stress function in the point A; b) value of the maximum shear stress

Caption: Fig. 9. Contour and the shape of the warping function w obtained by RFM

Caption: Fig. 10. Convergence diagram for the value of the maximum shear stress

Caption: Fig. 11. a) Geometry of the domain; b) Contour of the function [omega]

Caption: Fig. 12. Computed shearing stresses [[tau].sub.yz] along the horizontal axis of symmetry