Novel applications of BEM based Poisson level set approach

Engineering Analysis with Boundary Elements 36 (2012) 907–912

Contents lists available at SciVerse ScienceDirect

Engineering Analysis with Boundary Elements journal homepage: www.elsevier.com/locate/enganabound

Novel applications of BEM based Poisson level set approach H. Xia a,n, P.G. Tucker b, G. Coughlin b a b

Thermo-Fluid Mechanics Research Centre, University of Sussex, Brighton, BN1 9QT, UK Whittle Laboratory, University of Cambridge, 1 JJ Thomson Avenue, Cambridge CB3 0DY, UK

a r t i c l e i n f o

abstract

Article history: Received 27 April 2010 Accepted 28 July 2011 Available online 4 October 2011

Accurate and efficient computation of the distance function d for a given domain is important for many areas of numerical modeling. Partial differential (e.g. Hamilton–Jacobi type) equation based distance function algorithms have desirable computational efficiency and accuracy. In this study, as an alternative, a Poisson equation based level set (distance function) is considered and solved using the meshless boundary element method (BEM). The application of this for shape topology analysis, including the medial axis for domain decomposition, geometric de-featuring and other aspects of numerical modeling is assessed. & 2011 Elsevier Ltd. All rights reserved.

Keywords: Boundary element method Distance function Medial axis De-featuring

1. Introduction Distance level set, d, is a key parameter in many (numerical) simulation approaches [1], in peripheral applications incorporating additional solution physics [2] and also in mesh generation [3]. It has proved important in computer vision, solid modeling and other computational science. Especially, distance function can be helpful in constructing the medial axis transform (MAT) for a given geometry [4]. The latter is regarded as a key step in shape analysis and solid modeling [5,6]. In this paper, as a further study to [4], we propose a hybrid MAT approach based on the Poisson distance function, namely dP -MAT. The advantage of such a hybrid method being to extract a well approximated medial axis point cloud on a properly calculated distance field has been discussed in detail in Refs. [4,7]. A simple Laplacian or Hessian determinant criterion of the distance field can be applied to mark the medial axis, and optionally combined with the a-shape thinning/representation techniques to thin the marked area to mathematically thin curves. Results show that the Poisson distance is very competitive in predicting wall proximity and at the same time efficient. Another application of dP is for geometric de-featuring. In computer aided design (CAD), the model of a part is often used for analysis. These aesthetically pleasing models contain much manufacturing information that is superfluous to analysis requirements. It would be highly advantageous to easily alter the CAD model to suit analysis purposes. De-featuring, or feature suppression, involves removing small, detailed information from the CAD model. These small features increase the density of the mesh, thus

n

Corresponding author. E-mail address: [email protected] (H. Xia).

0955-7997/$ - see front matter & 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.enganabound.2011.07.013

increasing the complexity of both the mesh generation and simulation, without adding significant information to the solution. Also, it is an extremely time-consuming and difficult process to remove features by hand. There are a number of methods so far seen in the literature [8–10]. The crucial part here is the automation. The possibility of using distances for automatic geometrically shape feature removing is explored here with the Poisson distance function. In this study, the focus is on the solution of a Poisson distance function equation within the framework of the meshless boundary element methods. First, we will discuss the Poisson equation of an auxiliary scalar r2 c ¼ r, where c is linked to the distance function d by a quadratic equation. The attraction of applying BEM here is driven by efficiency. Without meshing and integrating in the interior elements, the computation can usually be reduced by one dimension. That is to say, if the original distance problem is fully 3D, with BEM we only need to solve a 2D problem. It is also advantageous in moving boundary problems, where other discretization methods would require the entire interior domain to be re-meshed. However, as noted by Ramsak and Skerget in their recent work [11], there are not many efficient BEM formulations for a generic Poisson equation in the literature. For example, Suciu et al. [12] used the Galerkin vector approach but limited the source function r to satisfy r2 r ¼ const rather than zero, i.e. harmonic. Fortunately, as will be shown later, in the present study the auxiliary Poisson equation of c in the distance function context can be actually reduced to a Laplacian equation due to the fact that the inhomogeneous term r can be chosen as a constant without affecting the distance solution. This paper is organized as follows. The Poisson equation of c is introduced in the next section followed by the numerical solutions and test cases. The dP -MAT approach for medial axes are

