Analysis-aware higher order smooth three-dimensional representations: creation, simulation and visualization

Volumetric parameterization is an emerging field in computer graphics, where volumetric representations that have a semi-regular tensor-product structure are desired in applications such as three-dimensional (3D) texture mapping and physically-based simulation. At the same time, volumetric parameterization is also needed in the Isogeometric Analysis (IA) paradigm, which uses the same parametric space for representing geometry, simulation attributes and solutions. One of the main advantages of the IA framework is that the user gets feedback directly as attributes of the NURBS model representation, which can represent geometry exactly, avoiding both the need to generate a finite element mesh and the need to reverse engineer the simulation results from the finite element mesh back into the model. Research in this area has largely been concerned with issues of the quality of the analysis and simulation results assuming the existence of a high quality volumetric NURBS model that is appropriate for simulation. However, there are currently no generally applicable approaches to generating such a model or visualizing the higher order smooth isosurfaces of the simulation attributes, either as a part of current Computer Aided Design or Reverse Engineering systems and methodologies. Furthermore, even though the mesh generation pipeline is circumvented in the concept of IA, the quality of the model still significantly influences the analysis result. This work presents a pipeline to create, analyze and visualize NURBS geometries. Based on the concept of analysis-aware modeling, this work focusses in particular on methodologies to decompose a volumetric domain into simpler pieces based on appropriate midstructures by respecting other relevant interior material attributes. The domain is decomposed such that a tensor-product style parameterization can be established on the subvolumes, where the parameterization matches along subvolume boundaries. The volumetric parameterization is optimized using gradient-based nonlinear optimization algorithms and datafitting methods are introduced to fit trivariate B-splines to the parameterized subvolumes with guaranteed order of accuracy. Then, a visualization method is proposed allowing to directly inspect isosurfaces of attributes, such as the results of analysis, embedded in the NURBS geometry. Finally, the various methodologies proposed in this work are demonstrated on complex representations arising in practice and research.

