Exact Boundary Derivative Formulation for Numerical Conformal Mapping Method

Conformal mapping is a useful technique for handling irregular geometries when applying the finite difference method to solve partial differential equations. When the mapping is from a hyperrectangular region onto a rectangular region, a specific lengthto-width ratio of the rectangular region that fitted the Cauchy-Riemann equations must be satisfied. In this research, a numerical integral method is proposed to find the specific length-to-width ratio. It is conventional to employ the boundary integral method (BIEM) to perform the conformal mapping. However, due to the singularity produced by the BIEM in seeking the derivatives on the boundaries, the transformation Jacobian determinants on the boundaries have to be evaluated at inner points instead of directly on the boundaries. This approximation is a source of numerical error. In this study, the transformed rectangular property and the Cauchy-Riemann equations are successfully applied to derive reduced formulations of the derivatives on the boundaries for the BIEM. With these boundary derivative formulations, the Jacobian determinants can be evaluated directly on the boundaries. Furthermore, the results obtained are more accurate than those of the earlier mapping method.


Introduction
The finite difference method (FDM) is a conventional numerical method commonly used in computational science because partial differential equations can be directly discretized [1][2][3]. However, when the computational region is irregular, boundary values must be evaluated through interpolation, extrapolation, or both, which results in greater computation costs and decreases accuracy [4,5].
A common solution to circumvent this problem is to transform the irregular region into a rectangular region and solve partial differential equations on the rectangular region [6][7][8][9][10]. According to a study by Thompson et al. [11], a system of elliptic partial differential equations can guarantee a one-to-one transformation between the physical and computational domains. When the transformation is conformal, the number of additional terms introduced by the transformation of the governing equation is the minimum. Thus, inaccuracy can be avoided.
By applying the conformal mapping method to build an effective grid-generation system, Tsay et al. [12] employed the boundary integral equation method (BIEM) (also called the boundary element method, BEM) to solve the Laplace equations. The BIEM is preferable to the Laplace-equationgoverned transformation process for a number of reasons: (1) it is not necessary to discretize the entire domain; (2) solutions can be evaluated exactly; (3) derivatives of variables can be evaluated directly without difference schemes; and (4) once the solutions at all the discretized boundary nodes are obtained, at any given point in the domain the solution as well as the partial derivatives can be easily found without remeshing. When the conformal mapping is performed, the governing equation of the physical domain will be transformed into a formulation in terms of derivatives with respect to the coordinates. Although the derivatives in the domain can be directly evaluated by using the BIEM, it suffers from the fact that the derivatives cannot be evaluated exactly on the boundaries [13,14]. This drawback makes one be only able to evaluate the derivatives close to the boundaries as the boundary derivatives.
A limitation of the Laplace-equation-governed transformation is that the region to be transformed should be a 2 Mathematical Problems in Engineering quadrilateral with right-angled corners, which is known as a hyperrectangular region [12]. When a hyperrectangular region is transformed into a rectangular region by using conformal mapping, the length-to-width ratio of the rectangular region is a particular value [15]. To find this ratio, Seidl and Klose [15] converted the Cauchy-Riemann equations into a set of two Laplace equations and proposed an iterative procedure with a finite difference solution. In this research, an explicit integral method was proposed to find this ratio, without the time consumption of solving a system of linear equations.
Although the applicability of the conformal transformation used by Tsay et al. [12] to arbitrary irregular regions could be extended by using a sequence of transformations [16], the applicability of the grid-generation system proposed by Tsay et al. [12,14,16] is still limited due to the abovementioned problem of inaccurate boundary partial derivatives. This problem not only precludes the derivatives and the transformation Jacobian determinants from being evaluated exactly on the boundaries, but also decreases the accuracy of the finite difference computation after the transformation.
In this paper, by applying the rectangular properties of the transformed region and the Cauchy-Riemann equations, the derivatives and the Jacobian determinants of the transformation can be evaluated on the boundaries. This evaluation significantly improves the transformation accuracy.