908

H. Xia et al. / Engineering Analysis with Boundary Elements 36 (2012) 907–912

discussed in Section 4. Finally, the application of distance functions for geometric de-featuring is discussed.

2. Derivation of Poisson distance function In contrast to the eikonal exact distance governing equation, other differential (and indeed integral) equations describing the ‘distance’ function are also possible. Several of these methods were specifically sought to find the wall proximity for turbulence modeling [13]. One of these approaches goes back to the work of Spalding [14], in which the Poisson equation of an auxiliary variable c is solved. Consider the ultimately simple case, 1D distance from the center point, and suppose there is a scalar function cðxÞ, such that d2 c 2

dx

¼ 1

and

cð0Þ ¼ 0

ð1Þ

After integration twice, it becomes dc þxC0 ¼ 0 dx

and

Although Eqs. (1) and (4) are exact, the extension (6) to multidimension is not. However, for most numerical modeling problems considered here, only the near wall accuracy matters. Another advantage of Eqs. (5) and (6) is that they overestimate d around sharp convex surfaces and underestimate it around concave. This in exactness has consistent traits with the solid angle based integral equation for distance function equation proposed by Launder et al. [13] and Spalart [15] intended for improving the modeling of turbulence model physics and potentially other aspects of numerical modeling.

3. Numerical solution methods 3.1. Fast multipole BEM solution

x2 c þ C0 x þ C1 ¼ 0 2

ð2Þ

Combining the above and the boundary condition in (1), we can arrive at the quadratic equation of x: 1 2 dc x þ xc ¼ 0 2 dx

where the source function r ¼ 1, and the distance can be approximated as ffi   sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  2 @c @c d  n ¼   7 þ2c ð6Þ @n @n

ð3Þ