ANALYSIS-AWARE HIGHER ORDER SMOOTH THREE-DIMENSIONAL REPRESENTATIONS: CREATION, SIMULATION AND VISUALIZATION by Tobias Martin A dissertation submitted to the faculty of The University of Utah in partial fulfillment of the requirements for the degree of Doctor of Philosophy in Computer Science School of Computing The University of Utah May 2012 Copyright c Tobias Martin 2012 All Rights Reserved The University of Utah Graduate School STATEMENT OF DISSERTATION APPROVAL The dissertation of has been approved by the following supervisory committee members: , Chair Date Approved , Member Date Approved , Member Date Approved , Member Date Approved , Member Date Approved and by , Chair of the Department of and by Charles A. Wight, Dean of The Graduate School. Tobias Martin Elaine Cohen 12/9/11 Robert M. Kirby 12/9/11 Richard F. Riesenfeld 12/9/11 Charles Hansen 12/9/11 Tom Lyche 12/9/11 Alan L. Davis School of Computing ABSTRACT Volumetric parameterization is an emerging field in computer graphics, where vol-umetric representations that have a semi-regular tensor-product structure are desired in applications such as three-dimensional (3D) texture mapping and physically-based simulation. At the same time, volumetric parameterization is also needed in the Isoge-ometric Analysis (IA) paradigm, which uses the same parametric space for representing geometry, simulation attributes and solutions. One of the main advantages of the IA framework is that the user gets feedback directly as attributes of the NURBS model representation, which can represent geometry exactly, avoiding both the need to generate a finite element mesh and the need to reverse engineer the simulation results from the finite element mesh back into the model. Research in this area has largely been concerned with issues of the quality of the analysis and simulation results assuming the existence of a high quality volumetric NURBS model that is appropriate for simulation. However, there are currently no generally applicable approaches to generating such a model or visualizing the higher order smooth isosurfaces of the simulation attributes, either as a part of current Computer Aided Design or Reverse Engineering systems and methodologies. Furthermore, even though the mesh generation pipeline is circumvented in the concept of IA, the quality of the model still significantly influences the analysis result. This work presents a pipeline to create, analyze and visualize NURBS geometries. Based on the concept of analysis-aware modeling, this work focusses in particular on methodologies to decompose a volumetric domain into simpler pieces based on appropriate midstructures by respecting other relevant inte-rior material attributes. The domain is decomposed such that a tensor-product style parameterization can be established on the subvolumes, where the parameterization matches along subvolume boundaries. The volumetric parameterization is optimized using gradient-based nonlinear optimization algorithms and datafitting methods are introduced to fit trivariate B-splines to the parameterized subvolumes with guaranteed order of accuracy. Then, a visualization method is proposed allowing to directly inspect isosurfaces of attributes, such as the results of analysis, embedded in the NURBS geometry. Finally, the various methodologies proposed in this work are demonstrated on complex representations arising in practice and research. iv To Sumathi and my parents CONTENTS ABSTRACT . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iii LIST OF FIGURES. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ix LIST OF TABLES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xix ACKNOWLEDGEMENTS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xx CHAPTERS 1. INTRODUCTION . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1 Outline. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.2 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 2. RELATED WORK . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 2.1 Model Completion, Parameterization and Meshing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 2.2 Model and Mesh Quality . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 2.3 Isosurface Visualization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 3. ANALYSIS-AWARE MODELING . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 3.1 Nomenclature. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 3.2 Outline. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 3.3 The Modeling to Analysis Pipeline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 3.4 Mathematical Formulation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 3.5 Parametric Representation of Geometry . . . . . . . . . . . . . . . . . . . . . . . . . 50 3.6 Studies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63 3.7 Summary and Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89 4. CYLINDRICAL-LIKE OBJECTS WITH SHAPE OVERHANGS . 90 4.1 Preliminaries and Notation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92 4.2 Framework Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95 4.3 Modeling the Exterior . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96 4.4 Modeling the Interior . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103 4.5 Tracing w-paths. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105 4.6 Possible Extensions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110 4.7 Experiments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111 4.8 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113 5. OBJECTS WITH HIGHER GENUS AND BIFURCATIONS . . . . . 114 5.1 Parameterization Strategies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118 5.2 Framework Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122 5.3 Midsurface Construction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123 5.4 Decomposition of Volume . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128 5.5 Volumetric Parameterization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133 5.6 Modeling Examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 135 5.7 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 137 6. MORE GENERAL OBJECTS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 140 6.1 Volumetric Representation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 143 6.2 Tetrahedralize the Interior . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 149 6.3 Simulation Studies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 153 6.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 157 6.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 159 7. GENERALIZED SWEPT MIDSTRUCTURE . . . . . . . . . . . . . . . . . . . 161 7.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 161 7.2 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 166 7.3 The GSM. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 167 7.4 Harmonic Functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 168 7.5 Computing Midstructure of Slices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 170 7.6 Matching Successive 2D Midstructure . . . . . . . . . . . . . . . . . . . . . . . . . . . 172 7.7 Results and Discussion. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 177 7.8 Applications . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 181 7.9 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 182 8. NONLINEAR OPTIMIZATION . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 184 8.1 Contributions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 185 8.2 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 187 8.3 Sobolev Preconditioning Background. . . . . . . . . . . . . . . . . . . . . . . . . . . . 187 8.4 Mathematical Framework . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 188 8.5 Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 193 8.6 Applications . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 197 8.7 Conclusions and Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 208 9. DATA FITTING . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 210 9.1 Preliminaries . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 213 9.2 Convergence . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 214 9.3 Simplification . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 217 10. DIRECT ISOSURFACE VISUALIZATION . . . . . . . . . . . . . . . . . . . . . 219 10.1Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 221 10.2Mathematical Formulation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 228 10.3Ray Frustum/Isosurface Intersection . . . . . . . . . . . . . . . . . . . . . . . . . . . . 233 vii 10.4Determining the Set of Intersection Patches . . . . . . . . . . . . . . . . . . . . . . 244 10.5Analysis and Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 245 10.6Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 250 11. CONCLUSION . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 251 11.1Contributions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 251 11.2Outlook . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 253 APPENDIX: PUBLICATIONS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 255 REFERENCES. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 257 viii LIST OF FIGURES 1.1 Input data. On the left: Triangle mesh boundaries of a femur bone; on the right: Triangle mesh boundaries of a pelvis bone. In both examples, the inner boundary separates cortical from trabecular bone material. . . . . 3 1.2 Overview of the pipeline introduced in this work. In the first step the volume of an input model is completed. Then, this parameterization is optimized and a higher order representation is generated from it. Simu-lation is applied to the volumetric shapre respresentation which is then visualized and inspected using direct isosurface visualization. . . . . . . . . . . 6 1.3 Sweeping can cause significant of parametric distortion, if the object is irregular. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 1.4 In (a), a harmonic function is computed through the volume of a configu-ration space. In (b), based on the gradient field of the harmonic function in (a), robot paths from start to corresponding end location are computed. 10 1.5 Orthogonal vector fields on a square domain. . . . . . . . . . . . . . . . . . . . . . . . 13 1.6 Parameterization on domains in 2D. On the left: Since the domain consists of right angle corners and straight lines, two orthogonal vector fields, resulting in a parameterization without any parametric distortion can be established. On the right: No orthogonal vector fields can be established, because the boundary of the domain is curved. . . . . . . . . . . . . . . . . . . . . . 13 2.1 Examples for mesh quality metrics in traditional finite elements. . . . . . . . 21 2.2 Discretization of domain from Figure 10.1c with 300k tetrahedra and application of marching tetrahedra (using ParaView). (a) isosurface; (b) scalar field on tetrahedra; (c) our approach on single triquintic NURBS patch . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 2.3 Isosurface from silicium data set (volvis.org), isovalue of 130 using March-ing Cubes (using ParaView), Afront (ρ = 0.3) and Direct visualization with our proposed method. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 3.1 F maps a point from the reference domain Θ to the physical domain Ω. Correspondingly, F−1 maps a point from Ω back to Θ. With F−1, it is possible to define basis functions on Ω as compositions of F−1 with basis function defined on Θ. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 3.2 Diagram presenting the classic concept-to-mesh pipeline (top branch) and the concept-to-model pipeline (bottom branch). A detailed discussion of the diagram is presented in the text. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 3.3 Boundary curves. (a) in 2D; (b) in 3D. . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 3.4 Coefficients of the geometry maps. The top image portrays the evenly spaced coefficients of Uk(x), k = 1, 2. The second and third rows are coefficients of Ik(x), the identity, for k = 1, 2. They depend on the knot vectors. The fourth row shows the coefficients of M1 are mostly the same as the coefficients of U1, except for the double control point. . . . . . . . . . . . 55 3.5 Filled boundary curves. (a) in 2D; (b) in 3D . . . . . . . . . . . . . . . . . . . . . . . 56 3.6 Three ways to represent a disc of radius one. Column (a) is mapped with D1, column (b) with D2, and column (c) with D3. Circles around black control points mean that more than one control point sits on the black control point location. Degeneracies are marked with red points. These representations are used to solve the drumhead problem. . . . . . . . . . . . . . 59 3.7 Input (left) and output (right) of the modeling methodology proposed in Chapter 4. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63 3.8 Given the knot vectors τ1 and τ2, the first column shows ψi,d,τk (U−1 k (x)), and the second shows ψi,d,τk (I−1 k (x)) (k = 1, 2). Note that ψi,d,τk(I−1 k (x)) are not stretched, but setting a and b appropriately causes the end func-tions to exhibit wider support and more resemble ψi,d,τk (Uk), k = 1,2 . . . . 66 3.9 Approximated eigenfunctions (exact solution in gray) for the modes k = 3, 7 using the nonuniform mapping U1 (first row), and the identity map-pings I1 (second row), I2 (third row). . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 3.10 Approximated eigenfunctions (exact solution in gray) for the modes k = 13, 18 using the nonuniform mapping U1 (first row), and the identity mappings I1 (second row), I2 (third row). . . . . . . . . . . . . . . . . . . . . . . . . . 68 3.11 Two cubic identity mappings constructed from different knot vectors. The particular choice of nonuniform knot vector yields a flatter normalized discrete spectra for the vibrating rod. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69 3.12 The cubic basis functions defined on Ω = [0, 1] for VM1,S1 and VU1,S1, respectively. The control points for the M1 are uniformly spaced, except the control point in the middle is duplicated. ψi,d,τ1 (U−1 1 (x)) are shown for comparison (bottom). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71 3.13 h-refinement: For both mappings (degenerate and nondegenerate), the normalized discrete spectra are slightly worse under uniform refinement. The uniformly refined space generated by M1 always has outliers. . . . . . . . 72 3.14 k-refinement: For both mappings (degenerate and nondegenerate), the normalized discrete spectra become flatter with elevated degree. However, the outlier exists throughout uniformly refining for the spaces based on M1. 73 3.15 Degenerate mapping (red) with control point 1/2 duplicated and corre-sponding inverse mapping (darker red). The inverse has a high slope and therefore large first derivative at t = 1/2 causing a negative impact on the stiffness matrix conditioning (results in larger eigenvalues). . . . . . . . . . . . . 74 x 3.16 A triple control point embeds a horizontal tangent at t = 1/2 for M2. Hence M−1 2 is singular at t = 1/2. However t = 1/2 is a knot so the singularity occurs only at an element boundary. The comparison with U1 (green) where slope of the inverse (dark green) is constant away from the boundary elements. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75 3.17 Normalized discrete spectra for different representations of a disc. . . . . . . 76 3.18 The disc model D3,2 with a nonuniform knot vector, where refining knots are chosen as to create τ2 ( Section 3.5.2). It yields a flatter discrete spectrum than the normalized spectrum showed in Figure 3.17. . . . . . . . . 77 3.19 Both NURBS models represent Ω exactly. S2 on the left is the identity map; the control points in S1 on the right are uniformly spaced. . . . . . . . . 78 3.20 Three exact representations of the boundary of a quarter annulus with an inner radius of one and an outer radius of two with different completions. In (b), control points in the interior are slightly perturbed. In (c), element boundaries are orthogonal where they cross. . . . . . . . . . . . . . . . . . . . . . . . 79 3.21 Poisson test 1 and 2. The plots on the left show the error curves for different mappings, on which the exact solution on the respective domain is shown on the right. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82 3.22 Poisson test 3 and 4. The plots on the left show the error curves for different mappings, on which the exact solution on the respective domain is shown on the right. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83 3.23 Poisson test 5. The plots on the left show the error curves for different mappings, on which the exact solution on the respective domain is shown on the right. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84 3.24 Linear elastic deformation of a femur. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88 4.1 Bimba input. (a) triangle mesh of Bimba statue; (b) corresponding trivari-ate B-spline, where interior represents material information used in simu-lation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91 4.2 A femur consists of two materials: The outer solid part, or cortical bone, represented by the volume between the input triangle meshes; and the inner soft part, or trabecular bone, represented by the volume of interior triangle mesh (red). These volumes are parameterized (middle) and a sin-gle trivariate tensor product B-spline is fitted against it (right), respecting the input triangle meshes in its parameterization. This makes it easier to specify respective material properties. Black isolines represent knotlines in the trivariate parameterization. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93 4.3 Critical paths end at the edge of a triangle, where one of its vertices is νmin or νmax. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97 4.4 X2 (left) is not bijective. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98 xi 4.5 Establishing the surface parameterization X2. On the left: the harmonic function uTB defined by two critical points is established over TB; middle: Based on uTB, the orthogonal harmonic function vTB is calculated. At this stage uTB and vTB define an injective transformation; on the right: Scaling and translation yields the parameterization X2. . . . . . . . . . . . . . . . 99 4.6 X maps a vertical line at u0 in parameter space onto a closed isopara-metric line on TB. Accordingly, X maps a horizontal line at v0 onto an isoparametric which starts at νmin and ends at νmax. . . . . . . . . . . . . . . . . . 100 4.7 Either there is one vertex in the rectangle defined by the points p0, p1, p2 and p3, or none. Crosses mark edge intersections. . . . . . . . . . . . . . . . . . . . 102 4.8 A cross section of an object with an exterior boundary and an interior isosurface representing geometry or attribute data. The skeleton and boundaries were used to establish ∇w. Isolines visualize the uw-scalar field used to trace w-paths from the exterior to the skeleton. . . . . . . . . . . 105 4.9 Due to the linear property of ∇u and ∇v and special cases during w-path extraction as discussed in section 4.5, adjacent paths may collapse. . . . . . 109 4.10 w-parameterizations using min/max points. . . . . . . . . . . . . . . . . . . . . . . . 110 4.11 w-parameterizations using min/max paths. . . . . . . . . . . . . . . . . . . . . . . . . 111 4.12 Isogeometric analysis: Elastostatics applied to the data reduced trivariate B-spline representation of the femur. Load is applied to the head of the femur. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112 5.1 Parameterization of a disc. Top left: Polar-like with one degenerate point (magenta); Top right: Polycube-like with four corner points (blue); Bottom: Blend of the two parameterizations. . . . . . . . . . . . . . . . . . . . . . . 117 5.2 Parameterization of a genus-3 torus. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118 5.3 Parameterization of a CAD propeller of genus-1 with a single midsurface (midsurface boundary in yellow). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119 5.4 2D input consisting of an exterior domain boundary. The inner boundary, which is the interface between the polar and square parameterization, is specified by the user. Boundaries are decomposed to divide the domain into subregions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120 5.5 Parameterization strategies on 2D input. . . . . . . . . . . . . . . . . . . . . . . . . . . 121 5.6 Simplified medial axis' of Olivier hand data set (left: tight cocone [46], right: our approach). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123 5.7 Olivier hand model (n = 1): midsurface M (yellow and blue area) and base surface (blue area). Since n = 1, M =M trim. . . . . . . . . . . . . . . . . . 125 5.8 Polycube like parameterization of midsurface M on hand model. (a) parameterization; (b) corresponding parameteric domain . . . . . . . . . . . . . 127 xii 5.9 A trim value m0 closer to 1 (right) can introduce parametric distortion, but favors the periodic parameterization. A smaller trim values favors the polycube parameterization. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128 5.10 Parameterization of midsurface M of kitten model. (a) Partial view of parameterized midsurface M ; (b) Corresponding parameter domain; (c) Scalar field f is used to trim M . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 129 5.11 Construction of P+ = {p+1 , . . . , p+n } and corresponding P− = {p+1 , . . . , p+n }. p+i and p− i , respectively, are evaluations of its corresponding γi(t) and not shown in this figure. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 130 5.12 Construction of S+ and S−. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 131 5.13 Decomposition of Olivier hand data set. . . . . . . . . . . . . . . . . . . . . . . . . . . 133 5.14 Second volumetric parameterization strategy for Olivier hand data set. (Parametric cut in u, showing a part of the interior material boundary.) . . 136 5.15 Different choices of corners (red), resulting in two different scalar fields on a midsurface. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 137 5.16 Parameterized pelvis. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 138 6.1 Mixed element representation: Tricubic NURBS elements (red) at the boundary. Linear tetrahedra (gray) in the interior of the object. The surface patches are constructed from a Catmull-Clark subdivision object. . 142 6.2 Effect of C(2) elements at domain boundary: The configuration of this 2D elasticity example, based on a Hookean material using Chauchy-Green strain, is illustrated on the left. The top row shows the state of maxi-mal impact, i.e., the state where the elastic potential energy reaches its maximum, when the shape hits the ground plane. The bottom row shows the final state of the object. By increasing the number of boundary layers and its thicknesses, the corresponding final states do not exhibit flipped elements. Note, due to the gravitational force the final shape is not circular.144 6.3 The input are NURBS surface patches, in this case generated from a Catmull-Clark subdivision model. A sampled midstructure such as a simplified medial axis is used to generate trivariate NURBS patches at the boundary of the volumetric representation. The remaining interior is filled up with unstructured tetrahedral elements. (Elements in (b), (c) and (d) are culled to visualize the interior.) . . . . . . . . . . . . . . . . . . . . . . . . 146 6.4 Paths following ∇w (gray) emanate at cl (green points) and terminate on isosurface w(x, y, z) = w0 which is approximated with S (red). The choice of the midstructure allows creation of elements in thinner tubular regions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 149 6.5 Illustration of C(2) NURBS surface and corresponding C(0) triangle mesh. c l (blue) for the NURBS surface and xl (green) for the triangle mesh boundary, represent the same degree of freedom. Due to geometric con-cavities, the representations are overlapping. . . . . . . . . . . . . . . . . . . . . . . . 150 xiii 6.6 Mixed element representation in 2D. Quadratic NURBS elements at the boundary and linear triangles in the interior. Quadratic NURBS elements were chosen at the boundary to demonstrate the worst-case scenario of the alignment with triangles in the interior. As discussed in the text, convergence is achieved even in this more general scenario. Using cubic NURBS elements, the triangles are more closely aligned with the B-spline element boundaries (see Figure 6.2). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 151 6.7 Convergence behavior of hybrid representation. Left: convergence plot comparing the mixed element representation with the representation which only uses triangles. For the error computation, we chose the l2 norm between the approximated and true solution. n is the degree of freedom for the refinement steps. Solutions along z-axis on mixed element rep-resentation (middle) and triangle only representation (right), during one h-refinement step. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 155 6.8 Methodology applied to four objects. Elements are culled to better vi-sualize the filled interior of the objects. Object boundaries and interior boundaries are shown in red and green, respectively. The input for (a), (b) and (c) are bicubic B´ezier patches. The input for (d) are bicubic NURBS surfaces. The wire-frame shows C(0) continuities. . . . . . . . . . . . . . . . . . . . 158 7.1 Sheet based simplification using GSM (left) vs. global simplification using DSA (right). Magnified regions compare triangulations in a region where surfaces meet. Only the GSM shows 3 distinct sheets (red, blue, yellow) while the DSA shows similar polygons. The yellow GSM sheet (top left) is removed (mid left GSM). To simplify DSA in this region requires removal of additional structures in other regions (mid right DSA). Bottom row illustrates removal of a large blue sheet on one side of the GSM, whereas simplification of the DSA results in removal of polygons from both sides. . 162 7.2 The tight co-cone is intrinsic to the object, consisting of accurate but complex medial topology. A GSM identifies sheets (blue, red and white) and because of its swept generation, it has a natural parameterization. . . . 163 7.3 Pipeline for the Generalized Swept Midstructure (GSM): (a) User choice of harmonic function. (b) Object is decomposed into (curved) slices, and a medial axis is computed for each slice. (c) Slices are iteratively matched into a GSM (colored surface sheets). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 164 7.4 Two different harmonic functions on kitten model result in different GSMs.169 7.5 Medial axis of a region bounded by a B-spline curve. Along with their maximal circles, end points (E) are shown in orange, junction points (J) in blue, and distance critical points (D) in red. The arrows point to the foot points. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 172 7.6 A triangulation is used to connect these matching pairs based on the samples along the original midstructures. . . . . . . . . . . . . . . . . . . . . . . . . . . 173 7.7 The discrete cutting will likely not capture the degenerate point. . . . . . . . 174 xiv 7.8 Illustration of matching algorithm. The top row shows two consecutive slices, Li and Li+1 and their midstructures, Mi and Mi+1. The foot points of the end points (orange dots) and junction points (blue dots) are highlighted on the boundaries. Each point in the midstructure and its foot points are linked through straight lines. The bottom figures illustrates the matching steps 1-4. Note that all the foot points in level Li have been projected to level Li+1. For illustration purpose, we overlap the midstructure Mi (skeleton with light colors) with Mi+1 (skeleton with dark colors). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 175 7.9 Handling of GSM bifurcations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 176 7.10 GSMs for rockerarm, fertility and pelvis models. Different GSM sheets are shown in different colors. The sweeping strategy is shown for each model on its boundary. For each GSM, the number of slices and the minimum and maximum are provided below the respective model. . . . . . . . . . . . . . 178 7.11 Comparison of midstructures. Top: GSM; middle: 1D curve skeleton; right: Discrete Scale Axis. Different GSM sheets are shown in different colors. Furthermore, the sweeping strategy is shown on its boundary. The number of slices and the minimum and maximum for the GSM are provided below the respective model. For the curve skeleton, Laplacian constraint scale and positional constraint weight are 2 and 1, respectively. For discrete scale axis, δ = 0.01 and s = 1.1. . . . . . . . . . . . . . . . . . . . . . . . 179 7.12 GSMs computed on a uniform vs. coarse feature-aware triangle mesh. Both use one minimum and one maximum. Curved slices created from the harmonic functions are swept from tail to nose and capture overhang regions consistently. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 180 7.13 Applications of the GSM. (a) Cut through a hexahedral mesh, where the mesh layout was determined by the GSM for the respective model; (b) Deformation based on GSM. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 182 8.1 A nonlinear, hinge-based bending energy [71] is optimized using gradient descent for a genus-0 triangle mesh. Depending on the minimization state, our algorithm (top row) chooses what feature size to minimize. High frequency features are rapidly removed in the initial optimization steps, and lower frequency features are removed in the later optimization steps. By contrast, a nonpreconditioned optimization (shown in the lower row) produces very little shape change as evident in snapshots taken at similar computational times. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 186 8.2 Each of the filters shown on the right can be used to select which eigen-modes of the Laplacian to keep in the gradient during the optimization process. The figures on the left, in which we applied the various filters to the L2 gradient of the Willmore energy, show that the eigenmodes correspond to surface features at different scales. . . . . . . . . . . . . . . . . . . . . 189 xv 8.3 Convergence of multiscale algorithm. Left: One cycle (k = 0, 1, 2, 3, 2, 1), where k is the filter index, to minimize bending energy on the input mesh (top left). The index is increased when the magnitude of changes to the shape falls below a given threshold in order to minimize more global features and then decreased to remove small-scale features introduced by the lower pass filters. Right: Multiscale minimizers based on our imple-mentation of the nonlinear conjugate gradient method and the L-BFGS method exhibit similar convergence rates. While L-BFGS requires fewer iterations, it requires more time per iteration. . . . . . . . . . . . . . . . . . . . . . . 196 8.4 The user-specified tangent constraints (red triangles) are intended to force the input shape (left) to inflate into a hemisphere. A linear method, such as solving the bi-Laplacian in x, y and z, is not able to reconstruct the hemisphere (middle). With our framework, the nonlinear Willmore energy is minimized in a few steps to create the desired shape (right). . . . . . . . . . 198 8.5 Starting from a C2-continuous B-spline surface generated from two input curves (left), we minimize the second-order Willmore energy and the third-order Mehlum-Tarrou energy subject to positional and tangential constraints. A traditional descent-based algorithm using the L2 gradient stalls after a few iterations (middle), whereas the Multiscale Gradient Filtering approach converges to a more optimal solution for both energies (right). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 199 8.6 Starting from the LSCM parameterization, the Multiscale Gradient Fil-tering approach is used to minimize texture stretch. By varying the filter, a minimum is achieved much faster than by just minimizing the energy using the L2 gradient. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 201 8.7 Parameterizations in 2D (top row) and 3D (bottom row). The left pa-rameterization in each row shows regions (circles) where the vector fields of the respective parameterization is not orthogonal. Resulting elements are skewed. The right parameterization of each row shows the improved parameterization. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 203 8.8 We minimize elastic potential energy based on nonlinear strain on an unstructured hexahedral mesh with a trilinear B´ezier basis. Multiscale Gradient Filtering quickly propagates the boundary conditions through the shape, whereas a minimizer based on the L2 gradient stalls. . . . . . . . . 205 8.9 Two animation sequences of an elongated bar. Top row: The elongated bar undergoes deformations based on the infinitesimal strain tensor. The lack of nonlinear terms causes large distortions and ghost forces in the shape. Bottom row: By filtering high frequencies from the forces based on the finite strain tensor, a larger timestep and nonlinear deformation behavior can be achieved. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 207 xvi 9.1 The approximation method proposed in this work constructs more regular curves without introducing additional features. The Frenet frames of the resulting tubes are more well-behaved (top). In comparison, the tubes computed based on a classic interpolation scheme (bottom) makes the tubes look more noisy and irregular. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 211 9.2 Iterative improvement of intput curve (yellow) to the final curve, approx-imating the input points (black). Intermediate curves are drawn in green. The red curve interpolates the input points. . . . . . . . . . . . . . . . . . . . . . . . 212 9.3 Different choices of λ achieve different qualities of approximations. For (a) λ = 0.1 and for (b) λ = 0.8 was used. In both cases = 0.01. . . . . . . . 214 10.1 Our method applied to four representative isosurface visualizations. (a) Vibration modes of a solid structure; (b) Solution to Poisson Equation; (c) Teardrop under nonlinear deformation; (d) Two Isosurfaces of the visible human data set . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 222 10.2 Cubic NURBS curve with nonuniform knot vector and open end conditions.223 10.3 2D analogy: Ray passing through a bivariate NURBS surface with color-coded attribute field A(u) intersecting isocontour at roots of f(t), where the red points refer to entry and exit points with the surface. . . . . . . . . . . 225 10.4 Piecewise trivariate cubic B´ezier patches results in black pixel artifacts, due to degenerate derivative at the B´ezier patch edges. . . . . . . . . . . . . . . . 227 10.5 Ray Frustum/Isosurface Intersection for pixel (s, t) shaded in magenta with adjacent pixels shaded in grey. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 229 10.6 Three ray frustum/isosurface intersection types: 1) Ray frustum and corresponding pixel is fully covered; 2) isosurface silhouette intersects ray frustum with ray intersecting isosurface; 3) Same as 2) but ray does not intersect isosurface. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 231 10.7 Transformation of ray and patch into (α, β, γ)-space. Left: A ray r(t), represented as the intersection of two planes, intersects the isosurface A(u) − ˆa = 0 of Vi(u). Right: Given Vi(u), Ai(u) and the two planes, a new set of coefficients Qk = (αk, βk, γk) are determined to construct Pi(u). The ray intersects the isosurface at uj where |Pi(uj)|∞ < . Pi(u) contains self-intersections and degeneracies depending on the number of intersections. The interior of Pi(u) is illustrated in wireframe. Parts of the (α, β, γ)-space boundary are formed by the interior of the parametric domain. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 235 10.8 OBB hierarchy of patches, referring to a ray/isosurface intersection. With growing subdivision level , the orientation of the OBBs get closer and closer to its parent's orientation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 239 xvii 10.9 Subdivision patches stored in L at subdivision level = 8. In this case, the ray glances the isosurface three times, as shown in Figure 10.7 involving more extensive subdivision and intersection tests. On the left, AABBs were used which result in |L| = 67. On the right, our OBB computation resulting in |L| = 7, significantly reducing subdivision work. . . . . . . . . . . . 240 10.10 S can contain duplicate solutions which can arise due to the scenarios I, II and III. The derivative of the scalar function f(t) is used to filter S to identify unique solutions and solutions representing the same root. . . . . . . 243 10.11 Number of subdivisions per pixel frustum using AABB and OBB for teardrop isosurface from Figure 10.1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 247 10.12 Additional data sets. (a) Unstructured hexahedral mesh (≈ 2.3 million elements) of a segmented torso. Isosurfaces representing voltages of the potential field (using a trilinear basis) are used to specify locations of electrodes to determine efficacy of defibrillation to find a good location to implant a defibrillator into a child. (b) Wake of a rotating canister traveling through a fluid (isosurface of pressure from spectral/hp element CFD simulation data as used in the work [142, 130]). The C(0) nature of the boundaries of the spectral/hp elements can be seen on the isosurface and is not an artifact of our proposed method. . . . . . . . . . . . . . . . . . . . . . 248 xviii LIST OF TABLES 3.1 Test cases for Poisson equation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80 8.1 Performance of the various optimizations of Willmore (W.), Mehlum- Tarrou (M.T.) and Elastic Potential (E.P.) energy, peformend on triangles (T.), B-spline surfaces (B.S.) and B-spline volumes (B.V.). Timings were taken on a quadcore machine. Datasets denoted with (∗) are shown in the supplementary video. In each case, optimization is stopped when shape reaches an acceptable state. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 208 10.1 Average image generation times using OBB and AABB, respectively. The table also shows the timings (in seconds) for each data set when PCA is used instead of our method to compute the OBBs. The degree column presents degrees for the geometry and attribute mapping (tl=trilinear, tc=tricubic, tq=triquintic); μ is the mean; and σ is the standard deviation. The image resolution is 512×512. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 249 ACKNOWLEDGEMENTS First and foremost, I would like to thank my advisor Elaine Cohen. It has been a real honor for me to be her Ph.D. student and academic son. I am very grateful and deeply appreciative of the fact that she was always patient with me and found time to discuss research. She also gave me advice and helped to solve problems related to my dissertation research, concrete feedback on the research papers I was working on, and also the freedom to pick new research problems. These things are not at all natural for a graduate student. Also, I am grateful that she funded me throughout my graduate studies (in part by the grant NSF IIS-1117997). Most importantly, the knowledge of geometry and math that she taught me opened doors to other research fields for me. I would like to thank Mike Kirby for his support and patience to teach me everything I know about the finite element method, which is one of the fundamental building blocks of this work. I thank Richard F. Riesenfeld for his professional advice and wisdom, and for giving me the opportunity to speak in German. I am very glad that I had Charles Hansen in my Ph.D. committee, and I am thankful for his important feedback on several of my research projects and for the opportunity to be a teaching assistant for his computer graphics class. Finally, I am proud to have Tom Lyche on one of my major Ph.D. relevant publications, and for being a member of my Ph.D. committee. Gratitude also goes to Pushkar Joshi and Nathan Carr, who gave me the opportunity to intern at Adobe Systems. The project that we worked on is a part of this dissertation. Many thanks also to the staff of School of Computing, especially to Karen Feinauer and Ann Carlstrom who helped resolve a lot of issues related to my study. I am very grateful that during my life as a graduate student, I found a very good friend in Mathias Schott. I believe that this friendship will last beyond my graduate studies. I am also very grateful that I was surrounded by Suraj Musuvathy and Guoning Chen. It was an honor to collaborate with them. I thank Pooyan Amini for being my gym buddy and introducing me to delicious Persian cuisine, and his brother Peiman Amini for the nice hikes. These friendships made the success of this work possible. Words cannot express my gratitude to my parents, who supported me not only through graduate school, but supported me and made sacrifices for me my whole life. And for Sumathi, who was always there for me even though she was geographically far away. Her love, constant encouragement, and never ending patience with me let me succeed. Thank you Sumi and Mama and Papa, this work is dedicated to you. xxi CHAPTER 1 INTRODUCTION The finite element method (FEM) [86] is a numerical technique to find approximate solutions to partial differential equations (PDE) on an object of interest. It is the de facto standard of many analysis packages which are part of Computer Aided Engineering (CAE) software. Given a Computer Aided Design (CAD) representation of the object of interest, in order to run analysis, the CAD model must be converted using a meshing technique into a CAE representation based, for instance, on hexahedral or tetrahedral elements. This conversion is often time consuming because generally, these CAD models are not watertight, nonmanifold or only represent the boundary of the object of interest and also often require manual work to repair the CAD representation due to these problems. See for instance Butlin et al. [26] for more discussion. A large amount of effort is expended in mesh generation [150] and optimization in also terms of mesh quality [108], a process which frequently takes much longer than the actual analysis. Once analysis has been performed, and the result of the analysis has been evaluated, the original CAD model typically must be modified and the procedure to generate a representation for FEM has to be iterated several times until a satisfactory solution is attained. Proposing to use the same Non-Uniform Rational B-Splines (NURBS) [37] bases that are used to define the parametric geometry as the underlying representation for approximating the solutions to the PDE, Isogeometric Analysis (IA) introduced by Hughes et al. [87] claims to break this time-consuming cycle. In IA which is based on the Galerkin [86] method, simulation parameters, such as forces or material parameters are represented as attributes of the NURBS representation [87]. One of the main advantages of the IA framework is that the user gets feedback directly as attributes of the NURBS 2 model representation, representing geometry exactly, thereby avoiding both the need to generate a finite element mesh and the need to reverse engineer the simulation results from the finite element mesh back into the model. Since IA works directly on the exact CAD representations, geometric refinement is unnecessary and an analyist can focuse only on refinement of the solution. However, as discussed below, one of the main challenges in IA is to create a volumetric shape representation suitable for analysis. By following [87] and [36], the IA pipepline consists of three stages: 1: Modification of the Shape Representation (SR) depending on analysis needs. 2: Analysis on the SR. 3: Investigation of the solution on SR through error analysis and visualization. Go back to stage 1 or 2 as necessary. Research in IA bridges aspects of approximation theory, modeling algorithms, sim-ulation, and visualization. Much research in IA has largely been concerned with issues of the quality of the analysis and simulation results assuming the existence of a shape representation (e.g., trivariate NURBS or T-spline model) that is appropriate for simulation. While Lipton et al. [116] demonstrate that analysis still can be performed successfully on inferior NURBS representations, it will be shown in this work that by carefully creating representations, numerical stability and convergence rates can be significantly improved. There are currently no generally applicable approaches to generating such a representation or visualization methodologies suitable to visualize the solution on the higher-order NURBS representations. One of the main foci of this work lies in proposing methodologies to create CAD representations that are suitable for analysis, so they, in turn, serve as CAE representa-tions as well. We shall call such a representation a CAD/E representation. As discussed below, there are two paths to create a shape representation. In the active modeling path, the CAD/E representation is created ab initio by a modeler using volumetric shape operations. This path is demonstrated, for instance, on a turbine blade in [4]. Already in the design process, tools and methodologies are necessary to produce models which do not require extensive manual repair, since CAD 3 systems in general do not create representations suitable to transform into volumetric representations. An important problem is to assist the user in generating a full volume shape representation and related attribute representation early in the design process that is appropriate for IA, and to maintain it as the design evolves. In the post facto modeling path, the representation is created from a point cloud, CT scan, triangle mesh boundaries, etc., depending on the application. The input as well can consist of multiple boundaries separating materials within the volume. Figure 1.1 illustrates the input for a femur and pelvis data set with inner and outer boundary represented by triangle meshes. Filling up the interior of the model, i.e., model completion, is a hard problem for surfaces, where the input is a set of boundary curves, and is even more difficult for volumes. As a first step, a volumetric parameter-ization is created from the input data. Then, given this parameterization, a NURBS representation is generated from it. Depending on the input, a global parameterization may not be possible, because of higher geometric complexity and topology. Therefore, decomposition strategies are needed to consistenly decompose the input object into a set of regions that are parameterized so that region parameterizations match along common boundaries. In the second step, given a parameterization of the object, data Figure 1.1. Input data. On the left: Triangle mesh boundaries of a femur bone; on the right: Triangle mesh boundaries of a pelvis bone. In both examples, the inner boundary separates cortical from trabecular bone material. 4 fitting techniques are required to create the final shape representation. For that, the parameterization must be sampled so that the resulting shape representions satisfy a specified error tolerance. Creating such a representation for IA is not an automatic approach and is therefore user-assisted. This is important because the user is aware of the parts of the object of interest that are important in simulation and require higher fidelity and so should be free of extraordinary points or degeneracies in the representation compared to other parts of the object. Therefore, various parameters can be specified by the user affecting the simulation result. For instance, given two models representing the model of interest, one may lead to a better rate of solution convergence relative to the other. Since there can be many representations for a single geometry, an appropriate representation should take into account the particular analysis being performed, as a study in [39] has demonstrated. Since in IA, the mesh generation pipeline is eliminated, issues affecting model quality must be considered. That is, model quality is a characterization of those representation properties of the model geometry that impact the model representing the object of interest. However, the longterm goal is to make more of the modeling steps automatic to reduce overall modeling time. Once simulation has been applied to a trivariate representation of a physical domain, in general, mesh extraction techniques or direct visualization techniques are needed to inspect the result of the simulation. This is the third stage of the IA framework. Techniques either directly visualize the solution using, e.g., a ray-based method, or the isosurfaces are extracted and represented as triangle meshes that are rendered. The majority of finite element visualization software, however, displays only linear shape elements. For instance the Marching Cubes technique [119] is based on a uniform grid and only approximates isosurfaces linearly. Although some extensions have generalized it to higher order surfaces, a represention used in IA does not necessarily share the uniform grid structure. A straightforward method is to sample the higher-order geom-etry and solution representation and create a unstructured tetrahedral mesh and use Marching Tetrahedra [33] to extract isosurfaces. However, finding the correct sampling to extract a topological correct isosurface is a complex issue. Furthermore, it violates the IA concept inasmuch as the geometric representation is changed. 5 As mentioned above, other approaches consider a per-ray basis such as [155, 142] and are generally based on classic root finding techniques that use Jacobians or interval arithmetic [134]. NURBS involve an additional nonlinear map where classic methods are difficult to adopt and often fail for complex geometry. Research on direct NURBS visualization has focused on parametric surfaces and trivariate attribute NURBS defined over uniform spatial grids, but there is little work to visualize NURBS isosurfaces for attributes of general NURBS geometry. 1.1 Outline Figure 1.2 gives an overview of the pipeline proposed in this work. Model completion, shape representation and analysis are the fundamental building blocks of this work. The Appendix gives an overview of the publications on which this dissertation is based. 1.1.1 Model Completion As discussed above, there are two long term goals in geometric modeling: The ab initio modeling path where the user creates a CAD/E representation using volumetric shape operations and the post facto modeling path where input boundaries are given and a volumetric representation is created based on these boundaries. Both are formidable problems. This work mostly focuses on the second goal, where in practice a closed triangle mesh in three space is given. This triangle mesh represents the boundary of the physical domain. Commonly, interior triangle meshes are given as well, separating different materials within the physical domain as seen in Figure 1.1. For IA, the goal is to create volume elements in the region enclosed by the multiple boundary triangle meshes. This process is referred to as model completion. The chapters outlined in Figure 1.2, where each one introduces a modeling methodology, are the core of this dissertation. The choice of appropriate modeling methodology strongly depends on the complexity of the input boundaries and the answers to the following questions: What is the genus of the boundaries, what is the size of the features, are the interior boundaries contained within each other, and so forth. 6 Input Shape Representation Model Completion Analysis Optimize Shape Data Fitting Solve Problem on Shape Rep. Direct Isosurface Visualization Modify Choices Output Chapter 8 Chapter 9 Chapter 10 Chapter 3 Generalized Cylinder Chapter 4 Higher Genus Chapter 5 General Object Chapter 6 Generalized Swept Midstructure Chapter 7 Figure 1.2. Overview of the pipeline introduced in this work. In the first step the volume of an input model is completed. Then, this parameterization is optimized and a higher order representation is generated from it. Simulation is applied to the volumetric shapre respresentation which is then visualized and inspected using direct isosurface visualization. 7 In the case where the object of interest exhibits genus-0 and cylindrical-like shape, Chapter 4 proposes a methodolgy to create a single and global volumetric parameteriza-tion of the object that could embody shape overhangs (see femur data set in Figure 1.1). This approach is often desired because decomposing of the object into subvolumes and subsequently recombining them is a complex process. Achieving continuity can be elusive, and decomposing in this manner often introduces extraordinary points. The next level of complexity involves creating a volumetric parameterization for an object with higher genus and structural bifurcations, shapes for which a single parameterization is very difficult, if not impossible, to achieve. For instance, if the object is of genus-1, a single parameterization can be established by sweeping a surface along a closed curve. However, due to the complex shape, the resulting parameterization may contain severe distortion as evidenced in the pelvis data set in Figure 1.1 and detailed in Figure 1.3. Chapter 5 proposes a methodology for modeling more complex objects that can exhibit higher genus and contain multiple bifurcations by consistently decomposing the objects into subvolumes. For the more general case, the input object could be any 2-manifold object, where the methods above are not applicable. In the general case it is very difficult to establish a volumetric parameterization consisting solely of NURBS elements. Therefore, Chapter 6 Figure 1.3. Sweeping can cause significant of parametric distortion, if the object is irregular. 8 introduces a mixed-element methodology that constructs a mixed-element volumetric representation of an arbitrary 2-manifold with high quality semistructured NURBS elements in selected regions and unstructured tetrahedral elements in the other regions. The methodologies discussed in these chapters are each based on a midstructure representation. Since these midstructures are very specific to the respective modeling methods, they are introduced in the respective section with the corresponding modeling method. Chapter 7 introduces a more general framework to compute a midstructure of a given object. We demonstrate that it is suitable for hexahedral meshing and other applications like medial-based shape deformation. 1.1.2 Shape Representation Once the model is completed, a higher order representation such as a trivariate B-spline is fit to an optimized version of the completed model. Chapters 8 and 9 discuss these two pipeline stages. 1.1.3 Analysis Chapter 3 introduces the concept of analysis aware modeling by discussing the classic FEM pipeline and the IA pipeline. It will be demonstrated that different representations result in different convergence behavior. This is central to this work, as the user significanly influences the analysis, i.e., each choice made by the user during the modeling process results in a different analysis result. 1.1.4 Visualization of Analysis Result Finally, once simulation has been applied to the representation, the solution is examined via isosurface visualization. Chapter 10 presents an approach which allows the user to directly inspect a solution without generating an intermediate representation of the volumetric object. The following section introduces core concepts such as parameterization and har-monic functions used throughout the rest of this document. Chapter 2 reviews related work in the fields of parameterization, visualization, mesh generation and mesh quality. 9 1.2 Background This chapter briefly presents core concepts required for this work. 1.2.1 Discrete Harmonic Functions In this work harmonic functions are used to establish a volumetric parameterization over a domain Ω respecting inner material attributes. In general, a harmonic function is a function u ∈ C2(Ω), u : Ω → R, with boundary ∂Ω satisfying Laplace's equation, ∇2u = 0, (1.1) where Ω ∈ Rd and ∇2 = d i ∂2 ∂x2i is the Laplace operator. In our case d = 2, 3, i.e., when d = 2 then ∇2 = ∂2 ∂x2+ ∂2 ∂y2, where ∂Ω is a polyline. When d = 3, then ∇2 = ∂2 ∂x2+ ∂2 ∂y2+ ∂2 ∂z2 where ∂Ω is a triangle mesh. u satisfies the maximum principle, i.e., it does not exhibit any local minima and maxima and is used in a variety of applications. For instance in [91], harmonic functions are used for path-planning of complex robots in complex geometric models. The configuration space represents all degrees of freedom of the robot in the scene, where the location of the robot is fully described by a point. Harmonic functions are computed in this space, where paths along the gradient field of the harmonic function are computed, guiding the robot through the scene without colliding with obstacles. Figure 1.4 shows a cut-through of the harmonic function computed on a tesselated configuration space example. The figure also shows the corresponding paths. The maximum principle makes harmonic functions suitable for tensor-products like parameterization as shown in [50, 196, 126]. In this work harmonic functions are utilized in order to fit a trivariate tensor product B-spline to a tetrahedral mesh generated from a set of triangulated isosurfaces. While in [143] harmonic functions are used for a robust Morse analysis [132] on triangle meshes, here they are used for geometric volume modeling. Chapter 4 uses harmonic vector fields, where the harmonic properties allow a consistent extraction of paths to construct a hexahedral representation (i.e., no self-intersection between adjacent paths) even when the boundary triangle mesh has overhangs (local concavities). 10 (a) (b) robot start positions robot end positions Figure 1.4. In (a), a harmonic function is computed through the volume of a configuration space. In (b), based on the gradient field of the harmonic function in (a), robot paths from start to corresponding end location are computed. Herein, the finite element method (FEM) [86] is used to discretize Equation 1.1. The domain Ω is exactly represented with a tetrahedral mesh (H, T , V, C), where H is the set of tetrahedra. T , the set of faces of the tetrahedra in H. V ⊂ R3, the set of vertices; and C specifies the connectivity of the mesh. V can be decomposed into the sets VB and VI, where VB is the set of boundary vertices that lie on the Dirichlet boundaries that may correspond to the exterior triangle boundary or any interior triangular boundary. VI is the set of vertices for which the solution is sought. The solution, then, is of the form, u(x, y, z) = vk∈VI ukφk(x, y, z) + vk∈VB ukφk(x, y, z), (1.2) where φi(x, y, z) is the linear hat function associated with vertex vi evaluating to 1 at vi and to 0 at the other vertices defining the respective tetrahedron. The weak Galerkin's 11 formulation is used to form the linear system S u = f, where S is the stiffness matrix and f is a right-hand-side function. Being positive definite [86], S lends the linear system to being solved efficiently using the preconditioned conjugate gradient method [12]. The value of u at a point p = (x, y, z) inside a tetrahedron is the linear combination u(p) = 4 j=1 ˆujφj(p), where ˆuj is associated with vertex vj of the respective tetrahedron. Note, that j is a local index. The gradient field ∇u over Ω is piecewise constant, i.e., ∇u defined over a tetrahedron is the linear combination ∇u(p) = 4 i=1 ˆuj∇φj(p), where ∇φj(p) is constant. 1.2.2 Parameterization Given domain Ω ∈ R3, a parameterization is a function F : Θ → Ω, (1.3) where F is a bijection and Θ ∈ R3 is the reference domain. In the discrete case, for instance, two tetrahedral meshes HΘ and H are isomorphic if there is a correspondence between their vertices, edges, triangles and tetrahedra such that corresponding edges join corresponding vertices, corresponding triangles join corresponding vertices and edges and corresponding tetrahedra join corresponding vertices, edges and triangles. As an example, the reference domain Θ of a trivariate B-spline is a single rectangular parallelepiped. On the other side, the domain of a polycube map [193] consists of many rectangular parallelepipeds, which is more similar to its corresponding physical domain Ω. Given such a parameterization, if Ω is to be represented by a higher-order trivariate B-spline, e.g., for purposes of applying physical analysis to Ω, its domain must be decomposed into a set of rectangular parallelopipeds where each one maps to some piece of Ω. Their images intersect only on parallelopiped boundaries. Since the collection of cubes is irregular, there can be extraordinary points on the surface so the T-spline generalization T-NURCCs [178] is more natural in this context, as it allows more complex reference domains. 12 In general, parameterization strategies search for ways to decompose Ω into sub-domains, where each can be parameterized independently. This makes the strategy more flexible in order to represent more complex geometries with reduced parametric distortion. However, establishing matching parameterizations and continuity across the boundary surfaces of the subdomains is often more difficult to achieve. For instance, in the surface case, Ray et al. [165] propose a global parameterization method from input vector fields where the cuts of Ω (i.e., the topology of the base complex) emerge simultaneously from the global numerical optimization process. In this work, functions satisfying Laplace's Equation, i.e., harmonic functions as discussed above are used, to decompose Ω (if necessary) into subvolumes and establish a parameterization of each subvolume so that adjacent subvolumes have a matching parameterization where they connect while respecting interior material attributes in the parameterization. 1.2.3 Orthogonality Two vector fields ∇u and ∇v are orthogonal, if the angle between two corresponding gradient vectors is π/2, i.e. ∇u,∇v = 0 everywhere on Ω. Orthogonality is important in the case of quadrilateral remeshing [196]. Two orthgonal vector fields result in rect-angular quadrilaterals that behave better in numerical simulation than quadrilaterals that are more distored or close to being degenerate. For instance, orthogonal vector fields can be established on a domain with right angle corners such as a cube. A 2D analog is shown in Figure 1.5. Let us call the sides of a discretized rectangle a, b, c and d. Side a and c are opposite and side b and c are opposite. In the following, u = 0 is assigned to a, u = 1 is assigned to c, v = 0 is assigned to b and v = 1 is assigned to d. Laplace's Equation is solved for u and v independently. The resulting vector fields ∇u and ∇v are orthogonal. Figure 1.6 shows orthogonality on more general domains in 2D. The left domain consists of right corners and therefore orthogonal vector fields can be established. Since the right domain in Figure 1.6 consists of a curved boundary, orthogonality cannot be established. 13 d a b c Figure 1.5. Orthogonal vector fields on a square domain. Figure 1.6. Parameterization on domains in 2D. On the left: Since the domain consists of right angle corners and straight lines, two orthogonal vector fields, resulting in a parameterization without any parametric distortion can be established. On the right: No orthogonal vector fields can be established, because the boundary of the domain is curved. CHAPTER 2 RELATED WORK This Section presents work that has been done in the fields of model parameterization (Section 2.1), isosurface visualization (Section 2.3), mesh generation and mesh quality (Section 2.2), relevant for this document. 2.1 Model Completion, Parameterization and Meshing Parameterization is a hard problem for surfaces and even more so for volumes. In addition to use in modeling and remeshing, surface parameterization techniques have a wide variety of applications including texture mapping, detail transfer, fitting and morphing. For a more detailed description, please refer to the surveys [181, 61, 83]. Surface parameterizing techniques such as [117, 70, 196] deal with surface related issues and are not designed to be extended to model volumes. For instance, the authors in [6] motivate anisotropic remeshing and align mesh elements using the principal direction of curvature of the respective triangle mesh. Their approach yields a high quality quadrilateral mesh that has no relationship to the interior. If one were to offset the mesh in the normal direction, one would quickly get self intersections among the elements. Even if this can be avoided, eventually the hexahedral elements have degeneracies and eventually touch each other without proper alignment. Requiring parameterization of the interior makes this problem even more difficult since it is prone to self intersecting offsets and has to deal with skewed and twisted parameterizations. The research on mesh generation is vast, and here, only the most relevant papers will be referenced. For fundamental algorithms and meshing methods, consult the survey [150]. In Chapter 6, a mixed element meshing methodology is proposed that is a combination of structured and unstructured mesh generation, based on a sampled 15 midstructure of the object. The method is related to volumetric parameterization, unstructured meshing, midstructure methods, surface offsetting, and mixed element representations, all of which are reviewed briefly in the following. Decomposing a domain of interest into a set of simpler subpatches involves "patch gluing" where a certain level of smoothness along the patch boundaries is desired. In [117], quadratic B-splines are generalized to fit arbitrary meshes creating triangular and rectangular surface patches. The lack of structure in the irregularities makes it clear that volumetric extensions do not immediately follow. Similarly, manifold splines [70] extend B-splines to surfaces of arbitrary topology, by modeling the domain of the surfaces with a manifold whose topology matches that of the polyhedral mesh. Then it embeds this domain into 2-space using a basis-function/control-point formulation. The domain of this technique is more complicated than the domain of a standard tensor product surface. As in [117], this approach also generates spline patches and glues them together by overlapping them, to get a "match" in the parameterization. The "glue" consists of mathematical operations such as control point constraints. In the case of a volume, patch boundaries are surfaces. Establishing smoothness and continuity between them is very difficult. In [144] the CubeCover algorithm is described which decomposes the volume enclosed by a triangle mesh into hexahedral elements. While the elements are of good quality, the structure depends on a user-defined metamesh. Therefore, in general, the output is an unstructured hexahedral mesh. While being able to reconstruct some types of real world objects, the generalized cylinder(GC) [20] approaches and the approach by [90] are limited, because they require planar cross-sections. Note, representations that rely on planar cross-sections fail for objects with overhangs such as the femur in Figure 1.1. GC is a subclass of our method (Chapter 4) since we are able to model objects with overhangs as the femur in Figure 1.1. Furthermore, interpolation introduces oscillations on the B-spline, which is amplified when dealing with volumes, as in our case. And, for instance, in the method presented in Chapter 4, once a certain stage of the process is reached, the rest proceeds automatically without further user input. Harmonic volumetric mappings between two solid objects with the same topology 16 have been used in a variety of instances. Using related techniques, [202] and [115] are performing a 3D time-variant harmonic deformation from one volume to another volume with the same topology. In a diffeomorphic way, [202] applies this method to brain data. Midstructures such as the medial axis and curve skeletons play an important role in mesh generation. For instance, Sheffer et al. [180] uses the embedded Voronoi graph of the object to decompose the object into simpler sweep-able subvolumes to automatically create a hexahedral mesh. Similarly, Armstrong et al. [10] shows that the medial axis of an object can be used to automatically identify features in mesh generation, dimensional reduction and detail removal. A computed simplified medial axis is then quadrangulated and extruded to construct a hex dominant mesh. Since the medial axis is part of the resulting mesh, it must have consistent topology and connectivity often requiring tedious manual efforts. The approach presented in Chapter 5 is also based on a midstructure. However, the midstructure can be represented by a set of points lying on it, significantly reducing the time to generate it. Furthermore, contrary to the approaches above, the method proposed here maintains the input surface representation in the resulting volumetric representation. In [170] a medial axis-based mesh generator is described. After the construction of a simplified medial axis, quad dominant meshes on the medial axis are generated and extruded to the boundary by advancing front schemes. In this case, advancing front schemes can cause inconsistency cases and for a more complex model, at least in part because it is difficult to control the front, e.g., the front may intersect itself resulting in degenerate elements. Furthermore, in addition to hexahedra, the resulting meshes contain prisms, pyramids and tetrahedrons. In a similar approach Sheffer et al. [180] propose a method for automatic hexahedral meshing based on the embedded Voronoi graph, containing the full symbolic information of the Voronoi diagram and the medial axis of the object. It is used to decompose the object into sweepable subvolumes. Their approach is tested on CAD models with relatively simple medial axes. A more general model with a much more complex medial axis containing many small features (e.g., Figure 1.1) might result in many subvolumes thereby reducing the structured volume behaviour of the resulting mesh. 17 Polycube-Maps were originally proposed by Tarini et al. in [193] to remove seams in texture mapping by making the texture topology of the mesh compatible with the texture domain. Polycubes have been successfully used to construct manifold and polycube splines [73, 200] where the main challenge is to construct the polycube surface mesh for a given object. Wang et al. [201] propose a user-controllable framework where the user directly selects the corner points of the polycubes on the original 3D surface and based on that choice create the polycube maps applying discrete ricci flow. Polycubes have been suggested for volumetric parameterization [79, 69] but have not been presented for a variety of volume models including those that can contain interior material boundaries to which data fitting has to be applied to construct a representation for simulation. Harmonic functions are holomorphic 1-forms as defined in Arbarello et al. [9] and were first introduced by Gu et al. [74] for use in surface parameterization. In this dissertation, harmonic functions are used in combination with appropriate midstruc-tures to attain a parameterization that is consistent with respect to inner material boundaries in the parameterization. Harmonic functions have been shown useful in various applications, especially in the modeling and meshing community. Dong et al. [50] uses them to remesh triangle meshes into quadrilateral meshes with arbitrary topology where the harmonic functions are generated from user specified critical points on the triangle mesh. In [49], critical points are automatically determined by using Laplacian eigenfunctions defined on the surface. Tong et al. [196] automates the surface parameterization problem with discrete differential forms. The aim of a volumetric parameterization method is to create as much interior structure as possible for a given input triangle mesh. The structure allows fitting smooth volumetric patches, e.g., trivariate B-splines to the geometry and simulation can be applied to the resulting volumetric representation. Volumetric parameterization is time consuming, especially because it can require significant manual user input. On the other hand, surface parameterization methods (e.g., see survey by [61] and the methods discussed above) are almost automatic and can produce well-shaped quadrilateral sur-face patches [49] and allow automatic B-spline fitting. A B-spline surface representation 18 created this way can be used as input to the mixed element parameterization method introduced in Chapter 6, creating semistructured trivariate B-spline patches at the boundary of the object and unstructured tetrahedra in the inside. It requires minimal user input and so reduces modeling time significantly. The user can specify the depth of the trivariate boundary patches and has therefore the control of how much of the volume can be filled up with these higher-order elements. It will be demonstrated in the Chapter 6 that simulation scenarios exist where such representations yield stable and high-quality simulation results. Mixed element methods, combining prisms and hexahedra with tetrahedra, have been used in various engineering scenarios such as aerospace applications [99], geo-physics applications, computation fluid dynamics applications [89, 207], just to mention a few. Given a mixed element meshing approach, if the input surface is a quadrilateral mesh, pyramids are used to transition from boundary hexahedra to tetrahedra. If the input is a triangle mesh [96, 203], boundary prisms can be directly connected to interior tetrahedra. Another relevant mixed element approach is the H-morph [151]. Based on an advancing front scheme, an input tetrahedral mesh is iteratively converted into a hexahedral mesh. Tetrahedral elements which cannot be converted into hexahedral elements are connected to its adjacent hexahedral elements using pyramids. In these cases, C(0) continuity between two element types is maintained. The approach presented in Chapter 6 does not enforce such continuity and in that respect is similar to, and motivated by approaches discussed in the introduction of Chapter 6.2.1, where the nodes defining the linear tetrahedra are linked to these higher-order elements. We show that the resulting representation yields convergence under mesh refinement. It will also be demonstrated that our mixed element approach has a similar convergence rate to a representation which uses only tetrahedral elements. However, as will be demonstrated in Chapter 6, the mixed element representation behaves more numerically stable in certain simulation scenarios. In a similar vein to the approach presented in Chapter 6 is that by Wang et al. [96], presenting an approach to generate a mixed element mesh consisting of prismatic and tetrahedral elements. The mixed element representation is constructed from a 19 triangulated surface representing the domain of interest. This triangulated surface is offset along a marching direction. A marching vector of a node is determined through local geometric constraints and a marching step size is chosen to reduce curvature of the previous marching surface. Similarly, [203] uses solutions to the Eikonal equation. Finally, tetrahedra are used to fill up the rest of the object. More recently, [76] first computes a distance field from an input triangle mesh. Based on a user specified thickness parameter, an isosurface of the distance field is extracted and a hexahedral shell mesh based on a polycube domain [193] and harmonic functions is constructed, i.e., the remaining volume in the object is void. Peng et al. [157] proposes a method to add volume textures to an input surface. The method offsets the input surface based on an ODE. The ODE is designed such that the limit offset surface is the medial axis of the input surface. Note, when the input surface contains sharp edges the medial axis touches the input surface. In this case volumetric element tend to have low quality and does not result in a thick shell as desired in our motivation. In the approach presented in Chapter 6, we offset inward a higher-order surface representation of the object, to construct thick higher-order trivariate NURBS elements at the boundary and fill up the interior of the object with linear tetrahedra. We make use of harmonic functions in combination with a midstructure representation to define an offset function for the exterior surface. The use of harmonic functions has several advantages. Firstly, they are flexible as they allow the user to specify any appropriate midstructure, e.g., a point-sampled simplified medial axis or a 1D curve skeleton of the object or an isosurface as used in [76]. Secondly, marching directions and step sizes are implicitly determined and guarantee not to introduce degenerate elements due to the maximum principle of harmonic functions making them easier to control. Thirdly, the use of harmonic functions allow offsetting the input surface deeper into the interior of the object of interest it represents. These advantages make the offsetting approach more robust while also simplifying implementation. Note that offset in this context means to offset based on a harmonic gradient field directions. 20 2.2 Model and Mesh Quality Mesh quality characterizes the geometric shape of the elements in a mesh represent-ing some physical domain. Mesh generation methods (e.g., [186]) producing mesh rep-resentations for traditional finite elements focus on quality measures for finite element meshes. These measures were originally derived from mathematics that characterizes the approximating properties and power of UF to a solution. Those properties depend on C, the geometric mapping F, the particular partial differential equation that must be solved, and the boundary conditions. Those same conditions are embedded in model quality issues of a CAD model suitable for analysis, being the key for analysis aware modeling for IA, however, they are manifested differently. Basic mathematical principles used to define mesh quality for finite element analysis are reviewed in the following. The typical map F for finite element analysis is a collection of mappings from the ideal element (triangle, tetrahedron or square, cube) to each element of the mesh (Figure 2.1), and S is typically the space of linear functions or bi- (tri-) linear functions over the ideal element. If Fi is the mapping from the ideal element to the i−th element, the space U = i UFi forms the approximating space, usually further constrained so that the approximating function space has only C0 elements. Since each Fi is usually affine, U is a collection of local affine linear or bilinear maps. Discretization error and stiffness matrix conditioning are two variables in mesh quality discussions and are therefore important issues for a successful analysis in the finite element method. The discretization error is the difference between the approx-imation computed by the finite element method and the true solution, in this case u(x, y)− u(x, y), measured in an appropriate norm. Stiffness matrix conditioning affects the time and accuracy of analysis and depends on the shape of the elements. A mesh which results in a small condition number better performs using an iterative solver. However, as discussed in [182], this does not imply that the discretization error is small. For instance, for triangles it is known that small angles can cause poor conditioning. However, depending on the PDE, these triangles cause a smaller initial discretization error. In traditional finite elements, it is well acknowledged that the initial mesh upon 21 element aspect ratio maximum angle shear minimum angle reference skew Figure 2.1. Examples for mesh quality metrics in traditional finite elements. one runs one's simulation and refinement algorithms greatly impact the quality of the results. In the best case, the quality of the mesh impacts the constants that exist in the asymptotic error estimates and determine the level of refinement at which the asymp-totic behavior begins. In the worst case, the mesh quality under adaptive refinement is a successively worsening ill-conditioned system which renders the computations useless for engineering practice. Thus, mesh quality is an important issue for finite element methods, when generated meshes will be used in engineering practice for analysis. Because of this reliance of analysis in practice on mesh quality, there is a large body of literature addressing the subject. In particular, the reader is referred to [154]. Given a mesh, the question has to be raised whether its quality is appropriate so that successful analysis can be applied to it. As in the case of B-Splines, h-refinement or knot insertion [37], respectively, do not improve the initial mesh quality. In order to improve the quality of the initial mesh, tremendous effort has been devoted to generate high quality meshes used in finite elements. Generally, schemes are used which relocate the vertex (or node) positions without altering the connectivity of the mesh. Shewchuk in [182] gives bounds on the extreme eigenvalues for Poisson's equation. For instance, λt max, which is the maximum eigenvalue of a local stiffness matrix for 22 triangle t, can be bounded by l2 1 + l2 2 + l2 3 8A ≤ λt max ≤ l2 1 + l2 2 + l2 3 4A , (2.1) where A is the area of t, and li is the length of edge i of t. Note that the maximum eigenvalue of the stiffness matrix can be bounded in terms of the maximum λt max and the maximum number of elements meeting at a single vertex [64]. This allows one to look at the elements t to make a conclusion about the global stiffness matrix conditioning. Therefore, local smoothing techniques also called node-movement strategies are usually used to improve a mesh which relocates each vertex according some objective function to improve the mesh quality in the neighborhood of that vertex. The most commonly used smoothing technique is Laplacian smoothing [59] which moves the current vertex to the geometric center of its incident vertices. However, Laplacian smoothing does guarantee improvement in element quality. This approach is improved by using an objective function which moves the vertex based on a quality metric such as sign of the volume, aspect ratio, angle between adjacent edges (best 90 degrees and worst 180 degrees in case of quadrilateral meshes), stretching orientation, Oddy metric [147] and quality measures for tetrahedral meshes [156]. See Figure 2.1 which shows some of these quality metrics for quadrilaterals. These measures are used to describe the quality of an element. The reader is referred to [107] where these metrics are discussed and applied to different input meshes. Quality metrics like orientation, volume, shape, length ratio and skew are embedded in the element's Jacobian matrix J, which is part of an affine map which maps a point from the reference element to its corresponding physical element (also see Figure 2.1), and can be accessed by applying the QR factorization of J. In the case of triangles, the reference element is the triangle defined by the vertices (0, 0), (1, 0) and (0, 1)). Then, linear maps like the Frobenius norm are used to convert J into a scalar which is then used to design a objective function optimizing certain quality metrics. The reader is referred to the work of Knupp [108] who develops a mathematical theory of geometric quality metrics applied to unstructured meshes. If a choice is made on the objective function, an optimization process is applied to minimize the objective function and 23 therefore mesh is improved. Newton search [107] for instance is an efficient way to minimize the objective function. Note that these method are proposed for a trilinear simulation and geometry basis. In the case of B-splines, it will be demonstrated that these quality measure are not sufficient and new and more appropriate metrics should be developed to describe the quality of the mapping F. Also, due to the tensor-product nature, modeling with NURBS has limitations. Sweeping [37] however, is a powerful technique to model a wide class of objects. To generate a swept volume, a NURBS surface is swept along a NURBS curve. Often, given a "sweepable" physical domain Ω, there are usually more ways for the sweep. Different choices result in degeneracies or more skewed elements at regions where this might be not desirable. As we will see on an 2D example in the results section, different choices affect the discretization error. They also affect the stability of linear systems and therefore eigenvalue problems. 2.3 Isosurface Visualization Visualization techniques are used in numerous engineering fields-including medical imaging, geosciences, and mechanical engineering-to generate a two-dimensional view of a three-dimensional scalar or vector data set. Additionally, they can visualize simulation results (e.g., generated with the finite element method). Consequently, the development of such visualization algorithms has received much attention in the research community. Techniques usually fall into three groups: (1) direct volume rendering, (2) isosurface mesh extraction followed by isosurface mesh rendering, and (3) direct rendering of isosurfaces. Techniques in category (1) typically involve significant computation, especially when dealing with arbitrary geometric topologies represented by high-order basis functions such as NURBS. In ray-based direct volume rendering methods (see [111, 127]), it is necessary to integrate each ray through the volume using sufficiently many integration steps. Each integration step requires an expensive root-solving due to the nonlinear mapping. Hua et al. [84] presented an algorithm to directly render attribute fields of tetrahedral-based trivariate simplex splines by integrating densities along the path of 24 each ray corresponding to a pixel. In the case of uniform grid data sets, accumulating slices aligned along the viewing direction (see [205]) is efficient and commonly used in practice, even though ray-based techniques offer a range of optimizations (e.g., empty space skipping). Methods in category (2) assume a regular grid of data and extract isosurfaces using Marching Cubes (MC) [119], resulting in a piecewise planar approximation of the isosurface. After isosurface mesh extraction, the faces of the isosurface mesh are rendered. Marching Tetrahedra (MT) [33] is applied to both structured and unstruc-tured tetrahedra-based grids. In both MC and MT, the corners of a hexahedral or tetrahedral element, respectively, are used to determine if the isosurface passes through the respective element. Then, the intersections between the element's edges and the isosurface are determined to create piecewise linear facets approximating the isosurface. Although these approaches are efficient and therefore widely used in practice, they approximate the isosurface by piecewise linear facets within an element with some ambiguity, and therefore do not guarantee topological correctness. As an example, Figure 2.2 shows a discretized domain with 300 000 linear tetrahedra. As will be shown in Chapter 10, this domain can be represented with a single triquintic NURBS element. As seen in Figure 2.2a, the respective isosurface extracted with MT has ambiguities in the topology, resulting from data that are known only at the corners of the elements and hence can miss isosurface features. Furthermore, the time to construct the respective (a) (b) (c) Figure 2.2. Discretization of domain from Figure 10.1c with 300k tetrahedra and application of marching tetrahedra (using ParaView). (a) isosurface; (b) scalar field on tetrahedra; (c) our approach on single triquintic NURBS patch 25 mesh representation can be computationally laborious. Schreiner et al. [173] propose an advancing-front method for constructing manifold isosurfaces with well-shaped triangles (Figure 2.3), although it has some difficulties when the front meets itself (the stitching problem). Meyer et al. [130] propose a particle system on high-order finite element mesh (arbitrary geometric topology), which applies surface reconstruction on the particles to construct the isosurface mesh; however, the visualization produced is not a water-tight surface. When the data are known only at the corners of a hexahedral mesh, our method constructs an approximation by filtering the data with a high-order approximating or interpolating trivariate B-spline filter (see [121]). The filter can be trilinear (only C(0)), tricubic (C(2)), or higher degree, as required by the user. Then, an isosurface of the high-order approximation is directly rendered with pixel accuracy. In category (3), the isosurface is rendered directly, i.e., for every pixel on the image plane, its corresponding point on the isosurface is determined (Figure 2.3, left). Once the point on the isosurface for a given pixel is known, the pixel can be shaded using the gradient as the normal for the given point. Another motivation to visualize specific isosurfaces is to color-code information, such as material density, to get a better Figure 2.3. Isosurface from silicium data set (volvis.org), isovalue of 130 using Marching Cubes (using ParaView), Afront (ρ = 0.3) and Direct visualization with our proposed method. 26 understanding through which materials the isosurface passes. Knoll et al. [106] use a trilinear reconstruction filter on a structured grid and a ray-based octree approach to render isosurfaces and achieve interactive frame rates. Nelson et al. [142] propose a ray-based isosurface-rendering algorithm for high-order finite elements using classic root-finding methods, but do not consider element curvature (i.e., the multiple entry and exit problem). Kloetzli et al. [104] construct a set of structured B´ezier tetrahedra from a uniform grid to approximate any reconstruction filter with arbitrary footprint. Given this reconstruction, generated from gridded input data (e.g., medical or simulation data), they directly visualize isosurfaces using the ray/isosurface intersection method presented by Loop [118]. The method proposed in Chapter 10 is most closely related to class (3) approaches, i.e., our proposed method directly visualizes an isosurface from a trivariate NURBS of arbitrary geometric complexity. However, instead of following only a ray-based scheme, our approach computes the intersection between a ray frustum and the isosurface. Fur-thermore, it is often desired to visualize the geometry represented by the NURBS. While approaches similar to the work in [2] can be used to render the object-surface geometry, our approach can be used to simultaneously visualize both the geometry represented by the NURBS and the visualization of isosurfaces of the attribute representation in a robust way. Intersecting a ray frustum with an object in the scene is related to the approaches that propose cone-tracing given in the work [7] and beam-tracing (see [80]) for more efficient anti-aliasing, soft shadows, and reflections. However, both of those techniques deal only with polygonal objects. For isosurfaces of algebraic functions, the dissertation [47] presents interval approaches to create intersection tests in the ray-tracing of implicit surfaces. In particular, it shows a ray sampling-based method to exploit the coherence of rays to accelerate the process of ray-tracing implicit surfaces, which can also be used for anti-aliasing isosurface silhouettes. CHAPTER 3 ANALYSIS-AWARE MODELING IA [87] has been proposed as a methodology for bridging the gap between Computer Aided Design (CAD) and Finite Element Analysis (FEA). Although both the traditional and isogeometric pipelines rely upon the same conceptualization to solid model steps, they drastically differ in how they bring the solid model both to and through the analysis process. The IA process circumvents many of the meshing pitfalls experienced by the traditional pipeline by working directly within the approximation spaces used by the model representation. This chapter demonstrates that in a similar way to how mesh quality is used in traditional FEA to help characterize the impact of the mesh on analysis, an analogous concept of model quality exists within IA. The consequence of these observations is the need for a new area within modeling - analysis-aware modeling - in which model properties and parameters are selected to facilitate IA. The concept of IA was first introduced by Hughes et al. in [87] as a means of bridging the gap between Computer Aided Engineering (CAE), including Finite Ele-ment Analysis (FEA), and Computer Aided Design (CAD). Those familiar with the application of FEA to CAD-based models are well-aware of the complications and frustrations which arise when one attempts to take a "solid model" (a term which we will define in the context of the modeling community in Section 3.3) as produced by a typical commercially available CAD system, generate a surface tessellation and corresponding volumetric representation in terms of meshing elements (e.g., triangles and quadrilaterals on the surfaces and tetahedra and hexahedra in the volume), and run an analysis. On this "preprocessing" side prior to the actually analysis step, a large amount of effort is expended in mesh generation and optimization (in terms of mesh quality), sometimes to the point of consuming more time than what is taken by the 28 actual analysis step. Once an analysis is run, solution refinement often requires mesh adaptation or in worst case regeneration, both of which require consultation with the original CAD-model. IA claims to break this common but insidious cycle by choosing an alternative route from the solid geometric model to analysis. In IA, one works with the functions used to generate the model directly by using the function space used for model generation as the approximating space in which field solutions are built (hence the name iso-geometric). By circumventing many of the pitfalls that one encounters during the mesh gener-ation process by working directly with the solid model, IA effectively eliminates the geometric error component of the analysis pipeline. Geometric refinement is no longer necessary; the analyst can focus attention solely on solution refinement. It is our thesis that although circumventing the mesh generation pipeline implies that one no longer needs to consider mesh quality, there are still issues of model quality that must be considered. In a similar way to how mesh quality is a geometric means of assessing the impact of a mesh on the function space which it induces in the classic finite element process, model quality is a characterization of those properties of the representation of the model geometry that impact the representation space (or trial space) used to approximate the fields of interest. The consequence of these observations is the need for a new area within modeling - analysis-aware modeling - in which model properties and parameters facilitate IA. 3.1 Nomenclature In this section we set up the environment for considering IA in the context of linear second-order partial differential equations (PDEs) with zero Dirichlet boundary conditions. We note that there exists a straightforward extension to nonzero boundary conditions and also to Neumann boundary conditions. The issues we raise are quite general and will arise in using the isogeometric concept in solving many types of partial differential equations: Some examples we will consider and upon which we will comment in the examples section will have nonzero boundary conditions and/or more complex partial differential operators (such as those found in the modeling of linear elasticity). 29 We do, however, set up here the nomenclature to illustrate these specific problems in an arbitrary number of space dimensions. Although some of these terms may appear obvious to either those familiar with ge-ometric modeling or those familiar with engineering analysis, we believe it is important for both the geometric modeling and finite element analysis communities to be overtly explicit during this time of confluence of ideas. Let Ω ⊂ Rs with s ∈ N be a bounded domain with boundary ∂Ω. Ω is the physical domain, often called the world space or physical space. Using the notation Dj = D1 j to denote the partial derivative with respect to the j-th variable, define Dα := Dα1 1 · · ·Dαs s , (3.1) a mixed partial derivative of total order |α| = α1 + · · · + αs. Then the column vector ∇f = [D1f, . . .,Dsf]T denotes the gradient of f. Further, let H1 = H1(Ω) := {f : Ω → R : Dαf ∈ L2(Ω), |α| ≤ 1} (3.2) V = H1 0 = H1 0(Ω) := {f ∈ H1(Ω) : f = 0 on ∂Ω}, (3.3) denote the Sobolev spaces of functions with values and first order partial derivatives in L2 = L2(Ω). Let {φi}ni =1 ⊂ H1(Ω) be linearly independent functions in H1(Ω). Moreover we define J := {j ∈ {1, . . . , n} : φj ∈ V } and |J| = #{j : j ∈ J}, 30 that is, the set of indices of those φj that vanish on ∂Ω. If s = 1 then typically J = {2, . . . , n − 1}. We define the space V q h := { j∈J cjφj : cj ∈ Rq, j ∈ J, q ∈ N}, (3.4) and note that Vh = V 1 h is a subspace of V . The index h is a flag indicating finite dimensionality and is often a measure of element diameter. The space V q h forms the space in which the approximation to the solution of the differential equation is made. Suppose {ψj}nj =1 is a set of real-valued linearly independent functions on a partition of the unit cube Θ = [0, 1]s in Rs, and the functions φj are given as φj(x) = ψj ◦ F−1(x), j = 1, . . . , n, (3.5) where F = (F1, . . . , Fs) : Θ → Ω is a bijection. Figure 3.1 illustrates a modification between the shape of a ψ and its corresponding φ induced by F−1. Moreover, we assume that F(Θo) ⊂ Ωo σ := F|∂Θ : ∂Θ → ∂Ω, (3.6) i.e., F maps interior to interior and boundary to boundary. If we use the same functions ψj to define both the φj's and F F = n j=1 γjψj , for some γj ∈ Rs, j = 1, . . . , n, (3.7) then the approach to solving the partial differential equation is called isogeometric. Typically analysis will be carried out on models whose definition is piecewise over multiple hypercubes Θ = {Θi = [0, 1]s : i = 1, . . .,K} for some finite K ∈ N, and F is defined piecewise in terms of the mappings from each Θi such that, for all i, j, i = j, 31 y x t s Figure 3.1. F maps a point from the reference domain Θ to the physical domain Ω. Correspondingly, F−1 maps a point from Ω back to Θ. With F−1, it is possible to define basis functions on Ω as compositions of F−1 with basis function defined on Θ. F(∂Θi) F(∂Θj) ⊂ ∂F(Θi) ∂F(Θj). Further F must be continuous on every nonempty intersection. Θ is called the para-metric domain or the reference space. CAD systems typically create representations of the model that define mappings σ from ∂Θ to ∂Ω, and can be written as a collection of mappings σi : [0, 1]s−1 → Rs, i = 1, 2 . . . , 2s that agree on their shared boundaries. When s = 3, this representation is sometimes suitable for performing isogeometric shell analysis, but in order to perform a full volumetric analysis, the model must be completed. That is, the representation must be extended to completely define the interior so that V q h is defined. We will discuss this further in what follows as this is often a nontrivial process. Typically, ψ's are tensor product B-splines, Non-Uniform Rational B-Splines, rect-angular subdivision surfaces or T-splines. 32 3.2 Outline This chapter is structured as follows. In Section 3.3 we present the modeling to analysis pipeline. We briefly describe the process of going from model conceptualization to the solid model, and then distinguish between the classic route of surface and volumetric mesh generation as done for classic finite element analysis from the boundary representation and volumetric model generation as done for IA. In Section 3.4 we provide the mathematical and algorithmic descriptions of the isogeometric methodology em-ployed in this work. In Section 3.5 we take a step back so as to appreciate the parametric modeling of geometry from the perspective of, and in the language of, isogeometric analysis. By doing so, we hope to demonstrate that in general there is not "a model" (a single ideal model) on which one does IA, but rather that designers are presented with a collection of modeling choices - some of which may inadvertently impact analysis. In Section 3.6 we present one-, two- and three-dimensional examples comparing different model completions and demonstrating the impact of model completion on quality of the solutions one obtains from an IA. We conclude in Section 3.7 with a summary of this work, some conclusions that we can draw and proposals of future work. 3.3 The Modeling to Analysis Pipeline In this section, we review the classic model to mesh to analysis pipeline as appre-ciated by most FEA researchers, and then provide the corresponding modification to the pipeline as introduced by IA. Note that we pay particular attention to the use of nomenclature in this section, as the confluence of concepts from two fields (modeling and analysis) has led to misconceptions in both fields as to what is being discussed. We will use Figure 3.2 as our visual guide through this process. 3.3.1 Conceptualization to Solid Model The stages of the pipeline from conceptualization to solid model are denoted by the left half of Figure 3.2. The designer has in mind a concept or ideal of what is to be designed, and uses a CAD modeling system to construct a collection of surfaces that are meant to represent the outer boundary of the object of interest. Note that the modeler is not working with three-dimensional manifold representations, but rather is working 33 Surface Representation "Solid" Model Volumetric Model Trim Surface Tesselation Volumetric Tetrahedralization Mesh Quality Model Quality Boundary Model Repair Repair Repair Figure 3.2. Diagram presenting the classic concept-to-mesh pipeline (top branch) and the concept-to-model pipeline (bottom branch). A detailed discussion of the diagram is presented in the text. 34 with surface subregions that are intended to bound the object of interest. These surfaces are often constructed one-by-one without regard for how they will connect, intersect or overlap with other pieces. The modeler then uses the CAD system to accomplish what is referred to as trimming, an attempt to connect (or stitch) the surfaces together to form a watertight model, that is, one that clearly delineates R into three regions, inside the model, outside the model, and on the boundary of the model. Within the shape modeling community, the term solid model is used to characterize such a representation. When the solid model is represented using pieces of bounding 2-manifolds, the represen-tation is called a boundary representation or b-rep. Whereas the surface representation (pretrimming) does not necessarily faithfully represent the geometric and topological properties (such as being a water-tight surface) on the conceptual object, the newly formed solid model should. At the conclusion of this process, a solid model is output from the CAD system. Although called a solid, it is a collection of pieces of surfaces and connectivity information that are intended to bound a three-dimensional object and that is intended to be watertight. Unfortunately the intersections between sculptured surface pieces that define the curves along which the pieces should be trimmed and stitched together cannot be exactly represented in the parameter spaces of the defining surfaces, but rather are defined implicitly. Hence, while CAD systems have different approaches to explicitly representing these curves, the trimmed surfaces and resulting b-rep models are all approximated along the trimming edges. Frequently it may be necessary to repair the model to make it suitable for later processes such as analysis or fabrication. 3.3.2 Traditional Meshing Pipeline Leading to Analysis The traditional meshing pipeline leading to analysis is denoted by the upper branch of the right half of Figure 3.2. In the figure, we have purposefully placed quotations around solid to draw the reader's attention to two things. First, as previously men-tioned, the solid model is not solid (in the sense of the term as used by analysts), but rather denotes the boundary of the object. Secondly, the solid model as produced by CAD systems is not always truly water-tight, but possibly only visually water-tight. 35 This issue has been the bane of many surface tessellation efforts which have devised schemes under the assumption that the solid model formed a mathematically water-tight (i.e., closed in the topological sense) representation. In going from the solid model to a surface tessellation, it is often necessary to invoke a repairing procedure. We mark the repair process as being the point of deviation between the traditional pipeline and IA as many of the repair procedures used in traditional mesh generation assume that the target representation is a piecewise linear tessellation. The result of the repair and surface generation process is a tessellation of boundary of the object of interest. This tessellation is an approximation to the true geometry (in this case, the CAD-model), where approximation decisions have been made both in terms of how repairs are done and in how finely the tessellation captures the features of the original model. If a three-dimensional analysis is desired, the next step is to generate a volumetric tessellation, normally by filling in the volume with elements of the appropriate type (hexahedra or tetrahedra) for the analysis of interest. In the classic finite element procedure, one generates a tessellation which approxi-mates the true geometry, and then uses this tessellation to induce a function space in which approximations will be made. For instance, in classic linear finite elements over triangles, the triangular tessellation induces a piecewise linear (in total degree) space which is C0 continuous along the edges of the triangles. As is well known by finite element practitioners, given two tessellations, both of which faithfully represent the geometry, one can get drastically different solutions due to the properties (or richness) of the approximating space that is induced. A natural feedback loop developed between analysts and mesh generation experts concerning the impact of meshes on solution quality. These metrics have commonly become known as mesh quality metrics. That is, they are geometric considerations (normally involving things like ratios of angles of elements, aspect ratios of edges, etc.) which help guide the development of meshes appropriate for analysis. Although these metrics have deep foundations within approximation theory, they are often abstracted away so that only geometric qualities of the mesh are discussed. We remind the reader, however, that maximizing mesh quality is, in its essence, an attempt to positively shape 36 the approximating function space induced by the mesh. 3.3.3 Isogeometric Pipeline Leading to Analysis While IA is still a young field, the authors hypothesize that isogeometric pipeline leading to analysis is denoted by the lower branch of the right half of Figure 3.2. As in the case of the traditional pipeline, repair is needed to ensure that the model being used for analysis meets the required topological constraints (such as closure) of the problem. In the case of IA, however, this repair process must be done keeping in mind the original and target representations. A starting point for isogeometric discussions in line with the finite element approaches is the boundary model, which should be a geometrically and topologically correct model of the bounding surfaces of the object. If a three-dimensional analysis is desired, volumetric representations must be generated prior to the analysis. The approximating space generated during an IA is dependent upon the boundary model (in 2D) or volumetric model (in 3D) that is used. Just as in the case of classic mesh generation, two different volumetric models generated from the same boundary model will create two different approximating spaces. Analogous to mesh quality impacting analysis, model quality impacts IA. A different starting point for IA is that consideration during the shape (usually boundary) modeling process should be given to create a representation that lends itself to IA. There are frequently many modeling operations that lead to different representations for either the exactly same or closely related boundary shapes. Some of those representations are better suited to analysis than others, and within those groups, some are better suited to certain types of analysis than others. Similar issues have been recognized in created representations for models suitable for the multitude of computer-aided manufacturing processes and techniques. However, progress has been made in developing CAD systems that develop representations that, while suitable for design and display, are fabrication-aware, thus enabling a smoother, faster transition between design and fabrication. Analysis-aware modeling in the context of IA may prove to be a key step towards that progression for design, engineering analysis, and simulation. Towards that end, this chapter raises several important issues through a combination of analysis and 37 demonstration in which the interaction between representation and analysis can either enhance or make the product evolution process difficult. Until such time as these issues have been quantified and embedded in analysis-aware modeling systems, the human modeler must be mindful of them. 3.4 Mathematical Formulation In this section, we first review the basic mathematical representational building blocks on which IA as well as many CAD and geometric modeling systems represent geometry. An overview of NURBS (Non-Uniform Rational B-Splines) can be found in [37]. All computational algorithms are presented there, so in the following section we discuss definitions of B-spline and NURBS functions and their combinations to define parametric mappings of global geometry. Note that this discussion provides the mathematical building blocks of modeling, but does not address how these building blocks are assembled as part of the modeling process. We will delve into the mind of the modeler in a subsequent section (Section 3.5). 3.4.1 The Framework Let Ω ⊂ Rs with s ∈ N be a bounded domain with boundary ∂Ω. The symbols Dα, H1, H1 0, and ∇ are defined as in Equations (3.1), (3.2), and (3.3), respectively. These Sobolev spaces have a norm given by u 21 = Ω u(x)2 + ∇u(x)T∇u(x) dx. Let a : H1 0 × H1 0 → R denote the bilinear form a(u, v) := Ω ∇u(x)T∇v(x)dx, (3.8) This bilinear form is positive definite on H1 0, and H1 0 is a Hilbert space with inner product a(u, v) and associated norm u = a(u, u). We let 38 (u, v) := Ω u(x)v(x)dx, u, v ∈ L2(Ω), be the usual L2 inner product. For a general review of Sobolev spaces we refer the reader to [175]. Discussions in this chapter will mostly be based on problems that arise in the relationship between geometry and analysis models. Most studies are focused on two prototypical mathematical model problems that arise in analysis, the Poisson problem and a corresponding eigenvalue problem, respectively. −∇2u = f on Ω, u= 0 on ∂Ω, (3.9) −∇2u = λu on Ω, u= 0 on ∂Ω. (3.10) Given f ∈ L2(Ω) the weak form of (3.9) is to find u ∈ V := H1 0(Ω) such that a(u, v) = (f, v), v∈V . (3.11) It is well known that (3.11) has a unique solution u, see [32]. The weak form of (3.10) is to find λ ∈ R and a nonzero u ∈ V such that a(u, v) = λ(u, v), v∈V . (3.12) Since this involves finding the eigenvalues and eigenfunctions of a compact, symmetric operator in the Hilbert space (V, a(·, ·)) there exists an increasing sequence of strictly positive eigenvalues 0 < λ1 ≤ λ2 ≤ · · · ≤ λk ≤ · · · 39 with lim λk = ∞ and associated eigenfunctions uk, which can be orthogonalized so that a(uj, uk) = λkδj,k, j,k≥ 1. Moreover, the eigenfunctions form a complete system, i.e.,the set of all linear combina-tions is dense both in V and L2(Ω), again see [32]. For V q h defined as in Equation 3.4, the Galerkin method for (3.11) consists in finding uh ∈ V q h such that a(uh, vh) = (f, vh), vh ∈ Vh. Writing uh = j∈J cjφj , we obtain a linear system for the unknown coefficients c. Sc = f, S = a(φi, φj) i,j∈J , f = (f, φi) i∈J . (3.13) Since the φ's are linearly independent and vanish on ∂Ω the stiffness matrix S is a symmetric positive definite matrix and (3.13) has a unique solution that is amendable to iterative methods like conjugate gradient [12]. The Rayleigh-Ritz method for (3.12) consists of finding λ and a nonzero uh ∈ V q h such that a(uh, vh) = λ(uh, vh), vh ∈ Vh. Writing uh = j∈J cjφj as before we obtain a generalized eigenvalue problem Sc = λMc, S = a(φi, φj) i,j∈J, M = (φi, φj) i,j∈J, (3.14) where both S and the mass matrix M are positive definite. 40 Thus the eigenvalues λkh of (3.14) that approximate the exact eigenvalues of (3.12) are positive 0 < λ1h ≤ · · · ≤ λ|J|h and the eigenfunctions uh can be chosen to be orthogonalized so that a(ujh, ukh) = λkhδj,k. Moreover limh→0 λkh = λk for 1 ≤ k ≤ |J|, provided limh→0 infvh∈Vh uk − vh = 0 for 1 ≤ k ≤ |J|. 3.4.2 Definition of Isogeometric Finite Element Analysis Suppose the basis functions φj(x) are given as in Equation (3.5). Moreover, we assume that (3.6) holds, i.e., F maps interior to interior and boundary to boundary, and that F is defined as in Equation (3.7). Then the Galerkin and Rayleigh-Ritz methods for (3.11) or (3.12) are called isogeometric. The elements sij of the stiffness matrix S can be expressed in terms of the gradients of the ψj basis functions. Let J = JF := ⎡ ⎢⎣ D1F1 · · · DsF1 ... ... D1Fs · · · DsFs ⎤ ⎥⎦ = ⎡ ⎢⎣ ∇FT 1 ... ∇FT s ⎤ ⎥⎦ (3.15) be the Jacobian of F. Note that the elements of J are functions defined on Θ. We assume that J(t) is nonsingular for all t ∈ Θ. Then sij = Ω ∇φi(x)T∇φj(x)dx = Θ ∇ψi(t)TN(t)∇ψj(t)dt, i, j = 1, . . . , n, (3.16) where 41 N = |det(J)|J−TJ−1. (3.17) Note that N(t) is positive definite for all t ∈ Θ. Explicitly, for s = 1, N(t) = 1 |F (t)| , and for s = 2 N = 1 |det(J)| ∇F2 22 −∇FT 1 F2 −∇FT 1 F2 ∇F1 22 . If K := |det(J)|−1/2J = UΣV T , Σ = diag(σ1, . . . , σs), UTU = V TV = I, with σ1 ≥ σ2 ≥ · · · ≥ σs > 0 is the singular value decomposition of K then N = UΣ−2UT is the spectral decomposition of N. Since N is positive definite the eigenvalues of N are the inverse square of the singular values of K and the orthonormal eigenvectors of N are the right singular vectors of K. 3.4.3 B-splines For integers n ≥ 1 and d ≥ 0 let τ = {τi} be a nondecreasing finite sequence of real numbers. We refer to τ as a knot vector and its components as knots. On τ we can recursively define degree d B-splines Bj,d = Bj,d,τ : R → R by Bj,d(t) = t − τj τj+d − τj Bj,d−1(t) + τj+d+1 − t τj+d+1 − τj+1 Bj+1,d−1(t), t∈ R, (3.18) 42 starting with Bj,0(t) = 1, if τj ≤ t < τj+1, 0, otherwise. Here we use the convention that terms with zero denominator are defined to be zero. We let Bd,τ = {Bj,d,τ}j . A B-spline Bj,d of degree d has the following properties: 1. It depends only on knots τj, . . . , τj+d+1 and is identically zero if τj+d+1 = τj . 2. For t ∈ (τj, τj+d+1), 0 < Bj,d(t) ≤ 1 and Bj,d(t) = 0 otherwise. The interval [τj, τj+d+1] is called the support of Bj,d. 3. Its derivative is DBj,d(t) = d Bj,d−1(t) τj+d − τj − Bj+1,d−1(t) τj+d+1 − τj+1 , again with the convention that terms with 0 denominator are set to 0. 4. If m of the τj, . . . , τj+d1 are equal to one value z, then DrBj,d is continuous at z for r = 0, . . . , d − m and Dd−m+1Bj,d is discontinuous at z. 5. Its integral is τj+d+1 τj Bj,d(t)dt = τj+d+1 − τj d + 1 . 6. It is affine invariant, i.e., for u, v, t ∈ R Bj,d,uτ+v(ut+v) = Bj,d,τ (t), where uτ+v := (uτj + v)j . Now suppose n, d are integers with 0 < d < n. We say that τ = {τi}n+d+1 i=1 is a (d + 1) extended knot vector on an interval [a, b] if 43 a = τd+1 < τd+2, τn < τn+1 = b, τi+d+1 > τi, i= 1, . . . , n. It is (d+1)-regular or (d+1)-open if in addition τ1 = a and τn+d+1 = b; it is (d+1)-regular uniform or (d + 1)-open uniform if τi+1 − τi = h for i = d + 1, . . . , n and h > 0. The knot vector is uniform if τi+1 − τi = h > 0 for i = 1, . . . , n + d. On the knot vector τ = {τi}n+d+1 i=1 we can define n B-splines of degree d. The linear space of all linear combinations of B-splines is the spline space defined by Sq d,τ := n j=1 cjBj,d | cj ∈ Rq for 1 ≤ j ≤ n , d≥ 0, q≥ 1. An element f = n j=1 cjBj,d of Sd,τ = S1 d,τ is called a spline function if q = 1 or just a spline of degree d with knots τ, and (cj)nj =1 are called the B-spline coefficients of f. For q > 1 the combination f = n j=1 cjBj,d is a spline curve. Suppose τ = {τi}n+d+1 i=1 is a (d+1)-open knot vector on [a, b]. A spline f : [a, b] → R is by definition continuous from the right. We define f(b) by taking limits from the left. Let f = n j=1 cjBj,d. Then the following properties hold: 1. B-splines (Bj,d)nj =1 are linearly independent on [a, b] and therefore a basis for Sq d,τ . 2. Partition of unity: n j=1 Bj,d ≡ 1, t∈ [a, b]. 3. Convex hull property: f(t) lies in the convex hull of [c1, . . . , cn]. 4. Smoothness: If z occurs m times in τ then f has continuous derivatives of order 0, . . . , d − m at z. 5. Locality: If t is in the interval [τk, τk+1) for some k in the range d + 1 ≤ k ≤ n then f(t) = k j=k−d cjBj,d(t), (3.19) 44 6. Affine invariance: If uτ + v := (uτj + v)n+d+1 j=1 then n j=1 cjBj,d,uτ+v(ut + v) = n j=1 cjBj,d,τ (t), t∈ [a, b], u,v ∈ R. (3.20) 7. Nodal representation: t = n j=1 τ ∗ j,dBj,d(t), τ∗ j,d = τj+1 + · · · + τj+d d , t∈ [a, b]. 8. Derivative of a spline: Df(t) = d n j=2 cj − cj−1 τj+d − τj Bj,d−1(t), t∈ [a, b], where terms with 0 valued denominator are set to 0. 9. Integral of a spline: τn+d+1 τ1 f(t)dt = n j=1 τj+d+1 − τj d + 1 cj . (3.21) 10. If z = τj+1 = · · · = τj+d < τj+d+1 for 1 ≤ j ≤ n then f(z) = cj . 11. Marsden's identity: (s − t)d = n j=1 j+d i=j+1 (s − τi)Bj,d(t), t∈ [a, b], s∈ R. (3.22) 45 3.4.4 Knot Insertion and Degree Raising (h- and p-refinement) Suppose τ is a knot vector. The distinct elements in τ are called break points. We define the multiplicity of z in τ as μτ (z) = #{τj ∈ τ : τj = z}, z∈ R. Notice that μτ (z) = 0 if z is not equal to one of the knots in τ. For k ≥ 0 we define the knot vector τ (k) to have the same break points as τ, and μτ(k)(ξ) = μτ (ξ) + k for all ξ ∈ τ . Thus we increase the multiplicity of each break point in τ by k. If t is another knot vector then we say that τ ⊂ t if each break point ξ in τ is also a break point in t and μτ (ξ) ≤ μt(ξ). Let d, e be integers, 0 ≤ d ≤ e, let τ = (τj)n+d+1 j=1 be (d+1) extended on [a, b] and let t = (ti)m+e+1 i=1 be an (e+1) extended knot vector on the same interval [a, b]. If τ (e−d) ⊂ t then Sd,τ ⊂ Se,t, and there is a matrix A ∈ Rm,n transforming the B-splines in Sd,τ into the B-splines in Se,t. Thus Bj,d,τ = m i=1 aijBi,e,t, j= 1, . . . , n, or BT d,τ = BT e,tA, where BT d,τ = [B1,d,τ , . . ., Bn,d,τ] and BT e,t = [B1,e,t, . . ., Bm,e,t] are row vectors. If f = n j=1 cjBj,d,τ then f = m i=1 biBi,e,t, where b = Ac, c = [c1, . . . , cn]T , b = [b1 . . . , bm]T . (3.23) The case where e = d is called knot insertion and corresponds to h-refinement in the finite element literature. The situation where e > d and t = τ (e−d) is called degree 46 raising or degree elevation and corresponds to what is commonly known as p-refinement or p-enrichment [175, 45, 98]. In the general case where τ (e−d) is a proper subset of t both knot insertion and degree raising occur. When this transformation is carried out with degree raising followed by knot insertion, Hughes [87] introduced the term k-refinement to the isogeometric literature. Although it is possible to do the transformation in opposite order, i.e., a knot insertion followed by a degree raising, as observed in [87] in their discussion of k-refinement, this ordering leads to more coefficients and less smooth functions. There are two algorithms for knot insertion. In Boehm's algorithm one knot at a t