Mathematical Description of the Coordinate Transformation
Based on the Cauchy-Riemann equations, the forward conformal transformation, from the physical domain ( -coordinates) to the computational domain ( -coordinates), can be reduced to a problem governed by two Laplace equations: with the boundary conditions where n is the outward normal vector of the boundaries and ∇ is the gradient operator. Here, the range of either or can be chosen arbitrarily while the range of the other coordinate is restricted by the satisfaction of the Cauchy-Riemann equations: That is, when satisfying the Cauchy-Riemann equations, there is always a specific value of 0 that corresponds to a given 0 , and vice versa. To find this specific value, Seidl and Klose [15] proposed an iterative scheme. Besides, Wu et al. [17] showed that the transformation is not conformal but still orthogonal if one simply chooses both values of 0 and 0 arbitrarily. In Section 4, we will propose a numerical integral method to find the specific ratio of 0 to 0 which satisfies the Cauchy-Riemann equations. By applying the divergence theorem, (1a) and (1b) can be represented as an integral equation [13]: where Φ represents or ; Γ is the bounding contour of the problem domain ; is the out-normal direction of Γ; represents a source point; represents a boundary point; is the distance from to ; = 2 when is located inside the domain; and is /2 or when is located on the corners or the straight boundaries, respectively, as shown in Figure 1. An elementwise local coordinate system ( -) is adopted to discretize the geometry as linear elements, as shown in Figure 2, where the subscripts and + 1 represent the starting and ending nodes of an element and the subscript represents the th source point. After discretizing the boundary and choosing the source points as the nodal points of the elements, the unknown function values, Φ, and their normal derivatives, Φ/ , can be formulated as simultaneous equation system and then be solved [13,16]. To ease the explanation, the above solution for the unknowns of boundary elements is named boundary-unknown solving step in this paper.
The inverse conformal transformation is also governed by Laplace equations, whose function values are and with independent variables of and : and the following Cauchy-Riemann equations should be satisfied: = . (6b) Since the values of and on the boundary are already known, the boundary condition of the inverse transformation is the Dirichlet type. The formulation of integral equations for (5a) and (5b) is the same as (4), but Φ represents or . A sketch of the transformation is illustrated in Figure 3.
After Φ and Φ/ on the boundary are solved, Φ inside the domain can be obtained by using (4) explicitly, where is now located either inside the domain or on the boundary. In this research, a hyperrectangular region is transformed into a rectangular region. It is easy to generate Mathematical Problems in Engineering     an orthogonal grid in the rectangular region by using (4), where Φ represents the coordinate of the rectangular region. Then, each grid node is mapped onto the hyperrectangular region, by using (4) again, where Φ represents the coordinate of the hyperrectangular region now. The above evaluation of grid nodes is named grid-generation step in this paper. The regular orthogonal grid in the rectangular region is very convenient for applying FDM. To perform FDM on the grid in the rectangular region, the governing equation usually has to transform into a formulation with Jacobian determinants [18][19][20]. For example, the Biharmonic equation, ∇ 4 ( , ) = 0, has to transform into ∇ 2 ( ( , ) ∇ 2 ( , ) ) = 0, where and are coordinates of the hyperrectangular region; and are coordinates of the rectangular region; ( , ) is the Jacobian determinant for the forward transformation.
The Jacobian determinants for the forward transformation and inverse transformation are defined, respectively, as below: With the help of Chain rule and these Jacobian determinants, the relationships between the partial derivatives of -and -can be built as follows: = 1 ( , ) .

