Updating meshes on deforming domains: An application ... - CiteSeerX
COMMUNICATIONS IN NUMERICAL METHODS IN ENGINEERING Commun. Numer. Meth. Engng (in press) Published online in Wiley InterScience (www.interscience.wiley.com). DOI: 10.1002/cnm.1013
Updating meshes on deforming domains: An application of the target-matrix paradigm‡ Patrick Knupp∗, † Optimization and Uncertainty Department, Sandia National Laboratories, M/S 1318, P.O. Box 5800, Albuquerque, NM 87185-1318, U.S.A.
SUMMARY A method for updating meshes in response to domain deformation is formulated within the target-matrix mesh optimization paradigm. By casting the problem within this paradigm, various formulations of the problem can be considered. Two local metrics are suggested for use in a multi-variable objective function to update meshes via numerical optimization. Both guarantee that if the deformation is null, the original mesh on the un-deformed domain will be replicated. Further, target matrices are constructed in a manner that ensures the updated mesh is as similar to the un-deformed mesh as possible. Numerical experiments confirm these properties on 2D examples. The update method has limitations, of course, the most important being that there is no guarantee that the deformed mesh will be invertible. Published in 2007 by John Wiley & Sons, Ltd. Received 11 May 2006; Revised 2 January 2007; Accepted 14 April 2007 KEY WORDS:
mesh optimization; deforming domains; shape optimization; deforming geometry; deforming meshes
1. INTRODUCTION The problem of updating meshes in response to domain deformation arises in many contexts including high energy deposition physics, projectile penetration studies, propellant burn, metals forging, fluid–structure interactions, and shape optimization for accelerator cavities. In this paper, the abstract problem is considered from a mesh optimization viewpoint. ∗ Correspondence
to: Patrick Knupp, Optimization and Uncertainty Department, Sandia National Laboratories, M/S 1318, P.O. Box 5800, Albuquerque, NM 87185-1318, U.S.A. † E-mail: [email protected] ‡ This article is a U.S. Government work and is in the public domain in the U.S.A. Contract/grant sponsor: Department of Energy’s Mathematics, Information, and Computational Science Program Contract/grant sponsor: Sandia Corporation; contract grant/number: DE-ACO4-94AL85000
Published in 2007 by John Wiley & Sons, Ltd.
P. KNUPP
The deforming-domain mesh-update problem is described as follows. One begins with an initial un-deformed domain !r (called the reference domain) and an original mesh Mr on the un-deformed domain that is assumed to be of good quality (i.e. it has the desired features for a particular application). Suppose the domain is deformed in some complex way, i.e. its shape, size, and or orientation changes, to create a deformed domain !d .§ How does one create a new mesh Md on !d that is also of good quality? Additionally, Md and Mr are often required to be ‘similar’ in some sense. Similarity might include, for example, identical mesh topology and/or the same size heterogeneity and anisotropy on both meshes. Similarity requirements can arise, for example, in domain deformations that occur between successive time steps within a simulation or between successive iterations in a shape-optimization procedure, wherein there is a need for smooth variation in the solution between deformations. The problem is, of course, very simple if one were given an ‘onto’ map X (!, t) from !r to !d at a specified time t. The mesh-update procedure in that case merely consists of evaluating the map at the vertices of the reference mesh, at the correct t. What if, however, such a map is unknown, not only on the interior, but even on the boundary? This is commonly the case in practice. One approach in this case would be to directly generate a mesh on the deformed geometry using whatever mesh generator is available to produce Md . Unless there were something special about the mesh generation scheme, however, there would be no guarantee in general of similarity between Md and Mr . In fact, this approach could result in Md " = Mr , even when !d = !r . To avoid that, one would be required to use exactly the same mesh generator and generation procedure to create both meshes. Even then, for unstructured meshes, preservation of topology would not be guaranteed when !d " = !r . Applications requiring automatic mesh updating may find a mesh generation approach problematic if the meshing requires human intervention to, for example, transfer mesh vertices on !!r to !!d . Even when the boundary map between the un-deformed and deformed domains is given, the mesh-updating problem is not entirely trivial due to the requirement of similarity. Methods for updating meshes on deforming domains when the boundary map is given have frequently involved ‘spring models’ related to the Laplace, variable diffusion, or biharmonic partial differential equations [1–5]. The requirement of fixed mesh topology strongly suggests mesh optimization procedures involving only node-movement be investigated. In this paper, a method for mesh updating on deforming domains where the boundary map is known is given within a mesh optimization framework. The method uses the idea of target matrices calculated from reference meshes. By placing the deformation method within the context of the target-matrix paradigm (TMP), all the methods of the paradigm are available to devise and study alternative formulations to improve techniques for addressing the mesh update problem on deforming domains.
2. BRIEF OVERVIEW OF THE TARGET PARADIGM A brief overview of the TMP for mesh optimization is described here; a more complete description of the paradigm is given in [6]. In TMP, optimal mesh vertex coordinates in #d , d = 2, 3, can be
§ In
the most general situation, the topology of !d could be different than !r . In this paper, it is assumed that the deformed domain has the same topology as the reference domain.
Published in 2007 by John Wiley & Sons, Ltd.
Commun. Numer. Meth. Engng (in press) DOI: 10.1002/cnm
UPDATING MESHES ON DEFORMING DOMAINS
calculated by minimization of a multi-variable objective function: F p,c,! {x} =
N 1 ! (cn !n ) p N n=1
where x = {xn ∈ #d |n = 1, . . . , N } ∈ #d N is a vector consisting of the mesh vertex coordinates. For a given mesh with coordinates x, M p,c,! {x} = (F p,c,! )1/ p is the ‘quality’ of the mesh. The quality is the power mean of the set of values "n = cn !n . In order for the power mean to be well defined, it is assumed that "n >0 for all n and that p " = 0. The parameter p in the mean determines the type of averaging: p = 1 yields the linear average of the "n , p = 2 is the root mean square (RMS), p = −1 is the harmonic average. In the limit, p = ∞ gives the maximum value over the set of "n , while p = −∞ gives the minimum. The scalar cn is found by evaluating a function c from #d to # at the point xn , i.e. cn = c(xn ). This scalar function is assumed to be known and serves as the mechanism for emphasizing or de-emphasizing terms in the sum defining the objective function. 2 The scalar !n is found by evaluating a function ! from #d to #; this scalar is the measure of the local mesh quality at particular locations within the mesh. Usually, ! is considered to be a scalar function of a d × d matrix. Let Md (#) be the set of real, d × d matrices and let T ∈ Md and assume ! = !(T ). The matrix T in the paradigm is the product of two matrices, A and W −1 , i.e. T = AW −1 . Matrix A is the Jacobian matrix of the map X from a point " within a master element to a point z within a given physical mesh element.¶ In particular, for d = 3, An = [!X/!"1 , !X/!"2 , !X/!"3 ] Matrix A is also referred to as the active matrix because it depends directly on values of local vertex (and possibly high-order node) coordinates; these coordinates are the coordinates of the mesh being optimized, i.e. they are components of x. The columns of A are the tangents to the map and thus can represent local mesh shape, size, and orientation. Matrix W is known as the target matrix because it represents (in matrix form), the target (or ideal) configuration of the local tangent vectors. In general, it is assumed that the ideal configuration arises from a non-singular map, thus det(W )>0 and so W −1 exists. The target matrices can, in general, be functions of space, thus Wn = W (xn ). The targets are assumed to be known a priori. Matrix T is known as the weighted-active or weighted-Jacobian matrix, and Tn = An (Wn )−1 . Finally, !n = !(Tn ). N is the total number of ‘sample points’ in the objective function, the sample points being the set of all the points within master elements at which ! is evaluated. Many variations on the ideas presented above are possible; some are described in [6], along with a complete analysis of the mathematical properties of the objects introduced here. The goal in the TMP is to shift focus away from absolute measures of mesh quality (based on geometrical properties of elements) towards measurement of quality relative (or with respect to) the set of target matrices. Ultimately, it is hoped that the target matrices, representing target or desired mesh quality, can be computed automatically using high-level user information such as physical solution Hessians, vector fields, PDE coefficients, surface curvature, and many other ¶ Note
that, even though local element maps are used in the TMP, it is not required that the application code use a finite element discretization.
Published in 2007 by John Wiley & Sons, Ltd.
Commun. Numer. Meth. Engng (in press) DOI: 10.1002/cnm
P. KNUPP
items not often used currently in mesh optimization methods. The targets may depend on the application goal (mesh smoothness, mesh shape improvement, mesh alignment, r -adaptivity) as well as the particular physics involved (elliptic, hyperbolic), and the precise details of the given application (magnitude of various quantities, the initial mesh, internal interfaces, etc.). In this paper, we describe the application of the TMP to the updating of mesh vertex coordinates in response to deformations of the domain boundary. In future papers, the TMP will be applied to entirely different application goals.
3. TMP APPLIED TO DEFORMING DOMAINS Assume that Mr and the boundary vertices of Md are known. Let Mdinit be the initial mesh on the deformed domain; this mesh consists of interior vertices having coordinates corresponding to the interior vertex coordinates of Mr , and boundary vertices having coordinates corresponding to the boundary vertex coordinates of Md . In general, Mdinit can have very poor quality and is not similar (at least on the boundary) to Mr ; this mesh can serve as the initial mesh in the proposed optimization procedure. Alternatively, any other mesh M init having the same mesh topology as Mr and the same boundary vertex coordinates as Mdinit can serve as the initial mesh. For example, M init could be the mesh resulting from applying Laplace smoothing to Mdinit with fixed boundary vertex coordinates. It is further assumed that the boundary vertices of M init form a boundary mesh that is ‘similar’ to the boundary mesh on !r . The assumption that the mesh boundary vertex coordinates of Md are known does not apply to many practical problems and must be dealt with; for the purpose of the present paper, these assumptions are made in order to focus on the problem of properly updating the interior vertex coordinates. A preliminary stab at the boundary problem is described in [7]. It is further assumed that boundary vertex positions are fixed in the optimization of the TMP objective function. To apply the TMP to the deforming domain application, targets W , local metrics !, scalar functions c, maps X, and sample points "n must be selected. For simplicity, assume cn = 1 for all n. Further assume the map X is affine for simplicial elements, bilinear for quadrilateral, and trilinear for hexahedral elements. It is then reasonable to select, as the sample points for quadrilateral or hexahedral elements, the ‘corners’ of the master element. For simplicial elements, there is one sample point, which can be located at any one corner of the master simplex.( The basic assumption of this method is that the original mesh on the undeformed domain is of good quality (as defined by the application) and thus, one wants to preserve this quality on the deformed domain. Therefore, targets will be calculated using the original mesh Mr as the reference. At each corner of the original mesh, one can compute the Jacobian matrix Aref assuming the maps above. These matrices then become the targets, i.e. Wn = (Aref )n and thus Tn = An (Aref )−1 n . A metric ! is then required such that, if the deformed geometry is identical to the undeformed geometry, the optimal mesh Md on the deformed geometry is the same as the mesh Mr on the undeformed geometry. In this case, An = (Aref )n for every n; hence Tn is the
( In
this case, the resulting matrix B in the affine map will be equal to the matrix T provided W is based on the master element.
Published in 2007 by John Wiley & Sons, Ltd.
Commun. Numer. Meth. Engng (in press) DOI: 10.1002/cnm
UPDATING MESHES ON DEFORMING DOMAINS
d × d identity matrix I . This observation suggests the following local quality metric: !def = |T − I |2F where F denotes the Frobenius matrix norm.∗∗ Note that !def !0, with !def = 0 if and only if T = I . If F p,1,!def is minimized, the optimal mesh will tend towards A = Aref , and thus the quality of the original mesh will be, to some extent, preserved in the optimal mesh. If no deformation of the boundary occurs, then the original mesh will be preserved exactly. It is possible to show that, except for degenerate cases, F p,1,!def is a strictly convex objective function since !def is a quadratic function of the local mesh vertex coordinates. The optimal mesh is therefore unique. Further, the local metric is well defined whether or not the local map is singular at a sample point. This is important because, except for very small deformations, the initial deformed mesh Mdinit (or even M init ) is likely to contain inverted elements. This feature is both an advantage and a disadvantage because, without a barrier, it cannot be guaranteed that the optimal mesh is non-inverted. This is an important limitation of this and other approaches. Mesh untangling methods (e.g. [8]), followed by optimization using barrier forms of !def , may sometimes solve the problem, but untanglers are not guaranteed to always work. Another important criticism of !def is that, although it can replicate the original mesh if there is no deformation of the domain boundary, it may not replicate if the domain boundary is re-scaled or rigidly rotated. To address this issue, one can consider the metric !)def , given by t 2 d =2 |T − (adj(T )) | F , ) !def = |T | F √ T − (adj(T ))t |2F , d = 3 3
where adj(T ) = det(T )T −1 . One can show that !)def = 0 if and only if T = s R where s>0 and R is an orthogonal matrix. Thus, when !)def = 0, A = s R Aref , i.e. the active-Jacobian matrix of the optimized deformed mesh will tend towards a scaled rotation of the original undeformed mesh Jacobian. 4. NUMERICAL EXAMPLES Examples are confined to d = 2 (two dimensions) in order to more easily visualize results. The case d = 3 (three dimensions) was used with success in [7, 9]. The case p = 1 was used for all of the examples. 4.1. Horseshoe domain To illustrate the basic ideas and how they work, the first numerical example uses a small structured mesh on a horseshoe geometry (see Figure 1). The mesh Mr on the initial un-deformed domain (northwest part of figure), is strongly biased, with smaller elements near the right domain boundary. ∗∗ |B|2 F
= tr(B t B).
Published in 2007 by John Wiley & Sons, Ltd.
Commun. Numer. Meth. Engng (in press) DOI: 10.1002/cnm
P. KNUPP
4
4
’ref_mesh.gpt’
3.5
3.5
3
3
2.5
2.5
2
2
1.5
1.5
1
1
0.5
0.5
0 -2
-1.5
-1
-0.5
0
4
0.5
1
1.5
2
0
3.5
3
3
2.5
2.5
2
2
1.5
1.5
1
1
0.5
0.5 -2
-1.5
-1
-0.5
0
0.5
1
1.5
-2
-1.5
-1
-0.5
0
0.5
4
’smoothed_mesh.gpt’
3.5
0
’ini_mesh.gpt’
2
0
1
1.5
2
’smo_mesh.gpt’
-2
-1.5
-1
-0.5
0
0.5
1
1.5
2
Figure 1. Biased-structured meshes on horseshoe domains: NW, undeformed; NE, initial deformed; SW, shape-optimized; and SE, optimized ! = !def .
The goal is to preserve this biasing and other properties under deformation. The northeast mesh was created by deforming the upper domain boundary upward, moving the left boundary vertices along with the deformation. This produces a mesh, Mdinit , that is not similar to Mr due to lack of smoothness. The southwest mesh was created by optimizing the NE mesh using a shape optimizer based on mean ratio (! = |A|2 /det(A)). Note that the result (SW) is still a non-smooth mesh and, in addition, that the biasing has been destroyed. The southeast mesh, Md , was created by optimizing the NE mesh using !def with target matrices calculated from the NW mesh as the reference. The result is smooth and properly biased, similar to the mesh on the un-deformed geometry. A similar experiment was performed on the horseshoe using a triangle mesh (see Figure 2). Once again, optimizing Mdinit (shown in the NE part of figure) with the mean-ratio shape-metric failed to smooth the mesh appreciably, although the biasing is preserved (SW). This is due to the fact that mean ratio uses targets W = I , indicating the ideal element is isotropic; the mesh elements in Figure 2, NW, are isotropic, whereas on the structured mesh (Figure 1, NW) they are anisotropic near the lower boundary. The mesh optimized with !def (Md ) is shown in Figure 2, SE; it preserves the properties of the NW mesh better than the SW mesh. Very large elements near the upper boundary in SE are still created, however. Published in 2007 by John Wiley & Sons, Ltd.
Commun. Numer. Meth. Engng (in press) DOI: 10.1002/cnm
UPDATING MESHES ON DEFORMING DOMAINS
4
4
’tri_horse12.gnu’
3.5
3.5
3
3
2.5
2.5
2
2
1.5
1.5
1
1
0.5
0.5
0
-2
-1.5
-1
-0.5
0
4
0.5
1
1.5
2
0
3.5
3
3
2.5
2.5
2
2
1.5
1.5
1
1
0.5
0.5 -2
-1.5
-1
-0.5
0
0.5
1
1.5
-2
-1.5
-1
-0.5
0
0.5
4
’smoothed_mesh.gpt’
3.5
0
’tri_horse12_dfn.gnu’
2
0
1
1.5
2
’smo_mesh.gpt’
-2
-1.5
-1
-0.5
0
0.5
1
1.5
2
Figure 2. Biased-unstructured meshes on horseshoe domains: NW, undeformed; NE, initial deformed; SW, shape-optimized; and SE, optimized ! = !def .
4.2. Airfoil domain with structured mesh and shape deformation A more realistic example is shown in Figure 3, where the mesh resulting from optimizing using !def is shown.†† On the northwest part of the figure, the optimized mesh around the deformed airfoil is shown. The northeast part of the figure shows the optimized mesh on the leading edge of the airfoil, while the southwest portion of the figure shows the optimized mesh on the trailing edge. The un-deformed domain was a RAE 2822 airfoil embedded within a large outer rectangle, at a zero-degree angle of attack. Unrealistically large deformations were created to stress the mesh optimization method. The original boundary layer around the airfoil is essentially preserved. Another example (not shown) on a two-airfoil problem was tried. In this case, it was possible to create inverted meshes by increasing the deformation sufficiently. The method in [1] also failed on this example [10].
††
Mesh and domain data courtesy of Tom Smith and Martin Berggren.
Published in 2007 by John Wiley & Sons, Ltd.
Commun. Numer. Meth. Engng (in press) DOI: 10.1002/cnm
P. KNUPP 0.3
0.04
’smo_mesh.gpt’
0.2
’smo_mesh.gpt’
0.02
0.1
0
0 -0.02 -0.1 -0.04
-0.2 -0.3 -0.2
0
0.2
0.4
0.6
0.8
0.1
1
1.2
-0.06 -0.1
-0.05
0
’smo_mesh.gpt’
0.05
0.1
’grid_eul_13725_mes.gnu’
0.4
0.08 0.06
0.2
0.04 0
0.02 0
-0.2
-0.02 -0.4
-0.04 0.8
0.85
0.9
0.95
1
1.05
1.1
-0.5
0
0.5
1
1.5
Figure 3. Optimized-structured mesh on deformed airfoil domain: NW, deformed mesh; NE, deformed leading edge; SW, deformed trailing edge (! = !def ); and SE, undeformed mesh.
4.3. Airfoil domain with unstructured mesh and rotational deformation In the example shown in Figure 4, a triangle mesh around an airfoil set within an ellipse is shown.‡‡ The deformation consists in this case of a simple change in the angle of attack of the airfoil with respect to the major axis of the ellipse; the un-deformed domain corresponds to the 0◦ angle of attack. The optimized meshes were created using !)def to account for the rotational deformation. In the northwest part of the figure, the optimized mesh for the 11◦ angle of attack is shown. In the northeast, the mesh at the trailing edge of the airfoil is shown for the 0◦ angle of attach (this is the reference mesh). In the southwest and southeast portions, respectively, the optimized mesh for the 11◦ and 22◦ angles of attack are shown. The triangles at the very tip of the trailing edge remain nearly equilateral as the angle of attack is increased. Meshes remained non-inverted up to a 45◦ angle of attack; data for higher angles of attack were unavailable.
‡‡
Mesh and domain data courtesy of Brian Dennis.
Published in 2007 by John Wiley & Sons, Ltd.
Commun. Numer. Meth. Engng (in press) DOI: 10.1002/cnm
UPDATING MESHES ON DEFORMING DOMAINS 6
0.1
’triangles.11norot.gnu’
’triangles.0norot.gnu’
4 0.05 2 0
0
-2 -0.05 -4 -6
-6
-4
-2
0
-0.1
2
4
6
-0.1 0.85
0.9
0.95
-0.25
’triangles.11norot.843RI.gnu’
1
1.05
’triangles.22norot.843RI.gnu’
-0.12 -0.14
-0.3
-0.16 -0.18
-0.35
-0.2 -0.22
-0.4
-0.24 -0.26
-0.45
-0.28 -0.3 0.75
0.8
0.85
0.9
0.95
1
1.05
-0.5 0.65
0.7
0.75
0.8
0.85
0.9
0.95
Figure 4. Optimized-unstructured meshes on deformed airfoil domains: NW, far-field, 11◦ ; NE, un-deformed; SW, rotation of 11◦ ; and SE, rotation of 22◦ .
5. CONCLUSIONS A method for optimizing meshes to update vertex coordinates in response to domain boundary deformations was presented within the context of the TMP. Two local metrics, !def and !)def , were proposed. In theory, the first is adequate for deformations not involving major re-scaling or rotation of the domain, while the second can work even when major re-scaling or rotation is involved. Limited numerical examples show that the method can work on domains with modest deformations. Larger deformations can result in inverted mesh elements. Future work should concentrate on adapting this approach to a method with guaranteed invertibility. This may be difficult since invertibility guarantees within the TMP context generally require local metrics with ‘barriers’ and, if the metric has a barrier,§§ then it cannot be used on an initially inverted mesh. For applications involving large deformations, the initial mesh is likely to be inverted. More challenging still, a basic assumption in this paper is that the boundary vertices of the mesh on the deformed domain are known; unless the deformation is given by a parametric description,
§§ See
[6] for a definition of barriers.
Published in 2007 by John Wiley & Sons, Ltd.
Commun. Numer. Meth. Engng (in press) DOI: 10.1002/cnm
P. KNUPP
the boundary vertex positions on the deformed geometry are unknown. This occurs, for example, when the un-deformed domain is given by a CAD geometry and the deformed domain by a second CAD geometry. This complication was encountered and dealt with in [7], using the objective function in this paper and a particular method for moving boundary vertices from the un-deformed to the deformed geometry. The method first re-meshed all the curves in the deformed geometry according to the curve spacing given on the un-deformed curves. Second, all surface meshes on the deformed geometry were optimized using !def while holding the positions of the bounding curve vertices fixed. Finally, the volume mesh on the deformed geometry was optimized using !def while holding the positions of the bounding surface vertices fixed. The deformations in that problem were relatively small, so it is unclear whether this approach to moving boundary vertices will work in general. Further research on this question is needed. ACKNOWLEDGEMENTS
This paper is SAND 2006-2885C. This work was funded by the Department of Energy’s Mathematics, Information, and Computational Science Program (SC-31) and was performed at Sandia National Laboratories. Sandia is a multiprogram laboratory operated by Sandia Corporation, a Lockheed Martin Company, for the U.S. Department of Energy under contract DE-ACO4-94AL85000. This work was performed by an employee of the U.S. Government or under U.S. Government contract. The U.S. Government retains a nonexclusive, royalty-free license to publish or reproduce the published form of this contribution, or allow others to do so, for U.S. Government purposes. REFERENCES 1. Berggren M. Edge-based Mesh Movement Strategies. Preprint, October 2003. 2. Helenbrook BT. Mesh deformation using the biharmonic operator. International Journal for Numerical Methods in Engineering 2003; 56:1007–1021. 3. Branets L, Carey GF. A local cell quality metric and variational grid smoothing algorithm. Engineering with Computers 2005; 21:19–28. 4. Shontz SM, Vavasis SA. A mesh warping algorithm based on weighted Laplacian smoothing. Proceedings of the 12th International Meshing Roundtable, Santa Fe, NM, September 2003; 147–158. 5. Stein K, Tezduyar T, Benney R. Mesh moving techniques for fluid-structure interactions with large displacements. Transactions of the ASME 2003; 70:58–63. 6. Knupp P. Formulation of a target-matrix paradigm for mesh optimization. SAND2006-2730J, Sandia National Laboratories, Albuquerque, NM, 2006. 7. Tautges T, Knupp P, Kraftcheck J, Kim H. Interoperable geometry and mesh components for SciDAC applications. Journal of Physics: Conference Series, 2005; 16:486–490. 8. Knupp P. Hexahedral and tetrahedral mesh untangling. Engineering with Computers 2001; 17(3):261–268. 9. Knupp P. Mesh quality improvement for SciDAC applications. Journal of Physics: Conference Series 2006; 46:458–462. 10. Berggren M. Personal communication, 2005.
Published in 2007 by John Wiley & Sons, Ltd.
Commun. Numer. Meth. Engng (in press) DOI: 10.1002/cnm
Updating meshes on deforming domains: An application of the target-matrix paradigm‡ Patrick Knupp∗, † Optimization and Uncertainty Department, Sandia National Laboratories, M/S 1318, P.O. Box 5800, Albuquerque, NM 87185-1318, U.S.A.
SUMMARY A method for updating meshes in response to domain deformation is formulated within the target-matrix mesh optimization paradigm. By casting the problem within this paradigm, various formulations of the problem can be considered. Two local metrics are suggested for use in a multi-variable objective function to update meshes via numerical optimization. Both guarantee that if the deformation is null, the original mesh on the un-deformed domain will be replicated. Further, target matrices are constructed in a manner that ensures the updated mesh is as similar to the un-deformed mesh as possible. Numerical experiments confirm these properties on 2D examples. The update method has limitations, of course, the most important being that there is no guarantee that the deformed mesh will be invertible. Published in 2007 by John Wiley & Sons, Ltd. Received 11 May 2006; Revised 2 January 2007; Accepted 14 April 2007 KEY WORDS:
mesh optimization; deforming domains; shape optimization; deforming geometry; deforming meshes
1. INTRODUCTION The problem of updating meshes in response to domain deformation arises in many contexts including high energy deposition physics, projectile penetration studies, propellant burn, metals forging, fluid–structure interactions, and shape optimization for accelerator cavities. In this paper, the abstract problem is considered from a mesh optimization viewpoint. ∗ Correspondence
to: Patrick Knupp, Optimization and Uncertainty Department, Sandia National Laboratories, M/S 1318, P.O. Box 5800, Albuquerque, NM 87185-1318, U.S.A. † E-mail: [email protected] ‡ This article is a U.S. Government work and is in the public domain in the U.S.A. Contract/grant sponsor: Department of Energy’s Mathematics, Information, and Computational Science Program Contract/grant sponsor: Sandia Corporation; contract grant/number: DE-ACO4-94AL85000
Published in 2007 by John Wiley & Sons, Ltd.
P. KNUPP
The deforming-domain mesh-update problem is described as follows. One begins with an initial un-deformed domain !r (called the reference domain) and an original mesh Mr on the un-deformed domain that is assumed to be of good quality (i.e. it has the desired features for a particular application). Suppose the domain is deformed in some complex way, i.e. its shape, size, and or orientation changes, to create a deformed domain !d .§ How does one create a new mesh Md on !d that is also of good quality? Additionally, Md and Mr are often required to be ‘similar’ in some sense. Similarity might include, for example, identical mesh topology and/or the same size heterogeneity and anisotropy on both meshes. Similarity requirements can arise, for example, in domain deformations that occur between successive time steps within a simulation or between successive iterations in a shape-optimization procedure, wherein there is a need for smooth variation in the solution between deformations. The problem is, of course, very simple if one were given an ‘onto’ map X (!, t) from !r to !d at a specified time t. The mesh-update procedure in that case merely consists of evaluating the map at the vertices of the reference mesh, at the correct t. What if, however, such a map is unknown, not only on the interior, but even on the boundary? This is commonly the case in practice. One approach in this case would be to directly generate a mesh on the deformed geometry using whatever mesh generator is available to produce Md . Unless there were something special about the mesh generation scheme, however, there would be no guarantee in general of similarity between Md and Mr . In fact, this approach could result in Md " = Mr , even when !d = !r . To avoid that, one would be required to use exactly the same mesh generator and generation procedure to create both meshes. Even then, for unstructured meshes, preservation of topology would not be guaranteed when !d " = !r . Applications requiring automatic mesh updating may find a mesh generation approach problematic if the meshing requires human intervention to, for example, transfer mesh vertices on !!r to !!d . Even when the boundary map between the un-deformed and deformed domains is given, the mesh-updating problem is not entirely trivial due to the requirement of similarity. Methods for updating meshes on deforming domains when the boundary map is given have frequently involved ‘spring models’ related to the Laplace, variable diffusion, or biharmonic partial differential equations [1–5]. The requirement of fixed mesh topology strongly suggests mesh optimization procedures involving only node-movement be investigated. In this paper, a method for mesh updating on deforming domains where the boundary map is known is given within a mesh optimization framework. The method uses the idea of target matrices calculated from reference meshes. By placing the deformation method within the context of the target-matrix paradigm (TMP), all the methods of the paradigm are available to devise and study alternative formulations to improve techniques for addressing the mesh update problem on deforming domains.
2. BRIEF OVERVIEW OF THE TARGET PARADIGM A brief overview of the TMP for mesh optimization is described here; a more complete description of the paradigm is given in [6]. In TMP, optimal mesh vertex coordinates in #d , d = 2, 3, can be
§ In
the most general situation, the topology of !d could be different than !r . In this paper, it is assumed that the deformed domain has the same topology as the reference domain.
Published in 2007 by John Wiley & Sons, Ltd.
Commun. Numer. Meth. Engng (in press) DOI: 10.1002/cnm
UPDATING MESHES ON DEFORMING DOMAINS
calculated by minimization of a multi-variable objective function: F p,c,! {x} =
N 1 ! (cn !n ) p N n=1
where x = {xn ∈ #d |n = 1, . . . , N } ∈ #d N is a vector consisting of the mesh vertex coordinates. For a given mesh with coordinates x, M p,c,! {x} = (F p,c,! )1/ p is the ‘quality’ of the mesh. The quality is the power mean of the set of values "n = cn !n . In order for the power mean to be well defined, it is assumed that "n >0 for all n and that p " = 0. The parameter p in the mean determines the type of averaging: p = 1 yields the linear average of the "n , p = 2 is the root mean square (RMS), p = −1 is the harmonic average. In the limit, p = ∞ gives the maximum value over the set of "n , while p = −∞ gives the minimum. The scalar cn is found by evaluating a function c from #d to # at the point xn , i.e. cn = c(xn ). This scalar function is assumed to be known and serves as the mechanism for emphasizing or de-emphasizing terms in the sum defining the objective function. 2 The scalar !n is found by evaluating a function ! from #d to #; this scalar is the measure of the local mesh quality at particular locations within the mesh. Usually, ! is considered to be a scalar function of a d × d matrix. Let Md (#) be the set of real, d × d matrices and let T ∈ Md and assume ! = !(T ). The matrix T in the paradigm is the product of two matrices, A and W −1 , i.e. T = AW −1 . Matrix A is the Jacobian matrix of the map X from a point " within a master element to a point z within a given physical mesh element.¶ In particular, for d = 3, An = [!X/!"1 , !X/!"2 , !X/!"3 ] Matrix A is also referred to as the active matrix because it depends directly on values of local vertex (and possibly high-order node) coordinates; these coordinates are the coordinates of the mesh being optimized, i.e. they are components of x. The columns of A are the tangents to the map and thus can represent local mesh shape, size, and orientation. Matrix W is known as the target matrix because it represents (in matrix form), the target (or ideal) configuration of the local tangent vectors. In general, it is assumed that the ideal configuration arises from a non-singular map, thus det(W )>0 and so W −1 exists. The target matrices can, in general, be functions of space, thus Wn = W (xn ). The targets are assumed to be known a priori. Matrix T is known as the weighted-active or weighted-Jacobian matrix, and Tn = An (Wn )−1 . Finally, !n = !(Tn ). N is the total number of ‘sample points’ in the objective function, the sample points being the set of all the points within master elements at which ! is evaluated. Many variations on the ideas presented above are possible; some are described in [6], along with a complete analysis of the mathematical properties of the objects introduced here. The goal in the TMP is to shift focus away from absolute measures of mesh quality (based on geometrical properties of elements) towards measurement of quality relative (or with respect to) the set of target matrices. Ultimately, it is hoped that the target matrices, representing target or desired mesh quality, can be computed automatically using high-level user information such as physical solution Hessians, vector fields, PDE coefficients, surface curvature, and many other ¶ Note
that, even though local element maps are used in the TMP, it is not required that the application code use a finite element discretization.
Published in 2007 by John Wiley & Sons, Ltd.
Commun. Numer. Meth. Engng (in press) DOI: 10.1002/cnm
P. KNUPP
items not often used currently in mesh optimization methods. The targets may depend on the application goal (mesh smoothness, mesh shape improvement, mesh alignment, r -adaptivity) as well as the particular physics involved (elliptic, hyperbolic), and the precise details of the given application (magnitude of various quantities, the initial mesh, internal interfaces, etc.). In this paper, we describe the application of the TMP to the updating of mesh vertex coordinates in response to deformations of the domain boundary. In future papers, the TMP will be applied to entirely different application goals.
3. TMP APPLIED TO DEFORMING DOMAINS Assume that Mr and the boundary vertices of Md are known. Let Mdinit be the initial mesh on the deformed domain; this mesh consists of interior vertices having coordinates corresponding to the interior vertex coordinates of Mr , and boundary vertices having coordinates corresponding to the boundary vertex coordinates of Md . In general, Mdinit can have very poor quality and is not similar (at least on the boundary) to Mr ; this mesh can serve as the initial mesh in the proposed optimization procedure. Alternatively, any other mesh M init having the same mesh topology as Mr and the same boundary vertex coordinates as Mdinit can serve as the initial mesh. For example, M init could be the mesh resulting from applying Laplace smoothing to Mdinit with fixed boundary vertex coordinates. It is further assumed that the boundary vertices of M init form a boundary mesh that is ‘similar’ to the boundary mesh on !r . The assumption that the mesh boundary vertex coordinates of Md are known does not apply to many practical problems and must be dealt with; for the purpose of the present paper, these assumptions are made in order to focus on the problem of properly updating the interior vertex coordinates. A preliminary stab at the boundary problem is described in [7]. It is further assumed that boundary vertex positions are fixed in the optimization of the TMP objective function. To apply the TMP to the deforming domain application, targets W , local metrics !, scalar functions c, maps X, and sample points "n must be selected. For simplicity, assume cn = 1 for all n. Further assume the map X is affine for simplicial elements, bilinear for quadrilateral, and trilinear for hexahedral elements. It is then reasonable to select, as the sample points for quadrilateral or hexahedral elements, the ‘corners’ of the master element. For simplicial elements, there is one sample point, which can be located at any one corner of the master simplex.( The basic assumption of this method is that the original mesh on the undeformed domain is of good quality (as defined by the application) and thus, one wants to preserve this quality on the deformed domain. Therefore, targets will be calculated using the original mesh Mr as the reference. At each corner of the original mesh, one can compute the Jacobian matrix Aref assuming the maps above. These matrices then become the targets, i.e. Wn = (Aref )n and thus Tn = An (Aref )−1 n . A metric ! is then required such that, if the deformed geometry is identical to the undeformed geometry, the optimal mesh Md on the deformed geometry is the same as the mesh Mr on the undeformed geometry. In this case, An = (Aref )n for every n; hence Tn is the
( In
this case, the resulting matrix B in the affine map will be equal to the matrix T provided W is based on the master element.
Published in 2007 by John Wiley & Sons, Ltd.
Commun. Numer. Meth. Engng (in press) DOI: 10.1002/cnm
UPDATING MESHES ON DEFORMING DOMAINS
d × d identity matrix I . This observation suggests the following local quality metric: !def = |T − I |2F where F denotes the Frobenius matrix norm.∗∗ Note that !def !0, with !def = 0 if and only if T = I . If F p,1,!def is minimized, the optimal mesh will tend towards A = Aref , and thus the quality of the original mesh will be, to some extent, preserved in the optimal mesh. If no deformation of the boundary occurs, then the original mesh will be preserved exactly. It is possible to show that, except for degenerate cases, F p,1,!def is a strictly convex objective function since !def is a quadratic function of the local mesh vertex coordinates. The optimal mesh is therefore unique. Further, the local metric is well defined whether or not the local map is singular at a sample point. This is important because, except for very small deformations, the initial deformed mesh Mdinit (or even M init ) is likely to contain inverted elements. This feature is both an advantage and a disadvantage because, without a barrier, it cannot be guaranteed that the optimal mesh is non-inverted. This is an important limitation of this and other approaches. Mesh untangling methods (e.g. [8]), followed by optimization using barrier forms of !def , may sometimes solve the problem, but untanglers are not guaranteed to always work. Another important criticism of !def is that, although it can replicate the original mesh if there is no deformation of the domain boundary, it may not replicate if the domain boundary is re-scaled or rigidly rotated. To address this issue, one can consider the metric !)def , given by t 2 d =2 |T − (adj(T )) | F , ) !def = |T | F √ T − (adj(T ))t |2F , d = 3 3
where adj(T ) = det(T )T −1 . One can show that !)def = 0 if and only if T = s R where s>0 and R is an orthogonal matrix. Thus, when !)def = 0, A = s R Aref , i.e. the active-Jacobian matrix of the optimized deformed mesh will tend towards a scaled rotation of the original undeformed mesh Jacobian. 4. NUMERICAL EXAMPLES Examples are confined to d = 2 (two dimensions) in order to more easily visualize results. The case d = 3 (three dimensions) was used with success in [7, 9]. The case p = 1 was used for all of the examples. 4.1. Horseshoe domain To illustrate the basic ideas and how they work, the first numerical example uses a small structured mesh on a horseshoe geometry (see Figure 1). The mesh Mr on the initial un-deformed domain (northwest part of figure), is strongly biased, with smaller elements near the right domain boundary. ∗∗ |B|2 F
= tr(B t B).
Published in 2007 by John Wiley & Sons, Ltd.
Commun. Numer. Meth. Engng (in press) DOI: 10.1002/cnm
P. KNUPP
4
4
’ref_mesh.gpt’
3.5
3.5
3
3
2.5
2.5
2
2
1.5
1.5
1
1
0.5
0.5
0 -2
-1.5
-1
-0.5
0
4
0.5
1
1.5
2
0
3.5
3
3
2.5
2.5
2
2
1.5
1.5
1
1
0.5
0.5 -2
-1.5
-1
-0.5
0
0.5
1
1.5
-2
-1.5
-1
-0.5
0
0.5
4
’smoothed_mesh.gpt’
3.5
0
’ini_mesh.gpt’
2
0
1
1.5
2
’smo_mesh.gpt’
-2
-1.5
-1
-0.5
0
0.5
1
1.5
2
Figure 1. Biased-structured meshes on horseshoe domains: NW, undeformed; NE, initial deformed; SW, shape-optimized; and SE, optimized ! = !def .
The goal is to preserve this biasing and other properties under deformation. The northeast mesh was created by deforming the upper domain boundary upward, moving the left boundary vertices along with the deformation. This produces a mesh, Mdinit , that is not similar to Mr due to lack of smoothness. The southwest mesh was created by optimizing the NE mesh using a shape optimizer based on mean ratio (! = |A|2 /det(A)). Note that the result (SW) is still a non-smooth mesh and, in addition, that the biasing has been destroyed. The southeast mesh, Md , was created by optimizing the NE mesh using !def with target matrices calculated from the NW mesh as the reference. The result is smooth and properly biased, similar to the mesh on the un-deformed geometry. A similar experiment was performed on the horseshoe using a triangle mesh (see Figure 2). Once again, optimizing Mdinit (shown in the NE part of figure) with the mean-ratio shape-metric failed to smooth the mesh appreciably, although the biasing is preserved (SW). This is due to the fact that mean ratio uses targets W = I , indicating the ideal element is isotropic; the mesh elements in Figure 2, NW, are isotropic, whereas on the structured mesh (Figure 1, NW) they are anisotropic near the lower boundary. The mesh optimized with !def (Md ) is shown in Figure 2, SE; it preserves the properties of the NW mesh better than the SW mesh. Very large elements near the upper boundary in SE are still created, however. Published in 2007 by John Wiley & Sons, Ltd.
Commun. Numer. Meth. Engng (in press) DOI: 10.1002/cnm
UPDATING MESHES ON DEFORMING DOMAINS
4
4
’tri_horse12.gnu’
3.5
3.5
3
3
2.5
2.5
2
2
1.5
1.5
1
1
0.5
0.5
0
-2
-1.5
-1
-0.5
0
4
0.5
1
1.5
2
0
3.5
3
3
2.5
2.5
2
2
1.5
1.5
1
1
0.5
0.5 -2
-1.5
-1
-0.5
0
0.5
1
1.5
-2
-1.5
-1
-0.5
0
0.5
4
’smoothed_mesh.gpt’
3.5
0
’tri_horse12_dfn.gnu’
2
0
1
1.5
2
’smo_mesh.gpt’
-2
-1.5
-1
-0.5
0
0.5
1
1.5
2
Figure 2. Biased-unstructured meshes on horseshoe domains: NW, undeformed; NE, initial deformed; SW, shape-optimized; and SE, optimized ! = !def .
4.2. Airfoil domain with structured mesh and shape deformation A more realistic example is shown in Figure 3, where the mesh resulting from optimizing using !def is shown.†† On the northwest part of the figure, the optimized mesh around the deformed airfoil is shown. The northeast part of the figure shows the optimized mesh on the leading edge of the airfoil, while the southwest portion of the figure shows the optimized mesh on the trailing edge. The un-deformed domain was a RAE 2822 airfoil embedded within a large outer rectangle, at a zero-degree angle of attack. Unrealistically large deformations were created to stress the mesh optimization method. The original boundary layer around the airfoil is essentially preserved. Another example (not shown) on a two-airfoil problem was tried. In this case, it was possible to create inverted meshes by increasing the deformation sufficiently. The method in [1] also failed on this example [10].
††
Mesh and domain data courtesy of Tom Smith and Martin Berggren.
Published in 2007 by John Wiley & Sons, Ltd.
Commun. Numer. Meth. Engng (in press) DOI: 10.1002/cnm
P. KNUPP 0.3
0.04
’smo_mesh.gpt’
0.2
’smo_mesh.gpt’
0.02
0.1
0
0 -0.02 -0.1 -0.04
-0.2 -0.3 -0.2
0
0.2
0.4
0.6
0.8
0.1
1
1.2
-0.06 -0.1
-0.05
0
’smo_mesh.gpt’
0.05
0.1
’grid_eul_13725_mes.gnu’
0.4
0.08 0.06
0.2
0.04 0
0.02 0
-0.2
-0.02 -0.4
-0.04 0.8
0.85
0.9
0.95
1
1.05
1.1
-0.5
0
0.5
1
1.5
Figure 3. Optimized-structured mesh on deformed airfoil domain: NW, deformed mesh; NE, deformed leading edge; SW, deformed trailing edge (! = !def ); and SE, undeformed mesh.
4.3. Airfoil domain with unstructured mesh and rotational deformation In the example shown in Figure 4, a triangle mesh around an airfoil set within an ellipse is shown.‡‡ The deformation consists in this case of a simple change in the angle of attack of the airfoil with respect to the major axis of the ellipse; the un-deformed domain corresponds to the 0◦ angle of attack. The optimized meshes were created using !)def to account for the rotational deformation. In the northwest part of the figure, the optimized mesh for the 11◦ angle of attack is shown. In the northeast, the mesh at the trailing edge of the airfoil is shown for the 0◦ angle of attach (this is the reference mesh). In the southwest and southeast portions, respectively, the optimized mesh for the 11◦ and 22◦ angles of attack are shown. The triangles at the very tip of the trailing edge remain nearly equilateral as the angle of attack is increased. Meshes remained non-inverted up to a 45◦ angle of attack; data for higher angles of attack were unavailable.
‡‡
Mesh and domain data courtesy of Brian Dennis.
Published in 2007 by John Wiley & Sons, Ltd.
Commun. Numer. Meth. Engng (in press) DOI: 10.1002/cnm
UPDATING MESHES ON DEFORMING DOMAINS 6
0.1
’triangles.11norot.gnu’
’triangles.0norot.gnu’
4 0.05 2 0
0
-2 -0.05 -4 -6
-6
-4
-2
0
-0.1
2
4
6
-0.1 0.85
0.9
0.95
-0.25
’triangles.11norot.843RI.gnu’
1
1.05
’triangles.22norot.843RI.gnu’
-0.12 -0.14
-0.3
-0.16 -0.18
-0.35
-0.2 -0.22
-0.4
-0.24 -0.26
-0.45
-0.28 -0.3 0.75
0.8
0.85
0.9
0.95
1
1.05
-0.5 0.65
0.7
0.75
0.8
0.85
0.9
0.95
Figure 4. Optimized-unstructured meshes on deformed airfoil domains: NW, far-field, 11◦ ; NE, un-deformed; SW, rotation of 11◦ ; and SE, rotation of 22◦ .
5. CONCLUSIONS A method for optimizing meshes to update vertex coordinates in response to domain boundary deformations was presented within the context of the TMP. Two local metrics, !def and !)def , were proposed. In theory, the first is adequate for deformations not involving major re-scaling or rotation of the domain, while the second can work even when major re-scaling or rotation is involved. Limited numerical examples show that the method can work on domains with modest deformations. Larger deformations can result in inverted mesh elements. Future work should concentrate on adapting this approach to a method with guaranteed invertibility. This may be difficult since invertibility guarantees within the TMP context generally require local metrics with ‘barriers’ and, if the metric has a barrier,§§ then it cannot be used on an initially inverted mesh. For applications involving large deformations, the initial mesh is likely to be inverted. More challenging still, a basic assumption in this paper is that the boundary vertices of the mesh on the deformed domain are known; unless the deformation is given by a parametric description,
§§ See
[6] for a definition of barriers.
Published in 2007 by John Wiley & Sons, Ltd.
Commun. Numer. Meth. Engng (in press) DOI: 10.1002/cnm
P. KNUPP
the boundary vertex positions on the deformed geometry are unknown. This occurs, for example, when the un-deformed domain is given by a CAD geometry and the deformed domain by a second CAD geometry. This complication was encountered and dealt with in [7], using the objective function in this paper and a particular method for moving boundary vertices from the un-deformed to the deformed geometry. The method first re-meshed all the curves in the deformed geometry according to the curve spacing given on the un-deformed curves. Second, all surface meshes on the deformed geometry were optimized using !def while holding the positions of the bounding curve vertices fixed. Finally, the volume mesh on the deformed geometry was optimized using !def while holding the positions of the bounding surface vertices fixed. The deformations in that problem were relatively small, so it is unclear whether this approach to moving boundary vertices will work in general. Further research on this question is needed. ACKNOWLEDGEMENTS
This paper is SAND 2006-2885C. This work was funded by the Department of Energy’s Mathematics, Information, and Computational Science Program (SC-31) and was performed at Sandia National Laboratories. Sandia is a multiprogram laboratory operated by Sandia Corporation, a Lockheed Martin Company, for the U.S. Department of Energy under contract DE-ACO4-94AL85000. This work was performed by an employee of the U.S. Government or under U.S. Government contract. The U.S. Government retains a nonexclusive, royalty-free license to publish or reproduce the published form of this contribution, or allow others to do so, for U.S. Government purposes. REFERENCES 1. Berggren M. Edge-based Mesh Movement Strategies. Preprint, October 2003. 2. Helenbrook BT. Mesh deformation using the biharmonic operator. International Journal for Numerical Methods in Engineering 2003; 56:1007–1021. 3. Branets L, Carey GF. A local cell quality metric and variational grid smoothing algorithm. Engineering with Computers 2005; 21:19–28. 4. Shontz SM, Vavasis SA. A mesh warping algorithm based on weighted Laplacian smoothing. Proceedings of the 12th International Meshing Roundtable, Santa Fe, NM, September 2003; 147–158. 5. Stein K, Tezduyar T, Benney R. Mesh moving techniques for fluid-structure interactions with large displacements. Transactions of the ASME 2003; 70:58–63. 6. Knupp P. Formulation of a target-matrix paradigm for mesh optimization. SAND2006-2730J, Sandia National Laboratories, Albuquerque, NM, 2006. 7. Tautges T, Knupp P, Kraftcheck J, Kim H. Interoperable geometry and mesh components for SciDAC applications. Journal of Physics: Conference Series, 2005; 16:486–490. 8. Knupp P. Hexahedral and tetrahedral mesh untangling. Engineering with Computers 2001; 17(3):261–268. 9. Knupp P. Mesh quality improvement for SciDAC applications. Journal of Physics: Conference Series 2006; 46:458–462. 10. Berggren M. Personal communication, 2005.
Published in 2007 by John Wiley & Sons, Ltd.
Commun. Numer. Meth. Engng (in press) DOI: 10.1002/cnm