and the solution x, i.e. the distance function d, can be expressed as sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  2 dc dc 7 d¼x¼ þ 2c ð4Þ dx dx where the positive solution is meaningful. It then can be extended to multi-dimension with dc=dx replaced by @c=@n where n is the outward unit normal, and the governing equation for cðxÞ becomes ( 2 r c ¼ r, x A O ð5Þ c ¼ 0, x A @O

The solution of the Poisson equation of the domain integral has been formulated using the Galerkin vector approach. It requires the knowledge of a fundamental solution and is valid only for harmonic source terms [12]. Fundamental solutions are known for some, but not all, differential equations. Other methods such as the dual reciprocity method (DRM) [16] and the multiple reciprocity method (MRM) [17] may be used for more complex problems. However, for the current special case, where r ¼ const, the Poisson equation (5) can be reduced to a Laplacian equation and solved using the efficient fast multipole method (FMM) by Liu and Nishimura [18]. Let f be a new function, such that

f ¼ c þ f ðx,yÞ Thus Eq. (5) becomes ( 2 r f ¼ 0, ðx,yÞ A O

f ¼ f ðx,yÞ, ðx,yÞ A @O

Fig. 1. The 2D d contours for three domains. Upper: Poisson. Lower: exact.

ð7Þ

ð8Þ

H. Xia et al. / Engineering Analysis with Boundary Elements 36 (2012) 907–912

where f ðx,yÞ is not unique and we choose the simplest form, for example f ðx,yÞ ¼ 14ðx2 þ y2 Þ. The above Laplacian equation can be solved using the conventional BEM with the FMM acceleration. The main idea of the FMM is to translate the node-to-node interactions into adaptive Cartesian cell-to-cell ones. GMRES is used for iteratively finding the resultant linear system’s solution. In general, the solution time is O(N) using FMM, where N is the number of nodes on @O. More specific details can be found, for example, in Ref. [18]. 3.2. Finite volume solution Secondly, for comparison, Eq. (5) can also be solved using the finite volume method (FVM). By introducing a pseudotemporal term @c=@t, one arrives at an analogous diffusion equation of c. The finite volume formula can then obtained after applying a further integral over an arbitrary control volume V bounded by @V: Z Z @c dV ¼ ðr2 c1Þ dV ð9Þ V @t V where t is the pseudotime. This is implemented in an existing inhouse FVM code for general CFD purposes, in which the multigrid and parallelism are available. 3.3. Test cases Three 2D domains are selected for numerical tests. Fig. 1 shows the 2D distance contours for the Poisson (upper frames)

909

and exact solutions for these (lower frames). As discussed earlier, the Poisson equation tends to smooth out some curvature features by increasing d at convex corners and decreasing d at concave corners. To better demonstrate this trend, 3D distribution of the distances are shown in Fig. 2. Comparing the upper row with the lower, we can see that the ‘roof top’ shapes from the exact solutions are smeared and smoothed by the Poisson equation. The overall trend is well captured. To help compare quantitatively the solution, Fig. 3a defines point locations from P1 to P5 and two slices at X ¼ 11; 16. As shown in Table 1, at P2 which is close to the boundary and far away from convex/concave features, the Poisson distance is accurate. As the point location moves into the midfield, such as P1 and P3, the Poisson distance is less accurate and under-predicted. When the point location moving towards the sharp corners, such P4 and P5, it over- and under-estimates the distance considerably. However, sometimes as noted earlier,

Table 1 Comparison of distance solutions at various locations defined in Fig. 3a. Location

Poisson distance

P1 P2 P3 P4 P5

Exact d

BEM

FVM

2.28 0.20 3.90 0.31 0.31

2.28 0.20 3.90 0.31 0.31

Fig. 2. The d distributions plotted in 3D. Upper: Poisson. Lower: exact.

6 Exact Possion

10 P5

5

+

4

P3

+

P1

X=16

d

Y

+

P2

0

+

+

P4

2 X=11

-5 0

5

10 X

15

20

0

0

2

4

6

8

Y

Fig. 3. Comparison of the Poisson and exact distances. (a) Slice and point locations. (b) Profiles along slices.

2.7 0.2 4.5 0.2 0.4

910

H. Xia et al. / Engineering Analysis with Boundary Elements 36 (2012) 907–912

this is quite desirable [19]. Fig. 3b compares distances at two X locations. Again, the Poisson solutions predict well in the vicinity of the wall boundary but smooth out the sharp roof top features of the distribution.

4. Medial axes using Poisson distance: dP -MAT The medial axis of a shape provides a compact representation of its features and their connectivity. As a result, researchers have discovered and are still exploring its use in many fields, such as topology recognition for grid generation [20]. The medial axis is defined when the shape is embedded in an Euclidean space and is endowed with a distance function. Therefore, an expedient route is to efficiently obtain d (as discussed in previous sections) followed by the medial axis construction. (Notice that this approach is different from the pure geometric Voronoi-diagram based approaches, in which the equal distance nature of the Voronoi-diagram is a key.) In space, a sphere is called medial if it meets S, the domain boundary, only tangentially at at least two points. The medial axis M is defined as the closure of the set of centers of all medial spheres. Informally, the medial axis of a surface is the set of all points that have more than one closest point on the surface. They are often called the medial axis transform (MAT) for that bounded domain. The key step for the MAT construction here is to detect the medial axis feature in a given d field. The work in Refs. [4,7] already addressed this for the exact d. In this study, our aim is to explore the use of the more efficiently evaluated Poisson distance, even though the Poisson d can compromise the accuracy in some regions, but broadly for topology recognition and domain decomposition this lower accuracy can be neglected. Notice that the unique property of a medial axis point is that it has equal distance to multiple boundaries. Hence, in x2d space, the medial axis represents the ‘local maxima’ or non-smoothness of the distance function dðxÞ. For example, in one dimension the medial axis is a single point. The Poisson MAT does give a reasonably good approximation of the exact MAT, as shown in Fig. 4. Another advantage of dP -MAT, as discovered by Xia and Tucker [4], is solution superpositioning to create customized MATs. This is also straightforward for Poisson distance functions.

Poisson distance function, i.e. over- and under-estimate the distance for concave and convex features. The more d increases from the solid boundary, the more geometric features are removed. Therefore, d contours in a domain can be seen as smoothing out features with increasing values. We can obtain a de-featured geometry by extracting a contour level and projecting it back by amount dsm, where dsm is the (nominal) distance in smooth regions, to preserve the original global dimension. The normal at a local point on one particular contour may be expressed as, n ¼ rd=JrdJ. This is used to determine the contour’s backward-projected new position by x0 ¼ xdn

ð10Þ

Strictly speaking, an ‘entropy preserving’ level set type approach should be taken into account for the backward projection. Since the offset dsm is usually small, the above direct offset approach Eq. (10) works well in practice. In Fig. 6, the original and Poisson based de-featured model for an idealized electronic board is shown. The figure clearly demonstrates that with higher value of dsm features tend to be smoothed out or suppressed. Fig. 7 shows the original and feature suppressed models for the inside of a complex mechatronic system. As can be seen from frame (a) the original geometry is extremely complex. Note the de-featuring in frames (b) and (c) are based on exact wall distance which is less aggressive than the Poisson. Notably, for both types of model, the role of such de-featuring is to suppress unnecessary complexity for CFD analysis aimed as establishing, for example, heat transfer. Such a de-featuring process is computationally efficient and, more importantly, automatic. Fig. 8 presents a cross section of the mechatronic system comparing de-featuring using Poisson and standard distance functions. The round and more aggressive de-featuring of Poisson distance is obvious, and also is much faster to evaluate. De-featuring can greatly reduce the complexity of the model and the number of elements required for generating a quality mesh. Fig. 9a shows the Cartesian mesh generated for a defeatured model containing only 75,000 grid elements whereas the full model mesh has 2.5 million elements. The flow solution is shown in frame (b) where the velocity contour and streamlines are highlighted.

5. Distance function for geometric de-featuring De-featuring is a useful tool for analysis as it removes small features from the model that would otherwise increase mesh density, and thus computation time, without adding significant information to the simulation. Here, it is proposed to use the Poisson distance function to remove superfluous features. Fig. 5 shows the schematic of solid boundary and two different levels of the distance d. It demonstrates the characteristics of the

Fig. 5. Schematic of different distance levels suppressing the geometric features.

Fig. 4. Poisson distance based medial axes (symbols) vs. exact solutions (dashed lines).

H. Xia et al. / Engineering Analysis with Boundary Elements 36 (2012) 907–912

911

Fig. 6. An idealized electronic board. (a) Original geometry. (b) De-featured using dsm ¼ 0:02. (c) De-featured using dsm ¼ 0:30.

Fig. 7. A real CAD model for the inside of a mechatronic system. (a) Original model. (b) De-featured using dsm ¼ 0:035. (c) De-featured using dsm ¼ 0:5.

Fig. 8. Comparison of Poisson and standard wall distance de-featuring on a 2D section of the mechatronic model. (a) Original. (b) Poisson dsm ¼ 0:005. (c) Poisson at dsm ¼ 0:035. (d) Standard dsm ¼ 0:035.

912

H. Xia et al. / Engineering Analysis with Boundary Elements 36 (2012) 907–912

flow Fig. 9. Flow solution for the de-featured model.

6. Conclusions In this work, we explore use of a Poisson PDE based distance function (level set) for medial axis transform, de-featuring and other aspects of numerical modeling. This equation is based on a constant inhomogeneous source term. Therefore, the Poisson equation can be further reduced to a Laplacian equation. Solutions using the fast multipole method based on the conventional BEM methods are explored and compared with finite volume approaches. On computational efficiency ground the results are encouraging. For the special application of medial axis transform, the Poisson distance is also found equivalent to the exact distance based dP -MAT approaches. As also shown, the Poisson distance function has useful properties for geometric de-featuring and suppression. Its potential for the mathematical modeling of turbulence is also noted.

References [1] Tucker PG. Differential equation-based wall distance computation for DES and RANS. J Comput Phys 2003;190(1):229–48. [2] Osher S, Sethian JA. Fronts propagating with curvature dependent speed: algorithms based on Hamilton–Jacobi formulations. J Comput Phys 1988;79: 12–49. [3] Sethian J. Level set methods and fast marching methods: evolving interfaces in computational geometry, fluid mechanics, computer vision, and materials science. Cambridge University Press; 1999. [4] Xia H, Tucker PG. Finite volume distance field and its application to medial axis transforms. Int J Numer Meth Eng 2010;82(1):114–34.

[5] Price MA, Armstrong CG. Hexahedral mesh generation by medial surface subdivision: part ii. Solids with flat and concave edges. Int J Numer Meth Eng 1997;40:111–36. [6] Quadros WR, Owen SJ, Brewer M, Shimada K. Finite element mesh sizing for surfaces using skeleton. in: 13th international meshing roundtable; 2004. p. 389–400. [7] Xia H, Tucker PG. Distance solutions for medial axis transform. In: 18th international meshing roundtable, Salt Lake City, Utah, USA; 2009. [8] Dey S, Shephard MS, Georges MK. Elimination of the adverse effects of small model features by the local modification of automatically generated meshes. Eng Comput 1997;13(3):134–52. [9] Lee Y-G, Lee K. Geometric detail suppression by the Fourier transform. Comput-Aided Des 1998;30:677–93. [10] Lee YK, Lim CK, Ghazialam H, Vardhan H, Eklund E. Surface mesh generation for dirty geometries by shrink wrapping using Cartesian grid approach. in: 15th international meshing roundtable. Springer-Verlag; 2006. p. 393–410. [11] Ramsak M, Skerget L. 3d multidomain BEM for a Poisson equation. Eng Anal Boundary Elem 2009;33(5):689–94. [12] Suciu R, De Mey G, De Baetselier E. Bem solution of Poisson’s equation with source function satisfying r2 r ¼ constant. Eng Anal Boundary Elem 2001;25(2):141–5. [13] Launder BE, Reece GJ, Rodi W. Progress in the development of a Reynoldsstress turbulence closure. J Fluid Mech 1975;68(3):537–66. [14] Spalding DB. Poisson equation for wall distance. in: 10th international heat transfer conference Brightion, UK; 1994. [15] Spalart PR, Allmaras SR. A one-equation turbulence model for aerodynamic flows. AIAA Paper 92-0439; January 1992. [16] Partridge PW, Brebbia CA, Wrobel LC. The dual reciprocity boundary element method. Computational Mechanics Publications; 1992. [17] Nowak AJ, Brebbia CA. The multiple-reciprocity method. A new approach for transforming BEM domain integrals to the boundary. Eng Anal Boundary Elem 1989;6(3):164–7. [18] Liu YJ, Nishimura N. The fast multipole boundary element method for potential problems: a tutorial. Eng Anal Boundary Elem 2006;30(6):371–81. [19] Tucker PG, Liu Y. Turbulence modeling for flows around convex features giving rapid eddy distortion. Int J Heat Fluid Flow 2007;28:1073–91. [20] Rigby DL. Topmaker: a technique for automatic multi-block topology generation using the medial axis. NASA/CR 2004-213044; April 2004.