The Derivatives and Jacobian Determinants on the Boundaries
The Jacobian determinant is a combination of partial derivatives. Since these derivatives are evaluated by the BIEM, singularities occur when the unknown derivatives are located on the boundaries. A previous research has suggested that these boundary derivatives can be approximated by evaluating the inner points that are very close to the boundaries [16]. In this study, the derivative formulations of the BIEM are carefully investigated. By applying the Cauchy-Riemann equations, a method of evaluating the boundary derivatives is proposed. Unlike conventional methods that just calculate Mathematical Problems in Engineering 5 the derivative close to the boundary as succedaneums, this method can produce accurate numerical results exactly on the boundaries. The partial derivatives of Φ can be found by differentiating (4), and then the following formulations are obtained [13]: The function value and its normal derivative on the boundary, Φ( ) and ( / )Φ( ), are already known after the evaluation of (4). 1 and 2 represent and , respectively, when the function value Φ is or . Conversely, 1 and 2 represent and , respectively, when the function value Φ is or . After the geometry discretization, the following derivative terms can be represented with the local coordinate system ( -): By applying (11a)-(11e), (10a) and (10b) can be represented with the following formulations, respectively: Liggett and Liu [13] have indicated that these cannot be evaluated on the boundaries because singularities will occur. However, they are actually applicable on the boundaries in the conformal mapping process, through the following investigation.
The discretization of these boundary integral equations is performed by linear approximations, and the potential and its normal derivative are subsequently approximated by for ≤ ≤ +1 , where the subscripts and + 1 represent the starting and ending nodes of an element. After substituting (13) into (12a) and integrating each term in the local coordinate system ( -), the discretized formulation of the integral equation for the derivative in the horizontal direction can then be represented as where represents the th discretized element on the boundary, is the number of the elements, and in which During the grid-generation step, when is located on the top or bottom boundary of the transformed rectangular region, where = 0, singularity will occur in (16b) and (16f). However, at the same time, the value of is either 0 or which makes sin = 0. This zero value of sin makes the singularities of (16b) and (16f) vanish. Besides, it should be emphasized that the value of the term (tan −1 ( +1 / ) − tan −1 ( / )) is also zero when = 0 [14]. After substituting = 0, sin = 0, and (tan −1 ( +1 / ) − tan −1 ( / )) = 0 into (16a)-(16j), the original discretized equations are reduced to the following formulations: (17c) This result not only is relatively simple, but can also be determined, except at the source point ( = 0 or +1 = 0). Similarly, (12b), the derivative formulation in the vertical direction, can be discretized to the following equation: Notice that the formulation of 2 11 , 2 12 , 2 11 , 2 12 , can be obtained by changing sin to cos and cos to sin in (16a)-(16j), and similar reduced formulations can be derived when is located on the left or right boundary of the transformed rectangular region: The Jacobian determinant of the inverse transformation, from the rectangular region to the irregular region, defined as (7b), is a combination of partial derivatives in the rectangular region. On the top and bottom boundaries of the rectangular region, the horizontal partial derivatives / and / can be evaluated using (14)-(17c). Instead of using the BIEM directly, high-accuracy vertical derivatives / and / can be evaluated by applying the Cauchy-Riemann equations for inverse transformation, (6a) and (6b). Similarly, on the right and left boundaries, the vertical derivatives / and / can be derived using (18)-(20c), and the horizontal partial derivatives / and / can be subsequently Figure 5: The integral paths of the numerical integral method for finding the length-to-width ratio to apply conformal mapping from a hyperrectangular region onto a rectangular region.
Present approach path Approach path by Tsay and Hsu [16] = sin(−0.5)cosh(y) − cos(−0.5)sinh(y) = − x (−0.5,y) | | Figure 6: The unit square with complicated boundary conditions. In order to evaluate the derivative on a source point, this research approached the point along the boundary, while Tsay and Hsu [16] approached the point from inside the boundary. evaluated by applying the Cauchy-Riemann equations. An illustration of this solution is provided in Figure 4. Once these partial derivatives on the boundaries are evaluated using the reduced formulations, high-accuracy Jacobian determinants of the inverse transformation, ( , ) , can be constructed directly on the boundaries using (7b). After constructing ( , ) , the derivatives of the forward transformation, / , / , / , and / , can be obtained using (8e)-(8h). Finally, the Jacobian determinant of the forward transformation, ( , ) , can be constructed by using (7a).
Taking advantage of the transformed rectangular geometry, most singular terms in the conventional derivative formulations vanish, and the derivatives can then be evaluated on the boundaries. Furthermore, most terms in the reduced formulations are equal to zero, which implies far fewer computational errors than using conventional formulations. Although the present reduced formulations still cannot evaluate the derivatives at the source points, adopting approximations with the present formulations can yield accurate results. An illustrated example is provided in Section 5.

