Abstract:
A method is disclosed for partitioning an input polytope into a conformal set of polytopes according to a subdivision pattern specified on lower-dimensional boundaries of the input polytope. A feature of this method is that it is applicable to any dimension. For input polytopes that are hexahedra isomorphic to a 3-dimensional cube, all members of the conforming partition will be hexahedra isomorphic to a 3-dimensional cube. A further feature of the invention is that it preserves convexity; a convex input polytope will result in a conformal partition comprising convex polytopes. Finally, the method may be used to process large conformal meshes of input polytopes in parallel and produce a globally conformal partition without the requirement of communication between processes.

Description:
CROSS REFERENCE TO RELATED APPLICATIONS 
     This application claims priority to prior co-pending U.S. Provisional Patent Application Ser. No. 61/058,802 originally filed Jun. 4, 2008 entitled “System and Method for Polytopic Mesh Refinement” from which benefit is claimed. 
    
    
     STATEMENT OF GOVERNMENT SUPPORT 
     The United States Government has a paid-up license in this invention and the right in limited circumstances to require the patent owner to license to others on reasonable terms as provided for by the terms of contract DE-AC04-94AL85000 awarded by the U.S. Department of Energy to Sandia Corporation. 
    
    
     TECHNICAL FIELD 
     This invention relates to geometric modeling; scalar, vector, and tensor field approximation including finite-element and finite-volume methods, and particularly including automatic spatial refinement of complexes of polytopes. The case where the polytopes are all hexahedra and refinement is specified by marking edges to be subdivided is the preferred embodiment. 
     BACKGROUND 
     A common problem in engineering is the design of a device whose operation is governed by physical phenomena too complex or whose shape is too irregular for a closed-form analytical expression describing its behavior to be developed. In these cases, it is typical for a finite element or finite volume approximation to be created by decomposing the geometric shape of the device and/or the space around the device into a number of primitive shapes over which a simple analytical form is used to approximate some quantity of interest (e.g., stress, strain, temperature). These primitive shapes may be simplices such as lines, triangles, or tetrahedra; or non-simplicial polytopes such as quadrilaterals, hexahedra, pyramids, wedges, or knives. While the analytical forms defined over each primitive are simple, their combination may exhibit complex patterns that closely approximate the solution to some engineering problem. Unfortunately, if the decomposition of the geometry into primitives is too coarse in some regions, the variations allowed by the analytical forms may not approximate the solution well, for example, the behavior of a region at the end of a narrow projection of a machine tool subjected to rapid, inelastic strain. 
     In these situations, it is desirable to refine the original decomposition in regions where primitives are too coarse. This is easily accomplished with simplicial primitives but not with other primitives because there are no algorithms that automatically create meshes from a boundary specification; this is still typically performed with human assistance. The refinement process must preserve two important properties of meshes to be useful for most simulation techniques. 
     One of the properties that refinement schemes must preserve is mesh conformality. A mesh of n-dimensional polytopes is conformal if every (n−1)-boundary (face for n=3) of a primitive is shared with at most one neighboring primitive. Furthermore, when an (n−1)-boundary is shared, it must be defined identically by both neighboring primitives; every point on the shared (n−1)-boundary of one primitive must appear on the corresponding shared (n−1)-boundary of the neighbor. Conformal meshes are called complexes. 
     Refinement schemes must also preserve—or at least gracefully degrade—the geometric quality of the primitives that compose the mesh. There are many different measures of geometric quality, but perhaps the minimum requirement is that: assuming all input primitives are convex then all output primitives shall be convex. Convexity is an important property in finite-element simulations because, for typical parameterizations of primitive polytopes, concavity implies that the Jacobian is negative inside the element. The geometric Jacobian appears in the quadrature formulation used to determine values of scalar, vector, and tensor fields and a negative Jacobian incorrectly weights the contribution of a region of the element to the whole, usually resulting in a divergent simulation. 
     Besides preserving properties of the mesh, any refinement scheme also requires a way for a simulation to indicate which regions require refinement. In 3-dimensional simplicial (i.e., tetrahedral) refinement schemes, this is typically accomplished by defining a target edge size at vertices (vertex-based refinement); by indicating a set of edges to be subdivided (edge-based refinement); or by indicating a set of faces to be subdivided (face-based refinement). We have shown (Pébay, P. P. and Thompson, D., “Communication-free streaming mesh refinement,”  Journal of Computing and Information Science in Engineering,  2005, v.5(4): pp. 309-316; and Thompson, D. and Pébay, P. P., “Embarrassingly parallel streaming mesh refinement,” Engineering and Computers, 2006, v.22(2): pp. 75-93, 13 th  International Meshing Roundtable) how tetrahedral edge-based refinement can be performed in parallel with little or no communication—depending on whether vertices shared between processes need to be cross-referenced. 
     RELATED ART 
     Nishigaki et al. (U.S. Pat. No. 6,578,189) discloses a hexahedral mesh generation method. However, their technique does not allow boundary specification, so it cannot be used to perform local mesh refinement. 
     Other related art in the literature involves template-based hexahedral refinement approaches where choices made for one input hexahedron must be propagated to arbitrarily many other hexahedra in the mesh since each of the neighboring hexahedra may need to alter its template in response to a change from a neighbor (Tautges, T. J. and Knoop, S. E., “Combinatorial aspects of dual-based hexahedral mesh modification,” SIAM Workshop on Combinatorial Scientific Computing, (2004); Zhang, Y. and Bajaj, C. “Adaptive and quality quadrilateral/hexahedral meshing from volumetric data,” Proceedings of the 13 th  International Meshing Roundtable (2004)). Similarly, there are mesh-generation approaches which require the propagation of constraints across regions (Staten, M., Owen, S., and Blacker, T., “Unconstrained paving and plastering: A new idea for all-hexahedral mesh generation,” Proceedings of the 14 th  International Meshing Roundtable, (2005), pp. 399-416; and Staten, M. L., Kerr, R. A., Owen, S. J., and Blacker, T. D., “Unconstrained paving and plastering, Progress update,” Proceedings of the 15 th  International Meshing Roundtable, (2006), pp. 469-486). 
     Finally, some have attempted hex-dominant mesh generation and refinement. In those techniques, hexahedra are preferred but other primitives (generally tetrahedra) are allowed when it is difficult or impossible to determine a configuration of hexahedra that matches a given geometric constraint (Tuchinsky, P. and Clark, B. W., “The hex-tet, hex-dominant automesher: An interim progress report,” Proceedings of the 6 th  International Meshing Roundtable (1997)). 
     Accordingly, what is desired, and what has not been developed heretofore, is a method for locally refining polytopic meshes according to some externally provided criterion in a way that preserves mesh conformality without requiring updates to neighboring primitives. 
     SUMMARY 
     It is an object of the present invention to provide a method for boundary-based refinement of a polytope that yields a conformal partition and preserves convexity. 
     It is an object of the present invention to provide a method for conformality-preserving, convexity-preserving, boundary-based refinement of non-simplicial meshes containing polytopic primitives. This is accomplished by a series of homothetic transformations of polytopes that isolate topological changes near the boundaries marked to be subdivided from the rest of the mesh. 
     Therefore, it is an object of this invention to provide a system where 3-dimensional, all-hexahedral, conformal meshes may be refined (remaining all-hexahedral) in parallel with little or no communication between a plurality of computer processes carrying out the refinement on some portion of the mesh. 
     It is also an object of this invention to provide a method for determining which boundaries of a polytope should be subdivided based on the value of an implicit function φ({right arrow over (x)}) and its gradient ∇φ({right arrow over (x)}) evaluated at a point or plurality of points on each boundary. 
     Both the foregoing general description and the following detailed description are exemplary and explanatory only and are intended to provide further explanation of the invention as claimed. 
    
    
     
       BRIEF DESCRIPTION OF THE DRAWINGS 
       The accompanying drawings are included to provide a further understanding of the invention and are incorporated in and constitute part of this specification, illustrate one or more embodiments of the invention and together with the description, serve to explain the principles of the invention. The drawings are only for the purpose of illustrating one or more preferred embodiments of the invention and are not to be construed as limiting the invention. In the drawings: 
         FIG. 1A  illustrates a perspective view of a single input hexahedron  1 . 
         FIG. 1B  illustrates a marked edge  2  indicating an edge to be subdivided on input hexahedron  1 . 
         FIG. 2A  shows a perspective view of a first configuration of homothecy showing an interior homothecy center  3  strictly interior to hexahedron  1 . 
         FIG. 2B  shows a perspective view of a second configuration of homothecy showing a face-bound homothecy center  5  strictly interior to a face of hexahedron  1 . 
         FIG. 2C  shows a perspective view of a third configuration of homothecy showing an edge-bound homothecy center  7  strictly interior to an edge of a hexahedron  1 . 
         FIGS. 3A and 3B  show a pair of perspective views of input hexahedron  1  with an interior homothecy applied. 
         FIG. 4  shows an orthographic view of an input hexahedron  1  subdivided into seven intermediate hexahedra  12 - 18 . 
         FIG. 5A , for each face of the input hexahedron  1  containing one or more marked edges  2  a homothecy with a face-bound center is applied to the corresponding hexahedron from  FIG. 3B , resulting in six hexahedra that replace one of the seven from  FIG. 3B . 
         FIG. 5B  shows an isolated view of the results of the face-bound homothecy applied to intermediate hexahedron (+X)  13 . 
         FIG. 5C  shows, for one of those four of the six resulting hexahedra that share an edge with the input hexahedron  1 , the result of a new edge-bound homothetic transform that is applied if that edge is marked. 
     
    
    
     DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS 
     Some applications of the finite element method—such as solid mechanics where shear of thin walls is involved—an provide more accurate results with an all-hexahedral mesh than with, for example, a tetrahedral mesh, by virtue of requiring fewer degrees of freedom. However, all-hexahedral meshing in general can be considered an open problem. Trivial approaches, such as using an automatic tetrahedral mesh generator and decomposing each output tetrahedron into four hexahedra, although automatic, yield meshes of such poor geometric quality that they cannot be used in practice. No practical, fully automatic hexahedral mesh generator is available to date; instead, existing generators require significant user intervention, or generate hybrid meshes, or fail with some inputs. 
     For these reasons, once a hexahedral mesh of the domain of interest has been obtained, it is important that the best computational use be made of it. In particular, finite element techniques known as h-adaptive methods have been developed in an effort to converge to a solution more rapidly and/or with better accuracy than previously possible. Finite element solvers that incorporate h-refinement have become increasingly popular, and it is therefore of interest to make this technique available for hexahedral meshes. Several approaches have been proposed; however, current techniques either yield non-conformal or hybrid meshes, or these methods require communication of refinement patterns to neighboring elements and iteration until all elements have a pattern that can be accommodated. Furthermore, as problem sizes grow, these tasks must be performed in parallel—which adds a penalty to methods requiring communication of boundary information between portions of a model residing in different processes. 
     To achieve parallel scalability, what is needed is a local, conformal refinement scheme, similar to that which we had developed previously for tetrahedral mesh refinement (Pébay P. P. and Thompson D. “Communication-free streaming mesh refinement,”  Journal of Computing and Information Science in Engineering,  2005, v.5(4): pp. 309-316). The keystone of this approach is to find hexahedral subdivision templates, given a set of edges to be refined (e.g., based on a target size from a posteriori error estimation), such that neighboring hexahedra always generate the same quadrilaterals on their shared boundary. Such templates ensure that no neighbor information is required for the refined mesh to be conformal everywhere. This global conformality obtained by means of a local criterion allows the method to be implemented in a streaming fashion, and in particular to achieve optimal parallel speedup under load-balanced conditions. 
     We have developed such a scheme using homothetic transformations to partition a hexahedron into subhexahedra according to some indicator function. The following detailed description will introduce the mathematical basis for the approach (e.g., using homothetic transformations, complexes, and indicator functions) and then demonstrate a simple scheme for local refinement that always produces convex subhexahedra given a convex input. 
     For the sake of generality, we will describe the application of the scheme based on polytopes in an arbitrary dimension. In particular, we have developed a method for partitioning a convex, non-simplicial, n-dimensional input polytope into a conformal set of output polytopes based on an indicator function. The indicator function takes values on over m-boundaries of the input polytope where 1≦m≦n−1. 
     Our main insight is that a series of homothetic transforms may be applied to an original element such as an input polytope, such that the (n−1)-boundaries of the input polytope, the (n−1)-boundaries of the transformed polytope, and (n−1)-planes joining pairs of (n−2)-boundaries (that are boundaries of the (n−1)-boundaries above) partition the input polytope into a set of convex output polytopes. For the case where n=3 and the input polytope is a hexahedron isomorphic to a 3-dimensional cube, all of the output polytopes are hexahedra isomorphic to a 3-dimensional cube as shown for example in  FIGS. 3 ,  4 , and  5 . 
     As part of the development of such a method, knowledge of indicator functions, implicit functions, and homothetic transforms will be required as indicated below. 
     An indicator function (or characteristic function) is an injective map χ from some set of interest S to the set {0, 1}. In our application, the domain of this map will be the set of m-boundaries (mε{(n−1), (n−2), . . . , 1}) of an n-dimensional input polytope. Elements of S which χ maps to 1 will be considered “marked” for subdivision while elements of S which χ maps to 0 must remain whole as the input polytope is partitioned into output primitives. 
     An implicit function is a function of the form ƒ({right arrow over (x)})=0. Where explicit functions provide a parametric enumeration of the members some manifold (i.e., an expression {right arrow over (x)}({right arrow over (t)}) which lists all {right arrow over (x)} for which ƒ({right arrow over (x)})=0 with {right arrow over (t)} in some simple domain), implicit functions must be inverted to identify members of the manifold where the equation holds. Typically implicit functions are smooth and we will require piecewise smooth functions for our application. Implicit functions, and specifically those which estimate the geometric signed distance to some manifold, are commonly used in computer graphics as indicator functions that identify the boundary of a surface to be rendered and are therefore well-known in the art. Our application will use implicit indicator functions to identify boundaries of polytopes where refinement should occur. If some point {right arrow over (x)} on an m-boundary (m&lt;n) of an input polytope satisfies |ƒ({right arrow over (x)})|&lt;ε, the boundary will be considered “marked” for subdivision. Here ε&gt;0 is a number chosen to normalize ƒ. The point {right arrow over (x)} may be defined in several ways. Examples include but are not limited to the case where {right arrow over (x)} is a reference point, such as a mid-point; where {right arrow over (x)} is chosen such that some property, such as ƒ, varying over the boundary is maximized or minimized; or where {right arrow over (x)} is dependent on the homothetic center {right arrow over (O)} defined below. An indicator function of the form: 
                     χ   ⁡     (   s   )       =     {         1                  f   ⁡     (       x   →     ⁡     (   s   )       )              ɛ   ⁡     (   s   )         &lt;   1             0       otherwise                   (   1   )               
may then be used to identify polytope boundaries that are near the surface ƒ({right arrow over (x)})=0. Note that where ƒ is smooth, ∇ƒ({right arrow over (x)}) is defined and is a vector normal to the surface. The indicator function may be modified to favor polytope boundaries that have some preferred orientation relative to this surface.
 
     A homothetic transform or “homothecy” is an affine transform that can be expressed solely as a translation and uniform scaling—no rotation or shear is allowed. Thus, an edge or face is parallel to its image under a homothetic transform. Another property of homothecies is that they have a unique fixed point. This fixed point (the origin {right arrow over (O)} about which scaling occurs) is known as the homothetic center. Mathematically, a homothecy can be written:
 
 h ({right arrow over ( x )})= {right arrow over (O)}+c ( {right arrow over (x)}−{right arrow over (O)} ), and  (2)
 
 h:P→T,   (3)
 
where {right arrow over (O)} is the homothetic center, {right arrow over (x)} is the untransformed point, h({right arrow over (x)}) is the transformed point, c is the contraction or dilation factor, P is the set of all points in some convex polytope, and T is the range of the map h. In those cases where c is positive and less than 1 the transformation moves points closer together while the opposite is true for those cases where c is greater than 1.
 
     When a homothetic transform with 0&lt;c&lt;1 is applied to all the points on the boundary ∂P of a polytope P, a smaller, transformed version ∂T of that polytopic boundary is produced. Denote P without boundary as P*=P\∂P and the closure of P as  P =P∪U∂P. When the homothetic center of the transform is strictly inside a polytope (i.e., {right arrow over (O)}εP*), we will call the transform an interior homothecy. When the homothetic center lies on a k-boundary within ∂P but is not in the boundary of this k-boundary (i.e., not on a (k−1)-boundary), we call the transform a k-bound homothecy. More specific terminology for the 3-dimensional case is as follows. When the homothetic center lies on a boundary face B of a polyhedron but is not in the boundary of the face (i.e., not on an edge of the face), we call the transform a face-bound homothecy. And when the homothetic center lies on an edge of ∂P but is not in the boundary of the edge (i.e., not on a corner vertex), we call the transform an edge-bound homothecy. Note that whenever 0&lt;c&lt;1 and where {right arrow over (O)}ε  P , then T⊂  P . Furthermore, when {right arrow over (O)}εB, then h(B)⊂∂P; that is, a face-bound homothecy will have a transformed element with one of its faces coincident to the corresponding face of polytope P. Similarly, an edge-bound homothecy will have a transformed element with two of its faces and their shared edge coincident to the two corresponding faces and shared edge of polytope P. 
     The generalized method then for an arbitrary n is: 
     If no boundaries  2  are selected by the indicator function then the output polytope is not transformed and it remains identical to the input. 
     If, however, boundaries  2  are selected by the indicator function, then a homothetic transform is applied with a center interior to the input polytope and corresponding (n−2)-boundaries of the input polytope and transformed polytope are joined with corresponding (n−1)-planes such that a level-0 partition is constructed; 
     for any member of the level-0 partition that shares an m-boundary with the input polytope that is selected by the indicator function apply (n−1)-bound homothetic transforms to obtain a level-1 partition; and, 
     for any member of the level-(l−1) partition that shares an m-boundary with the input polytope that is selected by the indicator function, for l =2, 3, . . . , n−m, apply (n−l )-bound homothetic transforms to obtain a level-l partition. 
     The output polyhedra are then all members of the level-(n−m) partition and any members of the previous partitions not covered by members of the level-(n−m) partition. 
     The primary application of the present method is the case where n=3, the input polyhedron is a hexahedron isomorphic to a 3-dimensional cube (see  FIGS. 1A-B ), and the indicator function is defined on 1-boundaries (edges).  FIG. 1  shows a l-boundary marked for subdivision. 
     As noted in the foregoing sections, an implicit function can be used in an indicator function defined on edges and the zero-set of this function will then specify where refinement should occur. Scaling this implicit function by the size of the m-boundary (the edge length in our primary application) will force long edges to be refined when they are further from this zero-set than short edges. (Example measures of the size of an m-boundary are the diameter, the radius of gyration, and the mth root of the m-dimensional volume.) Modifying the function so that m-boundaries with particular orientations are favored is another strategy of interest. As an example, consider in the present case: modulating the implicit function by the dot product of the gradient of the implicit function and the edge direction will force edges pointing directly towards or away from the zero-set to be subdivided while those edges which are parallel to the zero-set will not be subdivided. 
     Finally, the exact location of the zero off along an edge to be subdivided can be used to locate the center of a homothetic transform or the vertex of a transformed polyhedron. Thus, if we have an implicit function that describes a surface of interest, a volumetric mesh of that surface can be created by refining an input polyhedron containing the surface and discarding any output polyhedra that fall outside the surface. Thus refinement can also be used for mesh generation. 
     Once refined, the input mesh may be inspected or saved for future use but is most frequently used by one of any number of simulators to predict how the physical domain represented by the generated mesh will respond/react to hypothetical events or conditions present in the environment. Over time, simulators have used various means for providing the sought after predictions. In the past, simulators were actual scale models subjected to scaled conditions expected in the environment of a full scale structure. More recently, simulators were composed of analog electrical circuits whose properties mimicked the system being simulated. Most pertinent to the present application is the case where the simulator is a finite element or finite volume method embodied as instructions in a computer process—perhaps even the same process that performs the mesh refinement—which will utilize the generated mesh to predict and present results either graphically as a printed report or visually as an electronic file displayed on a computer monitor or similar device. 
     EXAMPLE 
     For the case where n=3 and input polytopes are hexahedra isomorphic to a 3-dimensional cube,  FIGS. 1-5  illustrate the steps of our method.  FIG. 1A  shows a perspective view of such an input hexahedron  1 .  FIG. 1B  shows marked edge  2 , which specifies that input hexahedron  1  is to be subdivided. Each of the three types of homothecy illustrated in  FIGS. 2A-2C  is used on different hexahedra created during the refinement process. The interior lines connecting adjacent vertices indicate how the vertices of input hexahedron  1 , and the transformed hexahedron  4 , transformed hexahedron  6 , or transformed hexahedron  8  (respectively), can be used to create a decomposition of a hexahedron into seven, six, or five hexahedra (respectively). 
     When at least one edge is marked for subdivision, an interior homothecy center  3  is used with 0&lt;c&lt;1 to create a transformed hexahedron  4  (shown in  FIG. 2A ), which is completely interior to the input hexahedron  1 . Each corresponding input hexahedron face  9  (shown in  FIG. 3A ) and transformed hexahedron face  10  are joined by a joining face  11  (shown in  FIG. 3B ) on each of their four edges to create a set of six hexahedra in addition to the transformed hexahedron  4 . An exploded orthographic view of the result of this process is shown in  FIG. 4  which shows the seven subdivided intermediate hexahedra. 
     Note that in  FIGS. 3A and 3B  input hexahedron face  9  and transformed hexahedron face  10  are parallel to one another by virtue of the homothetic transform applied to obtain transformed hexahedron  4 . This property guarantees that each of the corresponding edge pairs bounding the faces are contained by a joining face  11 . The faces of the input hexahedron  1 , the faces of the transformed hexahedron  4 , and joining faces  11 , therefore, divide the input hexahedron  1  into seven intermediate hexahedra. Because the seven hexahedra may each be formed as the intersection of six planar half-spaces, each will be convex. 
       FIGS. 5A-5C  show a series of three perspective views continuing the illustration of the refinement process begun in  FIGS. 3A and 3B . For any of these six hexahedra whose input face  9  contains at least one edge marked for further subdivision, this hexahedron will have a homothecy applied to it with a face-bound homothecy center  5  (shown in  FIG. 2B ), such that a face of the resulting transformed hexahedron  6  (shown in  FIG. 5A ) is coplanar with input face  9 . A face-bound homothecy results in a single hexahedron being divided into six intermediate hexahedra that replace one of the seven hexahedra shown in  FIG. 3B , that is, the hexahedron defined by faces  9 ,  10 , and  11  is divided into the transformed hexahedron  6 , the four “surrounding” hexahedra, and a “center” hexahedron behind the transformed hexahedron  6 , shown in  FIG. 5B : an isolated view of the results of the face-bound homothecy applied to intermediate hexahedron (+X)  13 . 
     Lastly, in the case where an edge of any of the four “surrounding” hexahedra shown in  FIG. 5B  that share an edge with input hexahedron  1 , is marked for further subdivision, then a final homothecy is applied. This time a homothecy with the edge-bound homothecy center  7  (shown in  FIG. 2C ) is applied so that an edge of the transformed hexahedron  8  (shown in  FIG. 5C ) lies on an edge of the original input hexahedron  1 . In this manner, any combination of edges to be subdivided may be handled.  FIG. 5C  shows, for one of those four of the six resulting hexahedra from  FIG. 5B  that share an edge with the input hexahedron  1 , the result of a new edge-bound homothetic transform that is applied if that edge is marked. For the input hexahedron  1 , with marked edge  2  (shown in  FIG. 1B ), intermediate hexahedron (+X)  13  and intermediate hexahedron (−Y)  14  would each have one face-bound homothecy and one edge-bound homothecy applied. Only intermediate hexahedron (+X)  13  is shown. 
     Note that none of the subdivisions performed specify the exact homothecy center and contraction factor used to generate the output primitives. While the homothecy center is constrained to lie in the interior of input hexahedron  1 , on a face of input hexahedron  1 , or on an edge of the input hexahedron  1 , the exact position of the homothecy center and the contraction factor may be chosen to maximize some geometric quality of the output hexahedra, provided the homothecy center and the contraction factor for faces and edges are chosen identically for all input hexahedra sharing the face or edge. 
     Having thus described exemplary embodiments of the present invention, it should be noted by those skilled in the art that the disclosures herein are exemplary only, and that various other alternatives, adaptations, and modifications may be made within the scope of the present invention. Accordingly, the present invention is not limited to the specific embodiments as illustrated herein, but is only limited by the following claims. 
     Finally, to the extent necessary to understand or complete the disclosure of the present invention, all publications, patents, and patent applications mentioned herein are expressly incorporated by reference therein to the same extent as though each were individually so incorporated.