Abstract
This paper proposes a simple surface interpolation attaining tangent-plane continuity. It is a natural extension of the local quadratic interpolator developed by the author (2005) in one of his works, which has already been applied successfully to diverse engineering problems. The methodology presented in this paper inherits most of the advantages possessed by the scheme. That is, (i) The algorithm is efficient and completely local requiring only the position vectors and normals given at the nodes of a patch, and hence it is suitable for parallel processing. (ii) It converges rapidly to the given surface with the increase in the number of nodes. (iii) Singular points (apexes, sharp edges, etc.) and nonmanifolds can be treated quite easily. (iv) Because of the minimization criteria assigned to the surface coefficients, it is rather robust and amenable to computational analyses. Validity and effectiveness of the proposed technique are demonstrated through numerical examples.
1. Introduction
There is a significant gap between the requirements on geometric models in the CAD and computational science communities. In the former, surface descriptions using few patches with large degree of freedom and high level of continuity (e.g., NURBS, Bรฉzier, Gregory) are considered to be desirable. On the other hand, models for numerical simulation (e.g., the finite element method) are almost always represented by fine meshes, where the major interests are in fast, stable, and accurate recovery of the surface information lost during the process of discretization. With this as background, the author proposed an interpolation scheme suitable for such analyses [1]. It has already been applied successfully to engineering problems including(A)high-precision machining [2],(B)simulation of elastoplastic mechanics [3, 4],(C)ray tracing of optical devices [5].
All those applications prohibited the usage of traditional sophisticated surface descriptions, due to severe tolerance as well as geometrical and physical complexity of the systems. In (A) the method readily yielded positional accuracy of around โm, which could not be attained by other interpolators showing poor or no convergence. Less than 10% of the linear surface patches turned out to be necessary to obtain the sufficient precision in (B). The technique reduced the number of the patches required for the ray tracing in (C) by a factor of . The interpolator is now referred to as Nagata patch by its users, and hence this paper follows this convention, though the name of the author is frequently dropped for simplicity. The characteristics of the patch listed below have enabled the aforementioned significant performance improvements. (i)The algorithm is completely local, requiring only the position vectors and normals given at the vertices of each patch. Therefore, it is suitable for parallel processing. (ii)The algorithm is simple, computationally inexpensive, and hence amenable to various geometrical and physical evaluations. (iii)Because the formulation accounts for discontinuity (multiplicity) of normals, sharp edges and singular points as well as non-manifolds can be treated quite easily. (iv)It has the minimum degree (two) of interpolation necessary for representation of the curvature. This property is desirable especially for ray tracing, contact problems, and so forth, which involve implicitization and inversion, since closed-form solutions may be obtained. (v)The interpolation assures the continuity, and converges to the original surface rapidly with the increase in the number of nodes, even in the presence of the singular features. Hence error in the normals can be sufficiently small using rather few patches. This implies asymptotic smoothness.
Other differences of the method from traditional approaches (smooth local interpolators in particular) are discussed extensively in [1], and hence are not repeated here.
Another merit of the patch is that it can easily be extended to establish tangent-plane continuity, as reported in this paper. The proposed interpolation is named Nagata patch here for convenience of reference.(The term does not imply that the interpolation always generates surfaces; it can reproduce singular points (apexes, sharp edges, etc.) or may yield cusps similar to the patch.) The new algorithm involves correction using a simple rational function retaining the boundary of the patch. Hence it is no longer quadratic, but it inherits all other desirable features from the patch, such as complete locality and capability of handling multiple normals. Similar to the patch, the present algorithm assigns minimization criteria (inherent in the generalized inverse) to the surface coefficients, and hence is rather robust. This is notable since accuracy and stability are not compatible with smoothness for most surface interpolators presently available.
This paper is organized as follows. The next section formulates the problem, briefly reviews the Nagata patch, identifies the conditions on its correction for the tangent-plane continuity, and gives a simple solution satisfying the requirement, thus deriving the interpolation scheme in Section 2.2.5. Section 3 assesses the convergence and accuracy of the interpolation through numerical examples, followed by Section 4 summarizing the results with a scope for future work.
2. Formulation
Assume that position vectors and normals are given at the vertices of a triangular mesh with general topology. The objective here is to recover the curvature and smoothness of the surface through interpolating each patch independently. The process consists of the following three steps:(A)to replace each edge of a patch with a curve orthogonal to the normals given at the endpoints of the edge,(B)to determine a parametric quadratic polynomial patch ( Nagata patch) reproducing the modified boundary to keep the continuity,(C)to correct the patch preserving the boundary in order to obtain tangent-plane continuity.
Algorithms for (A) and (B) were already proposed by the author in [1]. However, they are essential for the present approach and hence are summarized in the next section for completeness. The patch is then upgraded to enable local smooth interpolation.
2.1. The Patch
For the time being, assume a smoothsurface. Consider a curve segment on the surface as illustrated in Figure 1. Its endpoints and have the position vectors , and the unit normals , (In this paper, a bold typeface is used for vectors (column matrices), and bold symbols in square brackets denote matrices.). They are supposed to be known as input.
The step (A) mentioned at the beginning of Section 2 approximates the curve segment using the quadratic function to reproduce the position and normal vectors at the endpoints , . Here denotes the parameter for the curve, and is the vector connecting the endpoints. The coefficient in (2.1) introduces curvature to the segment, and hence is named the curvature parameter. It is determined from the boundary conditions at , as Here is the zero vector. For general surfaces which may have multiple normals at , , the curvature parameter can be extended as where , are constructed by (32) of [1] from the position and normal vectors at , . The generalized inverse of the former is explicitly written as where is the zero matrix, and , are coefficients defined by () and () of [1], respectively. Note that , , are dimensionless and should be treated as zeros when their absolute values are below a small positive threshold , for avoiding numerical instability.
The interpolator (2.1) is at most quadratic, and hence its approximation capability is rather limited. If the magnitude of the curvature parameter is too large, failure in the aforementioned singularity detection or insufficient resolution of the original data is suspected. This situation can be handled as follows. Note that the distance between and is (Figure 1) and the curvature parameter also has the dimension of length. If the latter increases, the interpolated curve deviates from the line segment . Therefore, can be used as a criterion for inadequate interpolation. Here is a dimensionless parameter and hence can be a constant regardless of the geometric scale. If the above holds, the next expression on the right-hand side of (2.5) is chosen for evaluating the generalized inverse irrespective of the decision based on the threshold , and the curvature parameter of (2.4) is recalculated. If it still satisfies (2.6), the same process is repeated. In the worst case, the generalized inverse (2.4) becomes the zero matrix resulting in making the curve segment linear. Since the case is obviously against (2.6), this iteration always terminates and is safe. Similar treatment should be taken also for (2.3).
With the above algorithm as a basic tool, the triangular patch in Figure 2(a) is interpolated. Its geometry can be described by where denotes the position vector to a point on the patch, and , are its parameters defined within the domain of Figure 2(b).
(a)
(b)
In Figure 2(a), the position vectors to the vertices , , : and the unit normals at the three points are assumed to be known. The patch approximates the surface using the quadratic polynomial The edges of the patch are first interpolated in the form of (2.1) as whose coefficients are given by Here the curvature parameter is evaluated by (2.3) for smooth surfaces or by (2.4) for general cases. The parametric representation (2.10) of the patch is uniquely obtained from the boundary curves of (2.11) as In the subsequent discussion, symmetric expressions play an essential role. For this purpose, barycentric coordinates are introduced to replace the surface parameters , . In the parameter space of Figure 2(b), the vertices of the triangular patch have the following position vectors: Hence the parameter vector for a point on the patch can be described as where , , are the barycentric coordinates for . They are in one-to-one correspondence with the surface parameters provided that Equations (2.15a)-(2.15b) can be solved for the barycentric coordinates as Equations (2.16) and (2.7) show that the following holds on the patch: Let be the vertex for the parameter vector and its opposite edge. Equations (2.15a)-(2.15b) give their representation in the barycentric coordinates asHere denotes the set of all cyclic permutations of the indices 0, 1, and 2; that is, Introducing the barycentric coordinate vector equations (2.18a)-(2.18b) can be rewritten aswhere are unit vectors defined by As (2.21b) implies, the direction along the edge in the space of barycentric coordinates is given by Application of (2.15b) and (2.16) to (2.13) yields the position vector in the physical space as the following homogeneous polynomial: Hereare constant vectors. Equation (2.24) is equivalent to the quadratic Bernstein form where is a multi-index, whose valid values form the set , and denote the Bernstein basis functions. Substituting the above into (2.26) and comparing the result with (2.24) yield the coefficients which are the quadratic Bรฉzier control points for the patch.
Symmetry of (2.24) greatly simplifies the formulation. For instance, a position vector on any edge can be expressed through (2.18b) as
2.2. The Patch
Despite its extreme simplicity, the aforementioned interpolation has high geometric accuracy. Hence it has already found many engineering applications as mentioned in Section 1. However, it generally gives discontinuous normals across patches, even if the original surface is smooth. (In this study, parallel normals are considered to be identical.) This may not be desirable when visual quality needs to be respected or the tangent-plane continuity is essential for the analysis. To overcome the drawback, the patch has been designed. It is obtained through modifying the patch to conform with the following conditions, which are necessary for the complete locality of the algorithm. (A)The and patches have the same boundary. (B)The normals of the patch on its edge depend only on the position and normal vectors at the endpoints of the edge.
After mathematical preliminaries, the normal of (B) is defined in Section 2.2.2. The patch is then determined from the boundary conditions.
2.2.1. Tangent Vector and Cross-Boundary Derivative along the Edges
For the subsequent formulation, directional derivative (see, e.g., [6, page 86]) becomes handy. Here is an arbitrary nonzero vector along the direction of the differentiation. Recalling (2.16) and (2.20), the partial derivatives with respect to the surface parameters and are represented by the above operator as Consider a triangular patch described by . Selecting the vector of (2.23) as in (2.30), the tangent vector along the edge can be written as Replacing the above with the vector for the next edge leads to the cross-boundary derivativeNote that and defined above are linear operators generating functions from their arguments. Vector product of (2.32a) and (2.32b) gives the normal of the patch along the edge .
For the patch of (2.24), its partial derivative with respect to the barycentric coordinate is Substituting this into (2.32a)-(2.32b) yields the tangent vector and the cross-boundary derivative aswhere the following constant vectors are introduced:The right-hand sides of (2.35a) are obtained through (2.25b). Equations (2.34a)-(2.34b) and (2.18a)-(2.18b) reveal that , represent the tangent vector and the cross-boundary derivative, respectively, at the vertex on the edge . Similarly, the tangent vectors of the edges and meeting at the vertex are and , respectively. Their normalized vector product gives the unit normal of the patch (2.24) at the vertex . By definition, it is perpendicular to the tangent vector (2.35a) and the cross-boundary derivative (2.35b) at the same location, that is,For ordinary cases where the patch reproduces the normals of the original surface at the vertices, Figure 2(a) indicates and hence it is unnecessary to evaluate (2.36) in practice. Moreover, if the surface has the unique normal at the vertex , it naturally becomes parallel to the one for the adjacent patch. In the subsequent discussion, such smooth surfaces are assumed to focus on the tangent-plane continuity.
2.2.2. Normals on the Edges
The most important feature of the patch summarized in Section 1 is complete locality; that is, positional continuity across shared edges is automatically assured without requiring the coefficients of the adjacent patches. This is an outcome of the fact that each edge curve is constructed only from the position and normal vectors at its endpoints. If the same information can also determine the tangent plane on the patch boundary, smoothness can be achieved retaining the locality.
According to the definition of (2.36), the unit normals at the endpoints and of the edge are and , respectively. These normals tend to be parallel when the edge becomes shorter, if the original surface is smooth. Therefore, the convergence is not lost by an assumption that the normal of the patch at an arbitrary point on the edge is a linear combination of and as where and are unknown functions. The above normal must be orthogonal to the tangent vector of (2.34a), and hence the following needs to be satisfied everywhere on the edge : Here (2.37a) is employed, and the constants are defined. Equation (2.40) can always be established if the unknown functions are selected as Here is an arbitrary constant, but the choice normalizes (2.42a) and hence is recommended. Based on (2.42a)-(2.42b) and (2.39), is adopted for the normal vector of the patch on the edge with the following coefficients: The special treatment of the case prevents the normal from vanishing. The corresponding coefficients make the normal of (2.43) equivalent to (2.39) with , . This satisfies the orthogonality condition of (2.40) and hence causes no problem.
If only one of and is zero, the normal of (2.43) may vanish at an endpoint of the edge. However, this does not imply difficulty because the normals at the vertices are predefined and hence unchanged after the interpolation; that is, the orthogonality is not necessary to be specified at the endpoints. Thus it is concluded that the proposed normal vector is adequate as long as the patch matches the boundary conditions.
2.2.3. Conditions on the Correction
The next step is to find a patch whose normal on its boundary is (2.43). Note that the equation depends only on the edge curve and the normals at its endpoints. Hence the normals automatically coincide on an edge shared by adjacent patches, and tangent-plane continuity is attained for ordinary cases where the original surface is smooth and its normals at the vertices are reproduced by the patch.
The patch in question is represented by the following function of the barycentric coordinate vector : On the right-hand side, is given by (2.24) describing the patch. Since the patch is generally not orthogonal to the normal of (2.43), the correction term is added. As mentioned in Section 2.2, this modification is carried out without changing the boundary curves. Therefore, the correction needs to vanish on the edges: This preserves the tangent vector of (2.34a), which is already perpendicular to because of (2.40). If the cross-boundary derivative is also normal to , that is, the tangent-plane continuity is achieved. The first term on the right-hand side can be simplified through (2.34b), (2.37b), and (2.43) as where the following constants are defined: The correction should be chosen to fulfill (2.46) and (2.48a)-(2.48b). For sufficient generality, an arbitrary rational function is assumed. Here , are polynomials, and the former is a vector function vanishing on the edges due to (2.46); that is, A directional derivative of (2.50) can be evaluated as and hence the cross-boundary derivative (2.32b) for the correction becomes where (2.18b) and (2.51) are employed. The above can be simplified further through applying (2.32b) as The above derivative is subject to the condition of (2.48a); that is, Multiplying the above by its denominator and substituting (2.48b) into the result yields which is necessary for the tangent-plane continuity, together with (2.51) demanding that vanish on the edges. Due to the boundary condition and (2.18b), must hold when at least one of the barycentric coordinates , , becomes zero. Therefore, can be written in the following form where is a polynomial: Once the polynomials and satisfying (2.56a)-(2.56b) are determined, the parametric representation of the patch is obtained through (2.50) and (2.45). Among an infinite number of such candidates, a simple function is designed in the subsequent consideration.
A monomial involved in the polynomial of (2.56b) has the partial derivative which gives the cross-boundary derivative through (2.32b) as Assuming that (2.56b) is at most quartic, the degree of on its right-hand side is one or zero. This gives where , , and are unknown coefficients. The cross-boundary derivative of the above polynomial is obtained through (2.58) as Substituting this and (2.43) into the condition (2.56a)-(2.56b) results in whose right-hand side is quartic or lower. Hence on the left-hand side is at most quadratic and can be described similarly to (2.24)as where and are constants. This leads to Substituting the above into (2.61) and comparing the coefficients arrives at the conditionwhich is required to reproduce the normal of (2.43). Advancing the suffixes of (2.64a) by one gives , whose ratio to (2.64b) is Unfortunately, the above equation cannot hold in general since the values of both sides are predefined by (2.44) and (2.49). This fact makes it difficult to realize (2.59) based on the hypothesis that is quartic or lower. Hence, a simple quintic polynomial is considered for in the next section.
2.2.4. Example of the Correction
Because of (2.58), a term in needs to be linear with respect to the barycentric coordinate to have nonzero cross-boundary derivatives on the edge . If the term is quadratic or higher in and , it does not affect the cross-boundary derivatives on the other edges. This observation motivates the choice of the quintic polynomial where , , and are constant unknown vectors. It is obvious that the above conforms with (2.56b). Equation (2.58) provides the cross-boundary derivative of (2.66) as Note that the above result has the coefficient only. Hence the corresponding term can exclusively adjust the normal on the edge . Substitution of (2.67) and (2.43) reduces the condition (2.56a) to Here the right-hand side consists of two scalar terms, and the unknown vector already has three components. Therefore, the polynomial on the left-hand side is arbitrary, provided that its structure is consistent with the other side. If (2.62) is assumed again, (2.63) and (2.18b) yield Substituting the above into (2.68) and comparing the coefficients lead to the following conditions: The value of in (2.70) is predefined by (2.49) and hence has no degree of freedom. Therefore, the following is generally required: Consequently, only the terms with the coefficient remain on the left-hand sides of (2.71). The coefficient is in proportion to the magnitude of on the right-hand side, and hence can have any nonzero value without loss of generality. The simplest choice is Substitution of (2.72a)-(2.72b) into (2.62) determines the denominator polynomial as and reduces (2.71) to the simple linear problem Since the above equation is underdetermined (two conditions for three unknowns), it is generally satisfied by an infinite number of candidates. Among them, there is always a unique minimum-error minimum-norm solution In general, minimization criteria assigned to redundant degrees of freedom have positive effects improving robustness and convergence; for example, (2.4) used to determine the curvature parameter is the key to the robustness of the patch [1]. Hence the above solution is adopted to complete the correction term through (2.73), (2.66), and (2.50). Formula () in [1] explicitly gives the generalized inverse on the right-hand side of (2.75) as Here (2.76b) relies on the fact that (2.36) is a unit vector. Note that , of (2.76c) are dimensionless due to the fact that , in (2.76b) are unit vectors, and that the coefficients , are normalized by (2.44) and (2.42b), and should be treated as zeros when their absolute values are below a small positive threshold , for avoiding numerical instability.
The proposed surface correction does not modify the boundary of the patch, and hence its geometric representation capability is rather limited. Therefore, too large magnitudes of the coefficient vectors , , imply failure in the aforementioned singularity detection or insufficient resolution of the original data. Since the coefficients have the dimension of length, such difficulties can be detected similarly to (2.6) by the condition where is a dimensionless constant. If the above holds, the next expression on the right-hand side of (2.76a) is used for evaluating the generalized inverse.
It is important to verify the consistency of the solution (2.75) with (2.74) for tangent-plane continuity. If the vectors , are linearly independent, the condition is fulfilled irrespective of the values on its left-hand side: this yields due to (2.76b), (2.76c) and hence corresponds to the first case in (2.76a). The solution always gives by (2.75) satisfying (2.74). The vectors , can be linearly dependent in the following two situations. (A).The tangent vector and the cross-boundary derivative at the vertex (see (2.34a)-(2.34b)) become perpendicular also to the normal at the vertex . The same is true if the vertices and are interchanged. Therefore, This reduces (2.75) to , which satisfies the condition (2.74).(B) and at least one of , is 0.Assume that . Then (2.49) yields and the first row of (2.74) gives . Since is nonzero, (2.74) can hold only if . Similarly, is required when .
The above discussion concludes that tangent-plane continuity can be attained except for the following two pathological cases (As mentioned at the end of Section 2.2.1, it is assumed that the patch reproduces the normals at its vertices.):These do not occur in meshes discretized properly. Assume (2.79a), for instance. Then (2.43) makes and parallel and hence the whole edge is on the tangent plane at the vertex . The plane, however, does not contain the cross-boundary derivative at the other vertex , due to (2.79a). This implies that the original curve of the edge has displacement away from the tangent plane near the vertex and cannot be approximated by the simple quadratic function (2.29). Interchanging and leads to the case of (2.79b). Hence the edge should be subdivided if (2.79a) or (2.79b) is encountered.
Due to (2.73), (2.17), and (2.15b), the denominator polynomial is positive on the patch excluding its vertices. Although it vanishes at the vertices (see (2.18a)), smoothness of the correction term can be maintained, as will be clarified in the following. Consider an arbitrary path approaching the vertex of (2.21a) in the space of barycentric coordinates. The position vector to a point on the path can be described as where is the distance from the vertex in the space. The unit vector is variable and arbitrary as long as the above satisfies (2.15b); that is, Substituting (2.80) into (2.73) yields which reveals that the denominator polynomial is first order with respect to . Similarly by (2.80), holds. If the ratio has a positive exponent of in the numerator, it vanishes as it approaches the vertex for irrespective of the path; that is (note that never vanishes in the denominator of (2.84a) due to (2.81)), The above with as well as (2.66) and (2.50) gives the limit of the correction term which does not depend on the path. Therefore, adopting it as the value of at the vertex to redefine (2.50) as preserves the continuity. This function is also smooth: substitution of into (2.52) yields the partial derivative whose limit vanishes at every vertex because of (2.84a)-(2.84b) with and the boundedness of (2.81); that is, This and (2.30) tell that every directional derivative also has the limit at the vertices. Furthermore, the limit coincides with the value obtained directly from the definition (2.30): where is an arbitrary vector. This clearly indicates that the correction never changes the normals at the vertices. Equation (2.89) concludes that the following extended definition of the directional derivative becomes continuous at every vertex: This assures the smoothness of the correction term of (2.86). The directional derivative on the right-hand side of (2.90) can be evaluated through substituting (2.87) into (2.30). The partial derivatives with respect to the surface parameters , can also be obtained similarly through (2.31).
2.2.5. Algorithm
Finally, the procedure for determining the patch is summarized below.(A)Give the position and normal vectors at the three vertices of the patch. (B)Apply the algorithm described in Section 2.1 to determine the coefficients of (2.13) for the patch and obtain the vectors through (2.25a)โ(2.25c). (C)Compute the constants of (2.35a)-(2.35b) and the unit normals of (2.36). (D)After calculating , through (2.44), (2.42b), and (2.41), evaluate (2.49) to get . If one of , is zero and , the tangent-plane continuity on the edge is not achieved. Then the edge should be subdivided since the case implies lack of mesh resolution (see the comment on (2.79a)-(2.79b)). (E)Determine the unknown vectors through (2.75) and (2.76a)โ(2.76c). Note that should be treated as zero if or (2.77) holds for positive dimensionless thresholds and . (F)Equations (2.73), (2.66), and (2.50) then give the correction . Adding it to the function of (2.24) for the patch yields the parametric representation of (2.45) for the patch.
3. Numerical Examples
Figure 3 shows a typical example of the proposed interpolation. Figure 3(a) is a triangular mesh of an ellipsoid with semiaxes of lengths 1, 2, and 3. It has only six vertices and hence is essentially an octahedron. Application of the patch to the mesh with the normals to the original surface at the vertices gives Figure 3(b). It is closer to the ellipsoid, but its tangent plane is not continuous. This problem is then solved by the patch algorithm of Section 2.2.5 yielding a completely smooth surface in Figure 3(c).
(a)
(b)
(c)
Next, a more general mesh of Figure 4(a) is tried. It is obtained through discretization of a solid model and has the position and normal vectors of the original surface at its vertices. The mesh, however, lacks all the other geometric information. The patch interpolates it as Figure 4(b). Note that the sharp features as well as smoothness of the original model are reproduced properly. This is due to the capability of the proposed approach inherited from the patch which enables handling an arbitrary number of normals at a vertex.
(a)
(b)
As already mentioned in Section 1, the present methodology has been developed mainly for computational science and engineering. Hence quantitative accuracy has the primary significance rather than visual effect. It is assessed in the following steps. First, meshes for a unit sphere, a right circular cone, and a cylinder with both unit base radius and height, and a torus whose major and minor radii are 2 and 1, respectively, are generated with diverse resolutions. They are processed by the proposed scheme, and geometric errors of the interpolation are evaluated. The result is shown in Figure 5. Both axes are in logarithmic scale. In the figure, the error is the maximum distance of the interpolated surface from the exact geometry, computed using a sufficient number of sample points (regardless of the mesh resolution: more than points for the cone; more than points for the cylinder and the sphere; more than for the torus, uniformly distributed on the parameter planes of the patches with all their boundaries involving the singular points, i.e., the apex of the cone and the sharp edges of the discal bases). The mesh scale is defined as the square root of the average patch area. The solid and hollow symbols represent the results by the and patches, respectively, both showing decrease in the error with the mesh scale , that is, convergence to the exact surface. The order of convergence for the patch varies between 2.2 and 4, depending on the types of the surfaces. Its decrease may be caused by the existence of isolated singular points (apexes) or saddle points. It is notable that the proposed algorithm maintains its stability even for extremely coarse meshes. This robustness may be attributed to the minimization criteria introduced in (2.4) and (2.75).
For most of the cases in Figure 5, the and patches attain the same accuracy. This implies that the maximum errors tend to occur on the patch boundaries shared by both. The discrepancies between the two methods observed for the cone are considered to be due to the existence of the apex. Since normals at singular points have arbitrariness, their definition can change the convergence rate. The performance of the interpolator can naturally be affected also by the mesh subdivision scheme to control the resolution.
Regular points on a surface can be classified into elliptic, parabolic, and hyperbolic points. Since the four models used in the error analysis cover all these types as well as typical singular features, the present technique is expected to have similar effectiveness for general geometries.
4. Concluding Remarks
The Nagata patch has been extended to establish tangent-plane continuity, retaining its complete locality and other merits. The proposed interpolation is able to stably obtain the same level of accuracy as that of the patch, which is already sufficient even for high-precision engineering, as illustrated in Section 1. Nevertheless, there still exist problems demanding smoothness preserved in discretized models. The patch will find applications in such areas. They include optical and contact analyses, where artificial normal discontinuities may give rise to unwanted jumps in refraction angles and friction forces, respectively.
As can be imagined from the varying nonintegral order of convergence observed in Figure 5, its theoretical estimation is difficult. Since the algorithm involves the generalized inverse, whose definition changes according to the rank of the matrix, the discussion would be fairly complex and hence deserves future study.
Similar to other interpolators, the performance of the present technique can be affected by the mesh density distribution. In this context, further research efforts will be directed also to adaptive mesh refinement with assured upper bound of the geometric errors due to the approximation by the proposed algorithm.