Length-to-Width Ratio of Conformal Modulus
To perform conformal mapping from a hyperrectangular region ( -) to a rectangular region ( -), a specific ratio of length to width have to be fit on the rectangular region. An iterative method has been presented by Seidl and Klose [15].
In this section, we propose a numerical integral method to find this ratio as follows.
Let the width of the rectangular region ( 0 ) be equal to a suitable value. Then / , / , and ( , ) can be obtained  by solving ∇ 2 = 0 and (10a) and (10b). The corresponding length ( 0 ) satisfying the Cauchy-Riemann equations can be evaluated by the following integral equation: where / = / and / = − / . By using difference method, (21) was discretized to the following equation: where the subscript represents the th node on the integral path; the subscript * means a middle location between and + 1; is the number of the nodes. The numerical integration is performed from the side mapped onto = 0 to the other side mapped onto = 0 . Theoretically, the integration is independent of integral path. Since the numerical error is unavoidable, several integral paths are chosen in the hyperrectangular region ( -) to reduce the numerical error, for example, Figure 5. Then, perform integration using (22), respectively. After the integrations are completed, the average of 0 s from these integral paths, 0 , is used as the length of the rectangular region that fitted Cauchy-Riemann equations.

Results and Discussion
This research proposes a numerical integral method to calculate the length-to-with ratio of a rectangular region when a conformal mapping is performed from a hyperrectangular region to the rectangular region. Besides, an improved scheme to directly solve for the partial derivatives on the boundaries when conformal mapping is performed using the BIEM is also proposed in this research. By applying the properties of rectangular geometry, the reduced formulations of BIEM discretization, (17a)-(17c) and (20a)-(20c), significantly improved the accuracy of the calculation results. In order to validate the improved boundary-derivativesolving scheme, three examples are provided in this section: the first is a transformation from a unit-square region into another unit-square region; the second is a transformation from an arc region into a rectangular region; the third is a transformation from a wave-block region (one side is a wave curve and other sides are straight lines) into a rectangular region.

Unit-Square Transformation.
The case of unit-square transformation and the error estimates at a source point are introduced first. When the calculation point is located at a source point, neither the conventional nor the present numerical scheme can provide the derivative. However, since the governing equation is the Laplace equation, the analytical solution should be continuous and smooth. These properties make it reasonable to approach the derivative at the source point using a numerical result from a nearby point. Tsay and Hsu [16] suggested that an approach distance of 10 −6 from the inner domain be used to perform the derivative calculation, while in this research the approximation approach along the boundaries is used. To assess which approach can provide results that are more accurate, a unit-square region with two complicated Dirichlet and two Neumann boundary conditions is considered. The Dirichlet boundaries were set on the top and bottom and the Neumann boundaries were In order to capture the complicated analytical solution, 100 elements were applied on each side of the computational region, and then a series of computations in the derivative with respect to near the source point (0.0, 0.5) were performed. Figure 7 provides the percentage errors between the analytical solutions at the source point (0.0, 0.5) and the numerical results near the source point. Figure 7 shows that the percentage errors from the approximations made by Tsay and Hsu and in the present study both decrease as the calculation point approaches the source point. However, the present study generates significantly fewer errors. In addition, when the approach distance is shorter than 10 −12 , Tsay and Hsu's numerical result becomes singular, whereas the present scheme still provides an accurate result.
Unfortunately, a numerical result cannot be achieved at the corners using this boundary approach. The reason for this is as follows. Either (17a)-(17c) or (20a)-(20c) can be applied to overcome the singularity that occurs when the local coordinate is equal to 0 for an element on one side close to a corner. However, the local coordinate of the element on the adjacent side is very small, while the original formulations in (16a)-(16j) are applied. This substitution of a small in the original formulations results in a singularity.
For example, in order to calculate the derivative value on corner C in Figure 8, the calculation should be performed a small distance away, along the boundary, at C because the corner is a source point. When the derivative calculation is performed using (18) and (19a)-(19d), the reduced formulations in (20a)-(20c) are applied to solve the singularity on the boundary elements on the left-hand side, . Original formulations in (16a)-(16j) are applied on the boundary elements at the bottom, . The small distance, of , results in another numerical singularity. To overcome this problem, the conventional forward and backward difference schemes are suggested, although they cannot provide the same accuracy as the present scheme does at other points.

Arc-Region Transformation.
The second example is a transformation from an arc region into a rectangular region, as shown in Figures 9(a) and 9(b), where is from /4 to 3 /4; out = 2; in = 1; and 0 = 1. The remaining unknown 0 can be determined by the numerical integral method proposed in this research, using 5 integral paths with path = 1.3, 1.4, 1.5, 1.6, and 1.7, respectively, as the dash lines shown in Figure 9(a). The development of 0 with the number of discretized elements of the integral path, , is shown in Figure 10(a). As is increasing, 0 approaches the analytical solution in (24a)-(24h). The standard deviation of these 0 s from different integral paths is observed to be independent of , when is greater than 200, as shown in Figure 10(b). The error evaluations shown in Figure 11 illustrate that the percentage errors decrease with . When is equal to 500, the percentage error is only 1.55 × 10 −5 (%). Since the analytical solution is available in this case, conformal mapping was performed using the analytical 0 = / ln 4 = 2.26618, which derived from (24g). To perform the conformal mapping with BIEM, either inner or outer curve boundary of the arc region is discretized into 200 elements, and each straight boundary is discretized into 100 elements. By using the governing equations and boundary conditions for the forward transformations shown in Figures 12(a) and 12(b), AB, BC, CD, and DA sections in the arc region are mapped onto AB , BC , CD , and DA sections in the rectangular region, respectively. Then, the inverse transformation is performed with the governing equations, ∇ 2 = 0 and ∇ 2 = 0, and all the Dirichlet type boundary conditions. After the conformal mapping is completed, a 13 × 5 regular grid is built in the rectangular region and then mapped onto the arc region, as shown in Figures 13(a) and 13(b).
The numerical boundary derivatives and Jacobian determinants were examined with their analytical values. The analytical solutions of the forward transformation and the related derivatives are listed below: = ln √ 2 + 2 + , (24b) where Then, the analytical solutions of the inverse transformation and related derivatives are as follows: Finally, the analytical Jacobian determinants can be derived using above equations and shown as below: which agrees with (9). Figures 14(a)-14(d), respectively, compare the errors of / , / , / , and / , on the boundaries evaluated using the scheme suggested by Tsay and Hsu [16] and the scheme provided in the present research.  Figure 14: Comparison of the percentage errors of (a) / , (b) / , (c) / , and (d) / between the present results and those from Tsay and Hsu [16]. After adopting the present numerical scheme, the derivatives on the boundary are much more accurate. along the boundary, excluding the corners. The vertical axis represents the relative errors, , defined as where represents analytical solutions and represents numerical solutions. Although from Tsay and Hsu's scheme is small, it is observed to be further smaller than in the present scheme. The errors of ( , ) and ( , ) on the boundaries by using the two methods are compared in Figures 15(a) and 15(b). Since ( , ) has a direct relation with ( , ) , as shown in (9) and (26), the trends of error distribution of ( , ) and ( , ) are similar. As can be seen, the present numerical scheme provides more accurate results. To quantify the improved accuracy, an average error was defined as = ∑ / , where is the percentage error (%) at each boundary node and is the number of the boundary nodes excluding the corners. Although the scheme of Tsay and Hsu [16] Figure 15: Comparison of the percentage errors of (a) ( , ) and (b) ( , ) between the present results and those from Tsay and Hsu [16]. After adopting the present numerical scheme, the transformed Jacobian determinants on the boundary are much more accurate. results with = 0.074% in both ( , ) and ( , ) , the average error farther sharply decreases to = 0.00068%, after using the present scheme.

Wave-Block Region Transformation.
The third example is a transformation from a wave-block region into a rectangular region, as shown in Figures 16(a) and 16(b), where = 1 and elements is equal to 25600, 0 was convergent to 4.966196, used in the conformal mapping process. was observed to decrease with the number of the discretized elements. Orthogonality of the grid is a criterion to examine whether 0 approached the analytical 0 . The following relative errors quantify the grid orthogonality [17]: ) .   The smaller the values of AREo and MREo are, the higher the orthogonality of the grid will be. A linear relation is observed between not only AREo and 0 , but also MREo and 0 , as shown in Figure 18. When 0 is equal to 4.966196, AREo and MREo are equal to 8.14 × 10 −6 and 3.23 × 10 −4 , respectively.
The small values of AREo and MREo indicate that 0 is very close to the analytical value.
To perform the conformal mapping with BIEM, the boundary is discretized into 1400 elements, where AB and CD sections are discretized into 300 elements and AD and BC sections are discretized into 400 elements. AB, BC, CD, and DA sections of the wave-block region are mapped onto A B , B C , C D , and D A sections of the rectangular region, respectively. In the forward transformation, the governing equations and boundary conditions are set as in Figure 19. In the backward transformation, the governing equations are ∇ 2 = 0 and ∇ 2 = 0 with all the Dirichlet type boundary conditions. After the conformal mapping was completed, a 50 × 30 regular grid is built in the rectangular region and then mapped onto the wave-block region, as shown in Figures  20(a) and 20(b).
The boundary derivatives of the conformal mapping are evaluated using the improved scheme proposed in this research. Since there is no analytical solution in this case, a conventional numerical method, FDM, is used to examine  the correctness of the proposed scheme. It is applied in the rectangular region ( -coordinates) to obtain the boundary derivatives / , / , / , and / . Then, (8e)-(8h) are used to obtain the boundary derivatives / , / , / , and / . To increase the accuracy, the second-order formulations of FDM are used, as shown in the following: where Δ = 0.01 in this case; the subscript represents the th point to evaluate its derivative; Φ ± Δ represents the function value ( or ) from the th point plus or minus Δ distance. Since the boundary-unknown solving step has been finished, any function value can be obtained by using (4)  On the other hand, the absolute values of numerical errors of FDM are observed to be about 5 × 10 −4 on AB and CD sections, much higher than the errors of the present scheme. Similar results take place on / and / , which are not shown in this paper for elegance. Despite the tiny numerical errors, the boundary derivatives are observed to be consistent by using the proposed scheme and FDM. These consistent results demonstrate that the proposed scheme can be applied to irregular regions that have no analytical solutions for the conformal mapping. Since the improved boundary derivative formulations proposed in this research are coupled with the Cauchy-Riemann equations, their accuracy are affected by the correctness of the length-to-width ratio of the rectangular region. The consistent results of the boundary derivatives validate not only the proposed boundary derivative formulations but also the integral method for the length-to-width ratio.

Conclusions
A numerical integral method to find the length-to-width ratio of the rectangular region that fitted Cauchy-Riemann equations is proposed in this research. An arc example demonstrates that the numerical results can approach the analytical solution with negligible error (1.55 × 10 −5 %). A wave-block example shows that by using the length-towidth ratio produced by the proposed integral method the conformal mapping can perform successfully.
By applying the geometric property of the transformed rectangular region and the Cauchy-Riemann equations, the derivatives and Jacobian determinant of the numerical conformal mapping method developed by Tsay and Hsu [16] can be evaluated on the boundaries, and more accurate results are obtained. The approach for calculating the Jacobian determinant at the source points adopted in the present numerical scheme can produce numerical results that are more accurate than the conventional approach. Although the derivatives and the Jacobian determinant at the four corners of the transformed rectangular region cannot be improved by the present scheme, this research was able to improve the accuracy of the numerical conformal mapping on the boundaries and makes the mapping method more reliable when solving problems that are sensitive to accuracy.