Feature Article
A Responsive Finite Element Method to Aid Interactive Geometric Modeling Nobuyuki Umetani and Kenshi Takayama ■ University of Tokyo Jun Mitani ■ University of Tsukuba Takeo Igarashi ■ University of Tokyo
I
nteractive numerical simulations offer powerful assistance in designing items that satisfy specific physical requirements, as Ivan Sutherland’s SketchPad demonstrated.1 (For further information on interactive design and simulation, see the related sidebar.) However, most numericalsimulation methods, generally called computeraided engineering (CAE) tools, are used to verify designs offline. In particular, they’re used to reject designs that don’t satisfy specific requirements; they aren’t often used to explore better designs. Real-time simulation is emerging, but typical realtime applications include simulating deformation in animation2 and virtual training.3 Real-time numerical simulation isn’t widely used to design physical items. In line with Sutherland’s vision, we’ve integrated a numerical-simulation method into geometric modeling. Our system runs a finite element method (FEM) simulation in real time that responds to dynamic user input during geometric editing (see Figure 1). Standard real-time FEMs for deformation are applied to a single given initial geometry. In contrast, our method continuously updates simulation results in response to user modifications of the initial geometry. Although our current implementation is limited to 2D problems with simplex first-order elements, the basic concept of responsive FEM is independent of dimensionality and element types. Published by the IEEE Computer Society
Real-time, responsive feedback such as our system offers can provide principles for better design. It can also help users find a satisfactory design while avoiding the many trial-and-error experiments necessary with offline simulations. In addition, it’s useful for educational purposes.
Our Algorithm To make the FEM framework responsive—that is, to provide immediate feedback during geometric editing—we maximize the reuse of intermediate computation results and carefully schedule the computation pipeline to provide the best user experience. Our algorithm was mostly built on the basic FEM framework described in the following subsection, but we reuse intermediate computation results to make the standard FEM more responsive.
FEM Background
Computer-aided engineering systems use numerical simulation methods mainly to reject undesirable designs—not to guide users toward better designs. However, integrating a realtime finite element method into interactive geometric modeling can provide user guidance. During interactive editing, real-time feedback from numerical simulation guides users toward improved design without tedious trialand-error iterations.
FEM finds an approximate solution of partial differential equations by spatially discretizing the field. FEM first constructs a mesh inside the boundary geometry. It then solves a linear system Ax = b that’s defined by the node relationships (for nonlinear problems,
0272-1716/11/$26.00 © 2011 IEEE
IEEE Computer Graphics and Applications
43
Feature Article
Related Work in Interactive Design and Simulation
R
design.4 However, such approaches don’t differ significantly from conventional CAD and CAE systems in that the analysis occurs after modeling. In contrast, our system performs the simulation during modeling. Researchers have proposed several systems that let users design physical objects, such as stuffed animals, with the aid of interactive simulation.5 Because these systems focus mainly on the user interface, they use simplified simulation methods.
esearchers have investigated using automatic optimization to design objects that satisfy physical constraints. For example, Jeffrey Smith and his colleagues applied optimization to design truss structures.1 However, automatic methods have many difficulties in practice—such as employing parameter spaces that are too large and the difficulty of explicitly specifying constraints. An interactive approach is advantageous because users employ their own preferences and judgment while considering other less tangible factors, such as aesthetics. Researchers have proposed various ways to make physical simulations interactive. One approach is to use precomputations. Doug James and Dinesh Pai presented an interactive physical simulation of deformable objects by precomputing the matrices of reference boundary value problems.2 Another approach is to approximate nonlinear elasticity to achieve large deformations quickly and stably.3 Both approaches employ user input as an external force in a continuously running simulation. However, in our system, the user directly modifies the initial geometry, and the system reruns the simulation with the new geometry. We aim to use physical simulation to aid design. Some CAD systems provide an embedded simulation tool for design evaluation. Such tools seamlessly switch between the design environment and offline simulation. Mark Masry and Hod Lipson presented a sketch-based 3D modeling interface that can perform FEM evaluation in the early stages of
References 1. J. Smith et al., “Creating Models of Truss Structures with Optimization,” ACM Trans. Graphics, vol. 21, no. 3, 2002, pp. 295–301. 2. D.L. James and D.K. Pai, “ArtDefo: Accurate Real Time Deform able Objects,” Proc. Siggraph, ACM Press, 1999, pp. 65–72. 3. M. Müller et al., “Stable Real-Time Deformations,” Proc. 2002 ACM/Eurographics Symp. Computer Animation, ACM Press, 2002, pp. 49–54. 4. M. Masry and H. Lipson, “A Sketch-Based Interface for Iterative Design and Analysis of 3D Objects,” Proc. 2nd Eurographics Workshop Sketch-Based Interfaces, Eurographics Assoc., 2005, pp. 109–118. 5. Y. Mori and T. Igarashi, “Plushie: An Interactive Design System for Plush Toys,” ACM Trans. Graphics, vol. 26, no. 3, 2007, article 45.
Design with responsive FEM
Actual manufacture
Interactive shape design Real-time analysis
724 Hz Musical tone
F
839 Hz F#
G
G#
A
989 Hz A#
B
Fish-shaped metallophone
C
Figure 1. Our finite element method (FEM) concept applied to metallophone design. Real-time FEM analysis occurs during interactive geometric design sessions, letting users create real-world objects with physically desirable properties.
we must solve such linear systems iteratively). Because matrix A is sparse, it’s compactly represented by the combination of the nonzero pattern Ap (which represents the location of nonzero elements) and the value list Av (which represents the values at the nonzero elements). Iterative methods are commonly used to solve sparse linear systems; preconditioners often improve their performance. 44
September/October 2011
So, a preconditioner B approximates the inverse of A, which isn’t necessarily sparse. B is usually represented as a sparse matrix, with its nonzero pattern Bp and nonzero value list Bv. These data (Av, Ap, Bv, and Bp) must be constructed before the linear system can be solved. The construction of the mesh and the matrices can be considered a precomputation. When solving
Table 1. Multilevel reuse of computations. User operation Idle
Dragging Relocation
Reconnection
Reconstruction
Computation Coefficient matrix A Preconditioner matrix B
Value list Av
Reusable
Must be recomputed
Must be recomputed
Must be recomputed
Nonzero pattern Ap
Reusable
Reusable
Must be recomputed
Must be recomputed
Value list Bv
Reusable
Reusable
Reusable
Must be recomputed
Nonzero pattern Bp
Reusable
Reusable
Reusable
Must be recomputed
a linear problem, FEM runs the entire precomputation only once; it finds a solution without changing the matrices and reuses them multiple times to solve a time-varying problem. In contrast, to solve a nonlinear problem, FEM solves the problem iteratively by updating Av and Bv each time. Traditional FEM frameworks reconstruct mesh and matrix computations all at once for every geometry change. When the user applies the same analysis to even a slightly modified geometry, such frameworks discard the results of all precomputations and construct the mesh and matrices from scratch. This is a waste of time because most computations are redundant and can be reused.
Multilevel Reuse In our system, the user modifies the boundary geometry by dragging vertices, edges, or regions, and the system continuously runs FEM analysis on the domain. For structural analysis, users modify the rest shape, not the deformed shape that emerges from the simulation. The challenge is to provide immediate feedback to the user while maintaining a certain level of accuracy. Making a system responsive isn’t the same as simply making it fast. To maximize the trade-off between speed and accuracy, we must be careful when distributing computation that corresponds to the degree of data change the user caused. We achieve this trade-off by reusing intermediate computation results instead of recomputing everything each time the user modifies the boundary geometry. To obtain accurate results, we reuse most of the previous intermediate computation results when the geometric modification is small. If accumulated geometric modification becomes too large, we stop reusing results and run the costly computation to maintain accuracy. We divide the computation into multiple stages and choose the amount of recomputation appropriate for the situation. Because the main FEM computation is divided into mesh construction and matrix computation, we must recompute the mesh and matrices
when modifying the boundary geometry. We’ve defined three recomputation levels: relocation, reconnection, and reconstruction; we choose the appropriate one to balance speed and accuracy (see Table 1). If the mesh stays the same, reusing computations will result in the same final solution as regular FEM analysis because we solve the same coefficient matrix iteratively using the same convergence criterion. Continuously updating a mesh during simulation is already used to solve problems involving geometry changes, such as fluid-structure interaction based on arbitrary Lagrangian-Eulerian methods. However, such off line methods don’t selectively apply different update procedures based on user input, the way our method does. When the boundary geometry modification is small, we change only the mesh nodes’ position (relocation); this doesn’t change the mesh topology. So, we need to update only Av while reusing Ap, Bv, and Bp. Because the nodes move only slightly, the solution (the field values at the nodes) doesn’t change much. This means we can reuse the FEM solution as an initial guess in running an iterative solver, which is faster than starting from an arbitrary guess. When the geometry modification is larger, node relocation doesn’t eliminate mesh distortions. So, we change the mesh’s topology locally to improve mesh quality (reconnection). In this case, we must update Ap and Av. We can still reuse Bv and Bp because we’re not adding or deleting any nodes and because the existing nodes move only slightly. Other methods have reused the preconditioner, but usually to solve nonlinear problems. As with relocation, we can reuse the FEM solution as the initial guess in the iterative solver. When the geometry modification is significantly large, even reconnection isn’t sufficient. In this case, we stop the incremental mesh updates and reconstruct the entire mesh from scratch. This might sound radical. However, when the distortion has accumulated or the geometry has changed IEEE Computer Graphics and Applications
45
Feature Article
Figure 2. A modeling sequence used for performance measurement. The user continuously dragged the hole from left to right. The mesh consisted of 1,952 elements. The mesh is changed mainly with relocation and reconnection to avoid costly computation of the reconstruction.
too much, global reconstruction is often faster and yields a better mesh than local optimization with node insertion and deletion. In this case, we recompute Av, Ap, Bv, and Bp. In addition, we can’t reuse the previous solution because new nodes have completely replaced the old ones. So, we must start with a new initial guess. We reuse FEM data to maximize analysis responsiveness by considering the cost required for each recomputation level. We rely mostly on the lightweight relocation recomputation and perform the expensive reconstruction recomputation only when necessary. The basic concept we’ve described applies to linear and nonlinear problems; however, the details differ slightly in nonlinear cases. Specifically, we can’t reuse Av and Bv when solving nonlinear problems because they must be updated each time. However, we can still reuse Ap and Bp, which significantly improves performance.
Implementation Details Multilevel reuse first determines what part of the mesh to reuse and then uses this to decide what part of the matrix computation to reuse. When the user edits the boundary geometry, the system first relocates the nodes to minimize mesh distortion. In cases where the system detects an inverted element after the relocation, it pushes the nodes back to the previous positions and applies reconstruction. If reconnection doesn’t occur, this means that reconnection doesn’t improve the mesh quality, and the system reconstructs the mesh. In cases where the system detects no inverted element after relocation, it checks for distorted elements and, if they exist, applies reconnection. We use a simple mass-spring system for node relocation. The spring’s rest length is zero, and we solve equilibrium iteratively. Because the nodes move only slightly each time, the computation converges quickly. To avoid oscillation, an even number of iterations is desirable; we currently perform two. The distortion metric is based on the ratio of an inscribed circle and the maximal edge length. We apply edge swapping for reconnection; we swap any edge that violates the Delaunay condition. Mesh reconstruction is based on Delaunay triangulation and local optimization. 46
September/October 2011
The conjugate-gradient method is used for solving symmetric matrices, whereas the Bi-CGSTAB (biconjugate gradient stabilized) method is used for asymmetric matrices. We improve these iterative methods’ convergence by using a preconditioner based on incomplete LU factorization with level of fill-in. Incomplete LU factorization computes a sparse matrix B that approximates the inverse of a sparse matrix A (the exact inverse of A generally isn’t sparse). The method takes an integer called the level of fill-in as a parameter specifying the approximation’s level. This level affects both the convergence’s improvement and the factorization’s cost. The higher the level, the more closely B approximates the inverse of A, leading to faster convergence but requiring more computation for the factorization. A preconditioner with a high-level of fill-in benefits more from multilevel reuse, which greatly decreases the number of preconditioner recomputations. However, the best level of fill-in depends heavily on the target problem; we experimentally chose an appropriate level for each application. (For example, we used a level of three for the vibration analysis and cantilever deformation and a level of zero for the fluid and thermal-fluid examples.)
Algorithm Performance The modeling sequence in Figure 2 demonstrates multilevel reuse’s effectiveness. In that sequence, relocation occurred 117 times, reconnection occurred 39 times, and reconstruction occurred only three times. Figure 3 shows the cost for each level. We measured the cost using the FEM problem we describe later in our vibration analysis example, which we tested with a 2.5-GHz CPU and 2.0 Gbytes of RAM.
Application Examples Here, we apply our responsive FEM framework to typical 2D design problems. In each case, the user interactively manipulates an object’s shape within a certain physical constraint, and the system returns the analysis result in real time. We believe these responsive simulations will be useful for both early exploration of new design problems and refinement of designs that are almost finished,
FEM
Mesh
2.1
+
3.1
+
2.2
Reconnection
2.1
+
3.1
+
2.2
+ 1.4
Reconstruction
2.1
+
3.1
+
2.2
+ 1.4 +
4.4
Av
Ap
Bv
Relocation
Misc.
Solve
+
+
+
4.4
+
1.4
2.0
3.8
= 8.8
= 10.8
= 21.4
Bp
Figure 3. The cost (in ms) for each level of recomputation. Reconstruction recomputations are more than twice as expensive as relocation recomputations, owing mainly to the cost of constructing the preconditioner’s value list (Bv) and nonzero pattern (Bp). Av and Ap are the coefficient matrix’s value list and nonzero pattern.
(a)
(b) Figure 4. Vibration analysis. (a) During design, the user moves a window. (b) During analysis, the system displays the deformation caused by ground movement. Resonance occurs when the user moves the top right window toward the bottom, leading to a large destructive deformation. We exaggerated the displayed deformation for the purpose of visualization.
and we present current implementations primarily as a proof of concept. They might not be useful for real-world applications; building practical applications based on these examples remains a subject for future research. Still, we believe that the current implementations are useful for some enduser design problems and for teaching the general principles of physical phenomena.
Vibration Analysis of a Structural Object In this example, a structural object is fixed to the ground, which shakes constantly at a certain frequency, causing the entire structure to deform (see Figure 4). Resonance appears when the user manipulates the object into a specific shape—one that would be predictable only through simulation. With more sophistication, this application might help in designing buildings that avoid collapse due to the resonance caused by earthquakes or wind. This analysis is based on a nonstationary 2D linear solid without gravity:
rü = ∇ ⋅ s s = l(tre)I + 2me, where r is the density, u is the displacement, s is the Cauchy stress tensor, e is the linearized strain tensor, and l and m are the elastic Lamé coefficients. We use I and tr as they are commonly used—indicating unit tensor and trace operation of the matrix, respectively. The time integration is based on the Newmark-b method.
Cantilever Deformation As Figure 5 shows, the leftmost part of a horizontal cantilever is fixed to a vertical wall, and the remainder is free. Gravity causes the whole cantilever to deform. The user edits the initial cantilever shape trying to fit the deformed shape to the target shape. This application might benefit airfoil designs in which designers are more interested in the shape’s hydrodynamic performance after gravity and wind pressure cause deformation than they are in the original, undeformed shape. Automatic IEEE Computer Graphics and Applications
47
Feature Article
… (a)
(b)
Figure 5. Cantilever deformation. The user (a) continuously manipulates the original, undeformed shape to (b) fit the goal shape (shown in red lines) after gravity causes deformation. Because the deformed shape changes instantly according to user manipulation, the user can design a proper cantilever shape.
(a)
(b)
(c)
Figure 6. Fluid around an object. The velocity field appears as line segments; the pressure appears as color contours. As the user manipulates the object shape, various phenomena are visible. (a) Boundary layer separation. (b) A Karman vortex street. (c) The desirable object shape obtained after shape editing.
optimization is often used for this kind of problem when the goal shape is clearly defined. However, users might have only vague ideas about the goal shape and wish to try various designs before deciding on one. In such cases, continuous feedback is very useful. Also, the designed shape can serve as an initial guess for automatic optimization. For this application, we solve the St. VenantKirchhoff material equation: S = l(trE)I + 2mE, where S is the second Piola-Kirchhoff stress tensor and E is the Green-Lagrange strain tensor.
Fluid around an Object Here, an object is inside an air-filled space, and wind of a certain velocity blows constantly from left to right, creating a complex flow around the object. Various kinds of phenomena are visible depending on the object shape the user manipulates. For example, boundary layer separation (see Figure 6a) can cause a stall, and a Karman vortex street (see Figure 6b) leads to potentially destructive oscillation. Figure 6c shows an example of appropriate design. This has potential utility in the design of any object constantly exposed to strong flow, including airfoils, car bodies, door mirrors, and air ducts. 48
September/October 2011
For this application, we solve an incompressible Newtonian flow with the SUPG-PSPG (streamline upwind Petrov-Galerkin and pressure-stabilizing Petrov-Galerkin) algorithm:4 r(Dv/Dt) = -∇p + m∇2v ∇ ⋅ v = 0, where v is the velocity, r is the density, p is the pressure, and m is the viscous modulus. We use D as it’s commonly used—indicating material derivative. We use the implicit method for time integration.
Thermal Fluid in an Object In Figure 7, some kind of fluid (in this case, water) fills an object (a teapot) whose bottom is heated while other boundaries are constantly cooled. We observe how the thermal fluid’s complex, nonstationary behavior changes according to the object shape that the user manipulates. Besides the design of a heat-efficient teapot, we expect this application could be useful for various problems related to thermal-fluid phenomena—such as the layout of room air conditioners or the design of a computer case. We solve the Navier-Stokes equation with buoyancy proportional to the temperature, which we compute through a convection-diffusion equation:
Figure 7. Thermal fluid in an object. The temperature appears as a color contour; blue and red correspond to low and high temperatures. The user can determine the shape of the teapot while viewing the heat transport inside.
Table 2. Performance of our application examples tested with a 2.5-GHz CPU and 2.0 Gbytes of RAM. Application
Linear
Stationary
FPS (reuse)
FPS (no reuse)
Vibration
Yes
No
1,962
105.0
41.8
Cantilever
No
Yes
990
37.6
7.1
Fluid
No
No
1,971
30.0
18.2
Thermal fluid
No
No
1,938
21.3
14.3
r(Dv/Dt) = -∇p + m∇2v + rgb(T - T0) ∇⋅v=0 DT/Dt = a∇2T, where T is the temperature, T0 is the reference temperature, a is the thermal diffusivity, b is the volumetric thermal expansion ratio, and g is the acceleration of gravity. We again use the implicit method for time integration.
Application Performance Table 2 summarizes these applications’ performance. The frame rate depends on user manipulation (for example, the rate decreases when the user quickly makes a large shape modification). So, we averaged the frame rates for moderate-speed user manipulations similar to the manipulation in the “Algorithm Performance” section. The table shows that with data reuse, our applications ran at a high frame rate, sufficient for interactive modeling. We observed that our reuse scheme was especially effective for linear or stationary problems.
Informal User Studies To verify responsive FEM’s utility, we performed two informal user studies: one with a professional artist and one with nonprofessional students.
Metallophone Design To show that our current implementation is already useful as a practical tool for designing real-world objects, we used it for the design and creation of an artistically shaped metallophone. Metallophones are usually rectangular because that shape is suit
No. of elements
able for analytically predicting the instrument’s tone. However, researchers have used computer simulation to simulate the sound from arbitrarily shaped objects. For example, Jeffrey Chadwick and his colleagues proposed a computationally efficient framework for synthesizing realistic sound from nonlinear, thin-shell vibrations.5 However, we believe that designing a metallophone with a desired artistic shape and tone is possible only with responsive FEM, because this kind of highly constrained modeling naturally demands the tight integration of design and analysis. We model the metallophone as a 3D plate with a predefined thickness extruded from the 2D shape that the user designs. We compute its tone, or frequency, through eigenvalue analysis. We assume the lowest nonzero eigenfrequency is the output tone because higher eigenfrequencies usually attenuate quickly. Along with computing the eigenfrequency, we compute the eigenmode (which tells how the metallophone oscillates), which is helpful in determining which positions on the metallophone should be fixed. For details on these computations, see the sidebar “Computing the Metallophone’s Tone and Mode.” Simulation parameters must be calibrated for each specific material. Figure 8 shows our metallophone design software. Users edit the shape in 2D in the left window while checking the oscillation’s eigenmode in 3D in the right window. During design, the software updates the tone in real time, with visual feedback in the status bar and aural feedback using beeps. We asked a professional artist to design a metallophone using our software. We put no time limit IEEE Computer Graphics and Applications
49
Feature Article
Computing the Metallophone’s Tone and Mode
A
ssuming that the metallophone floats in nongravity space without any external forces and fixed boundaries, we use finite element method (FEM) discretization to obtain this equation: Mu + Ku = 0,(A)
where u is the nodal displacement vector, M is the lumped mass matrix, and K is the positive, semidefinite stiffness matrix. We split u into the product of the spatially varying amplitude f and the harmonic oscillation at the angle rate w, as u(x, t) = f(x)eiwt, where t indicates time and x indicates node positions. We then substitute this product into Equation A, which yields this generalized eigenvalue problem: Kφ = λMφ ,(B) where the eigenfrequency is f = ω (2π ) = λ (2π ) . We calculate the smallest nonzero eigenvalue l and its corresponding eigenvector f to obtain the metallophone’s tone and mode. We perform Cholesky factorization on the lumped mass matrix as M = LLT and multiply L –1 and L –T to the left and right sides of Equation A. Here, L stands for a lower triangular matrix. We use T as it’s typically used—indicating the transpose operation of the matrix. Thus we obtain this eigenvalue problem: Ay = ly,
Figure 8. Metallophone design software. Users edit the shape in the left window; the right window shows the analyzed eigenmode in 3D. The software updates the tone in real time with both aural and visual feedback.
where A = L –1KL –T and y = LTf. We solve this using inverse iteration. We modify the original iteration procedure by adding a step to remove the component of all zero eigenvectors of A from the current solution. That way, we obtain the smallest nonzero eigenvalue and its corresponding eigenvector. Thanks to our problem setting, we already know the zero eigenvectors of K as f0i (i = 1 … 6): the translations along the three coordinate axes and the rotations around them. We apply the modified Gram-Schmidt orthogonalization to LT f0i (i = 1 … 6) to obtain the orthonormal basis i vectors y0 (i = 1 … 6) that span the kernel of A. We also define a projection P that maps a vector v to the compliment space of the kernel of A as P (v ) = v − y0i (y0i ⋅ v ). In each iteration step, we apply this projection to the solution vector and normalize it. We add a small positive number e to the diagonals of A to improve the condition number. Once we’ve computed the shifted eigenvalue l1′ and its corresponding eigenvector y1 of A, we finally obtain the nonzero smallest eigenvalue λ1 = λ1′ − ε and the eigenvector f1= L –Ty1 in Equation B. Reusing the previous time step’s solution significantly improves the inverse iteration’s convergence. The simulation accuracy is sensitive to the mesh density because we use linear tetrahedral elements to represent the bend of a 3D thin plate. We can alleviate this problem, called shear locking, by preventing the excessive distortion of tetrahedral elements. We keep the longest edge shorter than twice the shortest edge.
∑
2D view (shape design)
3D view (eigenmode)
Scale and frequency
50
September/October 2011
C
D
E
F
G
A
B
(a)
(b) Figure 9. Designing a metallaphone. (a) The 2D shapes the artist designed. (b) The analyzed eigenmode of the designed pieces in 3D. For visualization purposes, the displayed eigenmode is considerably exaggerated. Using our system, the artist successfully designed complicated shapes that have specific tones.
on the design and provided instruction on the software. Figure 9 shows designed shapes corresponding to musical notes ranging from C (523 Hz) to B (987 Hz). We used a wire-electrical-discharge machine to cut the shapes from 4-mm aluminum plate; we fixed the shapes onto a wooden board according to the analyzed eigenmodes (see the right side of Figure 1). The top three rows of Table 3 show that for most pieces, the target, simulated, and actual frequencies were quite close. To further improve the metallophone’s quality, we manually adjusted the actual pieces’ tones by trimming their edges (we didn’t have to manually adjust the F piece’s tone). We used our software during the manual adjustment because it can predict tone changes caused by trimming. The bottom row of Table 3 shows the adjustments. We found that the actual metallophone produced sounds of sufficient quality for a hobbyist. After the design process, we interviewed the artist to obtain subjective comments. The artist reported that the design took roughly five hours, with most of the time devoted to the C and D pieces. Because these lower tones required larger areas, the analysis response time was much slower than for the other pieces. The artist had difficulty balancing the overall shape among the pieces while keeping their tones true; this was a huge design constraint. Additionally, sometimes a small shape modification resulted in a large change in tone, which demanded a highly responsive analysis. The artist commented that designing such an artistically shaped metallophone without using responsive FEM would be almost impossible. This prototype system is still a proof of concept, and many problems remain to be solved in designing a more realistic metallophone. For example, whereas we compute only fundamental frequencies, obtaining nicely tuned harmonics from secondary modes is necessary to produce a comfortable sound. Addressing these issues remains the subject of future research.
Table 3. The metallophone’s frequencies for each note, showing the accuracy of our analysis. Note Type of frequency Targeted
C
D
E
F
G
A
B
523
587
659
698
783
880
987
Simulated
525
588
661
699
786
880
989
Measured
506
604
621
698
787
860
993
Adjusted
523
587
659
698
783
881
987
Bridge Design The other user study concerned designing a bridge. The aim was to show that responsive FEM provides better support than traditional, nonresponsive FEM for nonprofessionals designing objects with physically desirable properties. The task. We asked the users to design the 2D shape of a bridge to span a certain gap and support a certain weight at its center, using as small an area as possible. That is, the goal was to design a strong bridge using the least amount of material. Figure 10 shows the overview of our bridge design system. The system provides a set of tools—such as curve editing using virtual pushpins,6 curve smoothing, and hole creation. It always displays the bridge’s area. The software tests the bridge’s strength through FEM analysis using a physical model based on equivalent stress. The system displays the amount of stress applied to each region as color contours (blue and red correspond to low and high stress) and judges whether the bridge passes the test. FEM analysis is available in two modes: responsive and nonresponsive. In responsive FEM, the analysis always occurs during user interaction (for example, mouse dragging), and the result updates in real time. In nonresponsive FEM, the analysis occurs only after the user completes the design and presses a toolbar button. The analysis result immediately disappears when the user changes IEEE Computer Graphics and Applications
51
Feature Article
Figure 10. Testing a bridge design’s structural soundness using FEM analysis. Blue and red contours correspond to low and high stress. Users were asked to design a bridge that wouldn’t break down under a weight located in its center, using the least amount of material. Responsive FEM
Nonresponsive FEM
Area
30 20 10
1
2 3 Participant
Avg.
4
5 6 Participant
Avg.
Subject
Figure 11. The bridge design study’s results. Subjects who used responsive FEM generally achieved better designs, with passing bridges encompassing smaller areas than those designed with nonresponsive FEM.
the design. The nonresponsive mode simulates the typical use of FEM systems, in which design and analysis are completely separate. In such systems, switching between these two processes can involve much tedious work, such as exporting and importing files and setting FEM parameters. Note that the nonresponsive FEM we used is already much more efficient than current commercial FEM tools. The experimental setup. Six university students majoring in art and design participated; none was familiar with FEM techniques or material mechanics. We split them into two groups: group A started with responsive FEM; group B started with nonresponsive FEM. We performed separate experiments for each group. Each group underwent a 15-minute lesson on using the software and then a 30-minute main design session. Afterward, we asked the participants in each group to try the other FEM mode. We then collected their subjective feedback. 52
September/October 2011
Results. Figure 11 shows the smallest-area bridge that passed the strength test for each participant. Participants who used responsive FEM generally achieved better results. The most common feedback was that responsive FEM is useful when the user wants to see how a small adjustment affects the analysis. The participants also mentioned that the analysis results displayed during responsive FEM were sometimes too conspicuous, making large design changes difficult because users felt encouraged to optimize results locally with small edits. Some participants even claimed that shape design without responsive FEM might be more appropriate for initial design exploration. This simple study obviously doesn’t “prove” anything, but we believe it suggests responsive FEM’s potential as a design aid for nonprofessionals. We plan to perform more formal user studies.
Limitations We must extend our techniques to 3D to make them useful for many practical real-world problems. Solving linear systems in 3D is straightforward; our main challenge will be working the continuous mesh update in 3D. As François Labelle and Jonathan Shewchuk noted, methods for improving tetrahedral mesh quality by continuously moving nodes and changing connectivity have yet to guarantee sufficiently accurate simulations.7 Another limitation is that currently we can use only first-order elements. Higher-order elements are much more desirable for such problems as bending thin-walled structures. However, using such elements in our approach is problematic because reconnecting the mesh will most likely change relationships between nodes and prevent matrix reuse. Additionally, we’ve yet to try nonsimplex elements (for example, quadrilaterals in 2D and hexahedrons in 3D). Such elements might—depending on the problem—be more appropriate than simplex elements (for example, triangles in 2D and tetrahedrons in 3D). However, reconnecting nonsimplex meshes without adding and deleting points is difficult. Our approach isn’t applicable to history-dependent problems (such as plasticity processing). In such problems, the solutions must be computed sequentially from the initial state, and our data reuse scheme is inappropriate for that purpose.
O
ur current applications already show potential utility in areas as diverse as aerodynamics, air conditioning layout, and airfoil design. Additional research could make our method even
more globally beneficial. To start with, we intend to test reuse of data other than the preconditioner matrix; these tests might include node reordering, which would improve analysis responsiveness. We could also apply our research to create a system that uses simulation results to more actively guide design. For example, the system could help users design shapes that satisfy certain constraints (for example, stress limits in certain areas) by providing suggestions and instructions when the user makes a change that doesn’t satisfy the constraints. We could also design a system in which the user interactively controls simulation accuracy. We use fixed criteria for the trade-off between speed and accuracy, but users might want more explicit control over this trade-off (for example, during design refinement). We also assume homogeneous mesh density, but it would be useful for users to locally control simulation accuracy by manipulating the local mesh density. In this way, they could examine analysis results more closely in specific regions of interest; existing automatic error estimation methods make this difficult.
References 1. I.E. Sutherland, “Sketch Pad: A Man-Machine Graphical Communication System,” Proc. SHARE Design Automation Workshop, ACM Press, 1964, pp. 6.329–6.346. 2. J. Mezger et al., “Interactive Physically Based Shape Editing,” Proc. 2008 ACM Symp. Solid and Physical Modeling, ACM Press, 2008, pp. 79–89. 3. N. Chentanez et al., “Interactive Simulation of Surgical Needle Insertion and Steering,” ACM Trans. Graphics, vol. 28, no. 3, 2009, article 88. 4. T.E. Tezduyar et al., “Incompressible Flow Compu tations with Stabilized Bilinear and Linear EqualOrder-Interpolation Velocity-Pressure Elements,” Computer Methods in Applied Mechanics and Eng., vol. 95, no. 2, 1992, pp. 221–242. 5. J.N. Chadwick, S.S. An, and D.L. James, “Harmonic Shells: A Practical Nonlinear Sound Model for NearRigid Thin Shells,” ACM Trans. Graphics, vol. 28, no. 5, 2009, article 119. 6. T. Igarashi, T. Moscovich, and J.F. Hughes, “AsRigid-as-Possible Shape Manipulation,” ACM Trans. Graphics, vol. 24, no. 3, 2005, pp. 1134–1141. 7. F. Labelle and J.R. Shewchuk, “Isosurface Stuffing: Fast Tetrahedral Meshes with Good Dihedral Angles,” ACM Trans. Graphics, vol. 26, no. 3, 2007, article 57.
computer simulation and interactive design. Umetani has an MS in frontier science from the University of Tokyo. Contact him at [email protected]. Kenshi Takayama is a PhD student in computer science at the University of Tokyo. His research interests include computer graphics, volumetric modeling, user interfaces, texture synthesis, geometric modeling, and image editing. Takayama has an MS in computer science from the University of Tokyo. Contact him at [email protected]. Jun Mitani is an associate professor of computer science at the University of Tsukuba. His research interests include computer graphics, geometric modeling, and origami. He’s a member of the Japan Society of Graphics Science and the Society for Art and Science of Japan. Mitani has a PhD in precision engineering from the University of Tokyo. Contact him at [email protected]. Takeo Igarashi is an associate professor of computer science at the University of Tokyo. He’s also the director of the Japan Science and Technology Agency ERATO (Exploratory Research for Advanced Technology) Igarashi Design Interface Project. His research interests include user interfaces and interactive computer graphics. Igarashi has a PhD in information engineering from the University of Tokyo. Contact him at [email protected].
Call rticles for A IEEE seek s re
a
uting p m o sive C he latest peerPerva on t
le, u ccessib
d de viewe
u
ou s biquit
s e f ul
velop
pap e
m e nt
co m p
s
uting
rs
vasiv in per
. To
e, mo
clud pics in
ra s t r re inf
u c tur
bile, a
e har
nd
d war
e , rea
e
l-wor
ld
wa u te r , s of t co m p ology man u h , ding n rac tio s , inclu d inte ration n e a id g s s co n sensin . s y s te m rivac y , an d an d p n io t c urit y, c e s , intera y labilit t , sca ymen deplo es: or.htm de lin /auth r gui vasive r o e h /p t Au rg /mc uter.o .comp www
te chn
er Furth ils: det a p e r va co m p
sive@
uter.o
rg
. www ter. u c om p ive er vas or g / p
Nobuyuki Umetani is a PhD student in computer science at the University of Tokyo. His research interests include
IEEE Computer Graphics and Applications
53
A Responsive Finite Element Method to Aid Interactive Geometric Modeling Nobuyuki Umetani and Kenshi Takayama ■ University of Tokyo Jun Mitani ■ University of Tsukuba Takeo Igarashi ■ University of Tokyo
I
nteractive numerical simulations offer powerful assistance in designing items that satisfy specific physical requirements, as Ivan Sutherland’s SketchPad demonstrated.1 (For further information on interactive design and simulation, see the related sidebar.) However, most numericalsimulation methods, generally called computeraided engineering (CAE) tools, are used to verify designs offline. In particular, they’re used to reject designs that don’t satisfy specific requirements; they aren’t often used to explore better designs. Real-time simulation is emerging, but typical realtime applications include simulating deformation in animation2 and virtual training.3 Real-time numerical simulation isn’t widely used to design physical items. In line with Sutherland’s vision, we’ve integrated a numerical-simulation method into geometric modeling. Our system runs a finite element method (FEM) simulation in real time that responds to dynamic user input during geometric editing (see Figure 1). Standard real-time FEMs for deformation are applied to a single given initial geometry. In contrast, our method continuously updates simulation results in response to user modifications of the initial geometry. Although our current implementation is limited to 2D problems with simplex first-order elements, the basic concept of responsive FEM is independent of dimensionality and element types. Published by the IEEE Computer Society
Real-time, responsive feedback such as our system offers can provide principles for better design. It can also help users find a satisfactory design while avoiding the many trial-and-error experiments necessary with offline simulations. In addition, it’s useful for educational purposes.
Our Algorithm To make the FEM framework responsive—that is, to provide immediate feedback during geometric editing—we maximize the reuse of intermediate computation results and carefully schedule the computation pipeline to provide the best user experience. Our algorithm was mostly built on the basic FEM framework described in the following subsection, but we reuse intermediate computation results to make the standard FEM more responsive.
FEM Background
Computer-aided engineering systems use numerical simulation methods mainly to reject undesirable designs—not to guide users toward better designs. However, integrating a realtime finite element method into interactive geometric modeling can provide user guidance. During interactive editing, real-time feedback from numerical simulation guides users toward improved design without tedious trialand-error iterations.
FEM finds an approximate solution of partial differential equations by spatially discretizing the field. FEM first constructs a mesh inside the boundary geometry. It then solves a linear system Ax = b that’s defined by the node relationships (for nonlinear problems,
0272-1716/11/$26.00 © 2011 IEEE
IEEE Computer Graphics and Applications
43
Feature Article
Related Work in Interactive Design and Simulation
R
design.4 However, such approaches don’t differ significantly from conventional CAD and CAE systems in that the analysis occurs after modeling. In contrast, our system performs the simulation during modeling. Researchers have proposed several systems that let users design physical objects, such as stuffed animals, with the aid of interactive simulation.5 Because these systems focus mainly on the user interface, they use simplified simulation methods.
esearchers have investigated using automatic optimization to design objects that satisfy physical constraints. For example, Jeffrey Smith and his colleagues applied optimization to design truss structures.1 However, automatic methods have many difficulties in practice—such as employing parameter spaces that are too large and the difficulty of explicitly specifying constraints. An interactive approach is advantageous because users employ their own preferences and judgment while considering other less tangible factors, such as aesthetics. Researchers have proposed various ways to make physical simulations interactive. One approach is to use precomputations. Doug James and Dinesh Pai presented an interactive physical simulation of deformable objects by precomputing the matrices of reference boundary value problems.2 Another approach is to approximate nonlinear elasticity to achieve large deformations quickly and stably.3 Both approaches employ user input as an external force in a continuously running simulation. However, in our system, the user directly modifies the initial geometry, and the system reruns the simulation with the new geometry. We aim to use physical simulation to aid design. Some CAD systems provide an embedded simulation tool for design evaluation. Such tools seamlessly switch between the design environment and offline simulation. Mark Masry and Hod Lipson presented a sketch-based 3D modeling interface that can perform FEM evaluation in the early stages of
References 1. J. Smith et al., “Creating Models of Truss Structures with Optimization,” ACM Trans. Graphics, vol. 21, no. 3, 2002, pp. 295–301. 2. D.L. James and D.K. Pai, “ArtDefo: Accurate Real Time Deform able Objects,” Proc. Siggraph, ACM Press, 1999, pp. 65–72. 3. M. Müller et al., “Stable Real-Time Deformations,” Proc. 2002 ACM/Eurographics Symp. Computer Animation, ACM Press, 2002, pp. 49–54. 4. M. Masry and H. Lipson, “A Sketch-Based Interface for Iterative Design and Analysis of 3D Objects,” Proc. 2nd Eurographics Workshop Sketch-Based Interfaces, Eurographics Assoc., 2005, pp. 109–118. 5. Y. Mori and T. Igarashi, “Plushie: An Interactive Design System for Plush Toys,” ACM Trans. Graphics, vol. 26, no. 3, 2007, article 45.
Design with responsive FEM
Actual manufacture
Interactive shape design Real-time analysis
724 Hz Musical tone
F
839 Hz F#
G
G#
A
989 Hz A#
B
Fish-shaped metallophone
C
Figure 1. Our finite element method (FEM) concept applied to metallophone design. Real-time FEM analysis occurs during interactive geometric design sessions, letting users create real-world objects with physically desirable properties.
we must solve such linear systems iteratively). Because matrix A is sparse, it’s compactly represented by the combination of the nonzero pattern Ap (which represents the location of nonzero elements) and the value list Av (which represents the values at the nonzero elements). Iterative methods are commonly used to solve sparse linear systems; preconditioners often improve their performance. 44
September/October 2011
So, a preconditioner B approximates the inverse of A, which isn’t necessarily sparse. B is usually represented as a sparse matrix, with its nonzero pattern Bp and nonzero value list Bv. These data (Av, Ap, Bv, and Bp) must be constructed before the linear system can be solved. The construction of the mesh and the matrices can be considered a precomputation. When solving
Table 1. Multilevel reuse of computations. User operation Idle
Dragging Relocation
Reconnection
Reconstruction
Computation Coefficient matrix A Preconditioner matrix B
Value list Av
Reusable
Must be recomputed
Must be recomputed
Must be recomputed
Nonzero pattern Ap
Reusable
Reusable
Must be recomputed
Must be recomputed
Value list Bv
Reusable
Reusable
Reusable
Must be recomputed
Nonzero pattern Bp
Reusable
Reusable
Reusable
Must be recomputed
a linear problem, FEM runs the entire precomputation only once; it finds a solution without changing the matrices and reuses them multiple times to solve a time-varying problem. In contrast, to solve a nonlinear problem, FEM solves the problem iteratively by updating Av and Bv each time. Traditional FEM frameworks reconstruct mesh and matrix computations all at once for every geometry change. When the user applies the same analysis to even a slightly modified geometry, such frameworks discard the results of all precomputations and construct the mesh and matrices from scratch. This is a waste of time because most computations are redundant and can be reused.
Multilevel Reuse In our system, the user modifies the boundary geometry by dragging vertices, edges, or regions, and the system continuously runs FEM analysis on the domain. For structural analysis, users modify the rest shape, not the deformed shape that emerges from the simulation. The challenge is to provide immediate feedback to the user while maintaining a certain level of accuracy. Making a system responsive isn’t the same as simply making it fast. To maximize the trade-off between speed and accuracy, we must be careful when distributing computation that corresponds to the degree of data change the user caused. We achieve this trade-off by reusing intermediate computation results instead of recomputing everything each time the user modifies the boundary geometry. To obtain accurate results, we reuse most of the previous intermediate computation results when the geometric modification is small. If accumulated geometric modification becomes too large, we stop reusing results and run the costly computation to maintain accuracy. We divide the computation into multiple stages and choose the amount of recomputation appropriate for the situation. Because the main FEM computation is divided into mesh construction and matrix computation, we must recompute the mesh and matrices
when modifying the boundary geometry. We’ve defined three recomputation levels: relocation, reconnection, and reconstruction; we choose the appropriate one to balance speed and accuracy (see Table 1). If the mesh stays the same, reusing computations will result in the same final solution as regular FEM analysis because we solve the same coefficient matrix iteratively using the same convergence criterion. Continuously updating a mesh during simulation is already used to solve problems involving geometry changes, such as fluid-structure interaction based on arbitrary Lagrangian-Eulerian methods. However, such off line methods don’t selectively apply different update procedures based on user input, the way our method does. When the boundary geometry modification is small, we change only the mesh nodes’ position (relocation); this doesn’t change the mesh topology. So, we need to update only Av while reusing Ap, Bv, and Bp. Because the nodes move only slightly, the solution (the field values at the nodes) doesn’t change much. This means we can reuse the FEM solution as an initial guess in running an iterative solver, which is faster than starting from an arbitrary guess. When the geometry modification is larger, node relocation doesn’t eliminate mesh distortions. So, we change the mesh’s topology locally to improve mesh quality (reconnection). In this case, we must update Ap and Av. We can still reuse Bv and Bp because we’re not adding or deleting any nodes and because the existing nodes move only slightly. Other methods have reused the preconditioner, but usually to solve nonlinear problems. As with relocation, we can reuse the FEM solution as the initial guess in the iterative solver. When the geometry modification is significantly large, even reconnection isn’t sufficient. In this case, we stop the incremental mesh updates and reconstruct the entire mesh from scratch. This might sound radical. However, when the distortion has accumulated or the geometry has changed IEEE Computer Graphics and Applications
45
Feature Article
Figure 2. A modeling sequence used for performance measurement. The user continuously dragged the hole from left to right. The mesh consisted of 1,952 elements. The mesh is changed mainly with relocation and reconnection to avoid costly computation of the reconstruction.
too much, global reconstruction is often faster and yields a better mesh than local optimization with node insertion and deletion. In this case, we recompute Av, Ap, Bv, and Bp. In addition, we can’t reuse the previous solution because new nodes have completely replaced the old ones. So, we must start with a new initial guess. We reuse FEM data to maximize analysis responsiveness by considering the cost required for each recomputation level. We rely mostly on the lightweight relocation recomputation and perform the expensive reconstruction recomputation only when necessary. The basic concept we’ve described applies to linear and nonlinear problems; however, the details differ slightly in nonlinear cases. Specifically, we can’t reuse Av and Bv when solving nonlinear problems because they must be updated each time. However, we can still reuse Ap and Bp, which significantly improves performance.
Implementation Details Multilevel reuse first determines what part of the mesh to reuse and then uses this to decide what part of the matrix computation to reuse. When the user edits the boundary geometry, the system first relocates the nodes to minimize mesh distortion. In cases where the system detects an inverted element after the relocation, it pushes the nodes back to the previous positions and applies reconstruction. If reconnection doesn’t occur, this means that reconnection doesn’t improve the mesh quality, and the system reconstructs the mesh. In cases where the system detects no inverted element after relocation, it checks for distorted elements and, if they exist, applies reconnection. We use a simple mass-spring system for node relocation. The spring’s rest length is zero, and we solve equilibrium iteratively. Because the nodes move only slightly each time, the computation converges quickly. To avoid oscillation, an even number of iterations is desirable; we currently perform two. The distortion metric is based on the ratio of an inscribed circle and the maximal edge length. We apply edge swapping for reconnection; we swap any edge that violates the Delaunay condition. Mesh reconstruction is based on Delaunay triangulation and local optimization. 46
September/October 2011
The conjugate-gradient method is used for solving symmetric matrices, whereas the Bi-CGSTAB (biconjugate gradient stabilized) method is used for asymmetric matrices. We improve these iterative methods’ convergence by using a preconditioner based on incomplete LU factorization with level of fill-in. Incomplete LU factorization computes a sparse matrix B that approximates the inverse of a sparse matrix A (the exact inverse of A generally isn’t sparse). The method takes an integer called the level of fill-in as a parameter specifying the approximation’s level. This level affects both the convergence’s improvement and the factorization’s cost. The higher the level, the more closely B approximates the inverse of A, leading to faster convergence but requiring more computation for the factorization. A preconditioner with a high-level of fill-in benefits more from multilevel reuse, which greatly decreases the number of preconditioner recomputations. However, the best level of fill-in depends heavily on the target problem; we experimentally chose an appropriate level for each application. (For example, we used a level of three for the vibration analysis and cantilever deformation and a level of zero for the fluid and thermal-fluid examples.)
Algorithm Performance The modeling sequence in Figure 2 demonstrates multilevel reuse’s effectiveness. In that sequence, relocation occurred 117 times, reconnection occurred 39 times, and reconstruction occurred only three times. Figure 3 shows the cost for each level. We measured the cost using the FEM problem we describe later in our vibration analysis example, which we tested with a 2.5-GHz CPU and 2.0 Gbytes of RAM.
Application Examples Here, we apply our responsive FEM framework to typical 2D design problems. In each case, the user interactively manipulates an object’s shape within a certain physical constraint, and the system returns the analysis result in real time. We believe these responsive simulations will be useful for both early exploration of new design problems and refinement of designs that are almost finished,
FEM
Mesh
2.1
+
3.1
+
2.2
Reconnection
2.1
+
3.1
+
2.2
+ 1.4
Reconstruction
2.1
+
3.1
+
2.2
+ 1.4 +
4.4
Av
Ap
Bv
Relocation
Misc.
Solve
+
+
+
4.4
+
1.4
2.0
3.8
= 8.8
= 10.8
= 21.4
Bp
Figure 3. The cost (in ms) for each level of recomputation. Reconstruction recomputations are more than twice as expensive as relocation recomputations, owing mainly to the cost of constructing the preconditioner’s value list (Bv) and nonzero pattern (Bp). Av and Ap are the coefficient matrix’s value list and nonzero pattern.
(a)
(b) Figure 4. Vibration analysis. (a) During design, the user moves a window. (b) During analysis, the system displays the deformation caused by ground movement. Resonance occurs when the user moves the top right window toward the bottom, leading to a large destructive deformation. We exaggerated the displayed deformation for the purpose of visualization.
and we present current implementations primarily as a proof of concept. They might not be useful for real-world applications; building practical applications based on these examples remains a subject for future research. Still, we believe that the current implementations are useful for some enduser design problems and for teaching the general principles of physical phenomena.
Vibration Analysis of a Structural Object In this example, a structural object is fixed to the ground, which shakes constantly at a certain frequency, causing the entire structure to deform (see Figure 4). Resonance appears when the user manipulates the object into a specific shape—one that would be predictable only through simulation. With more sophistication, this application might help in designing buildings that avoid collapse due to the resonance caused by earthquakes or wind. This analysis is based on a nonstationary 2D linear solid without gravity:
rü = ∇ ⋅ s s = l(tre)I + 2me, where r is the density, u is the displacement, s is the Cauchy stress tensor, e is the linearized strain tensor, and l and m are the elastic Lamé coefficients. We use I and tr as they are commonly used—indicating unit tensor and trace operation of the matrix, respectively. The time integration is based on the Newmark-b method.
Cantilever Deformation As Figure 5 shows, the leftmost part of a horizontal cantilever is fixed to a vertical wall, and the remainder is free. Gravity causes the whole cantilever to deform. The user edits the initial cantilever shape trying to fit the deformed shape to the target shape. This application might benefit airfoil designs in which designers are more interested in the shape’s hydrodynamic performance after gravity and wind pressure cause deformation than they are in the original, undeformed shape. Automatic IEEE Computer Graphics and Applications
47
Feature Article
… (a)
(b)
Figure 5. Cantilever deformation. The user (a) continuously manipulates the original, undeformed shape to (b) fit the goal shape (shown in red lines) after gravity causes deformation. Because the deformed shape changes instantly according to user manipulation, the user can design a proper cantilever shape.
(a)
(b)
(c)
Figure 6. Fluid around an object. The velocity field appears as line segments; the pressure appears as color contours. As the user manipulates the object shape, various phenomena are visible. (a) Boundary layer separation. (b) A Karman vortex street. (c) The desirable object shape obtained after shape editing.
optimization is often used for this kind of problem when the goal shape is clearly defined. However, users might have only vague ideas about the goal shape and wish to try various designs before deciding on one. In such cases, continuous feedback is very useful. Also, the designed shape can serve as an initial guess for automatic optimization. For this application, we solve the St. VenantKirchhoff material equation: S = l(trE)I + 2mE, where S is the second Piola-Kirchhoff stress tensor and E is the Green-Lagrange strain tensor.
Fluid around an Object Here, an object is inside an air-filled space, and wind of a certain velocity blows constantly from left to right, creating a complex flow around the object. Various kinds of phenomena are visible depending on the object shape the user manipulates. For example, boundary layer separation (see Figure 6a) can cause a stall, and a Karman vortex street (see Figure 6b) leads to potentially destructive oscillation. Figure 6c shows an example of appropriate design. This has potential utility in the design of any object constantly exposed to strong flow, including airfoils, car bodies, door mirrors, and air ducts. 48
September/October 2011
For this application, we solve an incompressible Newtonian flow with the SUPG-PSPG (streamline upwind Petrov-Galerkin and pressure-stabilizing Petrov-Galerkin) algorithm:4 r(Dv/Dt) = -∇p + m∇2v ∇ ⋅ v = 0, where v is the velocity, r is the density, p is the pressure, and m is the viscous modulus. We use D as it’s commonly used—indicating material derivative. We use the implicit method for time integration.
Thermal Fluid in an Object In Figure 7, some kind of fluid (in this case, water) fills an object (a teapot) whose bottom is heated while other boundaries are constantly cooled. We observe how the thermal fluid’s complex, nonstationary behavior changes according to the object shape that the user manipulates. Besides the design of a heat-efficient teapot, we expect this application could be useful for various problems related to thermal-fluid phenomena—such as the layout of room air conditioners or the design of a computer case. We solve the Navier-Stokes equation with buoyancy proportional to the temperature, which we compute through a convection-diffusion equation:
Figure 7. Thermal fluid in an object. The temperature appears as a color contour; blue and red correspond to low and high temperatures. The user can determine the shape of the teapot while viewing the heat transport inside.
Table 2. Performance of our application examples tested with a 2.5-GHz CPU and 2.0 Gbytes of RAM. Application
Linear
Stationary
FPS (reuse)
FPS (no reuse)
Vibration
Yes
No
1,962
105.0
41.8
Cantilever
No
Yes
990
37.6
7.1
Fluid
No
No
1,971
30.0
18.2
Thermal fluid
No
No
1,938
21.3
14.3
r(Dv/Dt) = -∇p + m∇2v + rgb(T - T0) ∇⋅v=0 DT/Dt = a∇2T, where T is the temperature, T0 is the reference temperature, a is the thermal diffusivity, b is the volumetric thermal expansion ratio, and g is the acceleration of gravity. We again use the implicit method for time integration.
Application Performance Table 2 summarizes these applications’ performance. The frame rate depends on user manipulation (for example, the rate decreases when the user quickly makes a large shape modification). So, we averaged the frame rates for moderate-speed user manipulations similar to the manipulation in the “Algorithm Performance” section. The table shows that with data reuse, our applications ran at a high frame rate, sufficient for interactive modeling. We observed that our reuse scheme was especially effective for linear or stationary problems.
Informal User Studies To verify responsive FEM’s utility, we performed two informal user studies: one with a professional artist and one with nonprofessional students.
Metallophone Design To show that our current implementation is already useful as a practical tool for designing real-world objects, we used it for the design and creation of an artistically shaped metallophone. Metallophones are usually rectangular because that shape is suit
No. of elements
able for analytically predicting the instrument’s tone. However, researchers have used computer simulation to simulate the sound from arbitrarily shaped objects. For example, Jeffrey Chadwick and his colleagues proposed a computationally efficient framework for synthesizing realistic sound from nonlinear, thin-shell vibrations.5 However, we believe that designing a metallophone with a desired artistic shape and tone is possible only with responsive FEM, because this kind of highly constrained modeling naturally demands the tight integration of design and analysis. We model the metallophone as a 3D plate with a predefined thickness extruded from the 2D shape that the user designs. We compute its tone, or frequency, through eigenvalue analysis. We assume the lowest nonzero eigenfrequency is the output tone because higher eigenfrequencies usually attenuate quickly. Along with computing the eigenfrequency, we compute the eigenmode (which tells how the metallophone oscillates), which is helpful in determining which positions on the metallophone should be fixed. For details on these computations, see the sidebar “Computing the Metallophone’s Tone and Mode.” Simulation parameters must be calibrated for each specific material. Figure 8 shows our metallophone design software. Users edit the shape in 2D in the left window while checking the oscillation’s eigenmode in 3D in the right window. During design, the software updates the tone in real time, with visual feedback in the status bar and aural feedback using beeps. We asked a professional artist to design a metallophone using our software. We put no time limit IEEE Computer Graphics and Applications
49
Feature Article
Computing the Metallophone’s Tone and Mode
A
ssuming that the metallophone floats in nongravity space without any external forces and fixed boundaries, we use finite element method (FEM) discretization to obtain this equation: Mu + Ku = 0,(A)
where u is the nodal displacement vector, M is the lumped mass matrix, and K is the positive, semidefinite stiffness matrix. We split u into the product of the spatially varying amplitude f and the harmonic oscillation at the angle rate w, as u(x, t) = f(x)eiwt, where t indicates time and x indicates node positions. We then substitute this product into Equation A, which yields this generalized eigenvalue problem: Kφ = λMφ ,(B) where the eigenfrequency is f = ω (2π ) = λ (2π ) . We calculate the smallest nonzero eigenvalue l and its corresponding eigenvector f to obtain the metallophone’s tone and mode. We perform Cholesky factorization on the lumped mass matrix as M = LLT and multiply L –1 and L –T to the left and right sides of Equation A. Here, L stands for a lower triangular matrix. We use T as it’s typically used—indicating the transpose operation of the matrix. Thus we obtain this eigenvalue problem: Ay = ly,
Figure 8. Metallophone design software. Users edit the shape in the left window; the right window shows the analyzed eigenmode in 3D. The software updates the tone in real time with both aural and visual feedback.
where A = L –1KL –T and y = LTf. We solve this using inverse iteration. We modify the original iteration procedure by adding a step to remove the component of all zero eigenvectors of A from the current solution. That way, we obtain the smallest nonzero eigenvalue and its corresponding eigenvector. Thanks to our problem setting, we already know the zero eigenvectors of K as f0i (i = 1 … 6): the translations along the three coordinate axes and the rotations around them. We apply the modified Gram-Schmidt orthogonalization to LT f0i (i = 1 … 6) to obtain the orthonormal basis i vectors y0 (i = 1 … 6) that span the kernel of A. We also define a projection P that maps a vector v to the compliment space of the kernel of A as P (v ) = v − y0i (y0i ⋅ v ). In each iteration step, we apply this projection to the solution vector and normalize it. We add a small positive number e to the diagonals of A to improve the condition number. Once we’ve computed the shifted eigenvalue l1′ and its corresponding eigenvector y1 of A, we finally obtain the nonzero smallest eigenvalue λ1 = λ1′ − ε and the eigenvector f1= L –Ty1 in Equation B. Reusing the previous time step’s solution significantly improves the inverse iteration’s convergence. The simulation accuracy is sensitive to the mesh density because we use linear tetrahedral elements to represent the bend of a 3D thin plate. We can alleviate this problem, called shear locking, by preventing the excessive distortion of tetrahedral elements. We keep the longest edge shorter than twice the shortest edge.
∑
2D view (shape design)
3D view (eigenmode)
Scale and frequency
50
September/October 2011
C
D
E
F
G
A
B
(a)
(b) Figure 9. Designing a metallaphone. (a) The 2D shapes the artist designed. (b) The analyzed eigenmode of the designed pieces in 3D. For visualization purposes, the displayed eigenmode is considerably exaggerated. Using our system, the artist successfully designed complicated shapes that have specific tones.
on the design and provided instruction on the software. Figure 9 shows designed shapes corresponding to musical notes ranging from C (523 Hz) to B (987 Hz). We used a wire-electrical-discharge machine to cut the shapes from 4-mm aluminum plate; we fixed the shapes onto a wooden board according to the analyzed eigenmodes (see the right side of Figure 1). The top three rows of Table 3 show that for most pieces, the target, simulated, and actual frequencies were quite close. To further improve the metallophone’s quality, we manually adjusted the actual pieces’ tones by trimming their edges (we didn’t have to manually adjust the F piece’s tone). We used our software during the manual adjustment because it can predict tone changes caused by trimming. The bottom row of Table 3 shows the adjustments. We found that the actual metallophone produced sounds of sufficient quality for a hobbyist. After the design process, we interviewed the artist to obtain subjective comments. The artist reported that the design took roughly five hours, with most of the time devoted to the C and D pieces. Because these lower tones required larger areas, the analysis response time was much slower than for the other pieces. The artist had difficulty balancing the overall shape among the pieces while keeping their tones true; this was a huge design constraint. Additionally, sometimes a small shape modification resulted in a large change in tone, which demanded a highly responsive analysis. The artist commented that designing such an artistically shaped metallophone without using responsive FEM would be almost impossible. This prototype system is still a proof of concept, and many problems remain to be solved in designing a more realistic metallophone. For example, whereas we compute only fundamental frequencies, obtaining nicely tuned harmonics from secondary modes is necessary to produce a comfortable sound. Addressing these issues remains the subject of future research.
Table 3. The metallophone’s frequencies for each note, showing the accuracy of our analysis. Note Type of frequency Targeted
C
D
E
F
G
A
B
523
587
659
698
783
880
987
Simulated
525
588
661
699
786
880
989
Measured
506
604
621
698
787
860
993
Adjusted
523
587
659
698
783
881
987
Bridge Design The other user study concerned designing a bridge. The aim was to show that responsive FEM provides better support than traditional, nonresponsive FEM for nonprofessionals designing objects with physically desirable properties. The task. We asked the users to design the 2D shape of a bridge to span a certain gap and support a certain weight at its center, using as small an area as possible. That is, the goal was to design a strong bridge using the least amount of material. Figure 10 shows the overview of our bridge design system. The system provides a set of tools—such as curve editing using virtual pushpins,6 curve smoothing, and hole creation. It always displays the bridge’s area. The software tests the bridge’s strength through FEM analysis using a physical model based on equivalent stress. The system displays the amount of stress applied to each region as color contours (blue and red correspond to low and high stress) and judges whether the bridge passes the test. FEM analysis is available in two modes: responsive and nonresponsive. In responsive FEM, the analysis always occurs during user interaction (for example, mouse dragging), and the result updates in real time. In nonresponsive FEM, the analysis occurs only after the user completes the design and presses a toolbar button. The analysis result immediately disappears when the user changes IEEE Computer Graphics and Applications
51
Feature Article
Figure 10. Testing a bridge design’s structural soundness using FEM analysis. Blue and red contours correspond to low and high stress. Users were asked to design a bridge that wouldn’t break down under a weight located in its center, using the least amount of material. Responsive FEM
Nonresponsive FEM
Area
30 20 10
1
2 3 Participant
Avg.
4
5 6 Participant
Avg.
Subject
Figure 11. The bridge design study’s results. Subjects who used responsive FEM generally achieved better designs, with passing bridges encompassing smaller areas than those designed with nonresponsive FEM.
the design. The nonresponsive mode simulates the typical use of FEM systems, in which design and analysis are completely separate. In such systems, switching between these two processes can involve much tedious work, such as exporting and importing files and setting FEM parameters. Note that the nonresponsive FEM we used is already much more efficient than current commercial FEM tools. The experimental setup. Six university students majoring in art and design participated; none was familiar with FEM techniques or material mechanics. We split them into two groups: group A started with responsive FEM; group B started with nonresponsive FEM. We performed separate experiments for each group. Each group underwent a 15-minute lesson on using the software and then a 30-minute main design session. Afterward, we asked the participants in each group to try the other FEM mode. We then collected their subjective feedback. 52
September/October 2011
Results. Figure 11 shows the smallest-area bridge that passed the strength test for each participant. Participants who used responsive FEM generally achieved better results. The most common feedback was that responsive FEM is useful when the user wants to see how a small adjustment affects the analysis. The participants also mentioned that the analysis results displayed during responsive FEM were sometimes too conspicuous, making large design changes difficult because users felt encouraged to optimize results locally with small edits. Some participants even claimed that shape design without responsive FEM might be more appropriate for initial design exploration. This simple study obviously doesn’t “prove” anything, but we believe it suggests responsive FEM’s potential as a design aid for nonprofessionals. We plan to perform more formal user studies.
Limitations We must extend our techniques to 3D to make them useful for many practical real-world problems. Solving linear systems in 3D is straightforward; our main challenge will be working the continuous mesh update in 3D. As François Labelle and Jonathan Shewchuk noted, methods for improving tetrahedral mesh quality by continuously moving nodes and changing connectivity have yet to guarantee sufficiently accurate simulations.7 Another limitation is that currently we can use only first-order elements. Higher-order elements are much more desirable for such problems as bending thin-walled structures. However, using such elements in our approach is problematic because reconnecting the mesh will most likely change relationships between nodes and prevent matrix reuse. Additionally, we’ve yet to try nonsimplex elements (for example, quadrilaterals in 2D and hexahedrons in 3D). Such elements might—depending on the problem—be more appropriate than simplex elements (for example, triangles in 2D and tetrahedrons in 3D). However, reconnecting nonsimplex meshes without adding and deleting points is difficult. Our approach isn’t applicable to history-dependent problems (such as plasticity processing). In such problems, the solutions must be computed sequentially from the initial state, and our data reuse scheme is inappropriate for that purpose.
O
ur current applications already show potential utility in areas as diverse as aerodynamics, air conditioning layout, and airfoil design. Additional research could make our method even
more globally beneficial. To start with, we intend to test reuse of data other than the preconditioner matrix; these tests might include node reordering, which would improve analysis responsiveness. We could also apply our research to create a system that uses simulation results to more actively guide design. For example, the system could help users design shapes that satisfy certain constraints (for example, stress limits in certain areas) by providing suggestions and instructions when the user makes a change that doesn’t satisfy the constraints. We could also design a system in which the user interactively controls simulation accuracy. We use fixed criteria for the trade-off between speed and accuracy, but users might want more explicit control over this trade-off (for example, during design refinement). We also assume homogeneous mesh density, but it would be useful for users to locally control simulation accuracy by manipulating the local mesh density. In this way, they could examine analysis results more closely in specific regions of interest; existing automatic error estimation methods make this difficult.
References 1. I.E. Sutherland, “Sketch Pad: A Man-Machine Graphical Communication System,” Proc. SHARE Design Automation Workshop, ACM Press, 1964, pp. 6.329–6.346. 2. J. Mezger et al., “Interactive Physically Based Shape Editing,” Proc. 2008 ACM Symp. Solid and Physical Modeling, ACM Press, 2008, pp. 79–89. 3. N. Chentanez et al., “Interactive Simulation of Surgical Needle Insertion and Steering,” ACM Trans. Graphics, vol. 28, no. 3, 2009, article 88. 4. T.E. Tezduyar et al., “Incompressible Flow Compu tations with Stabilized Bilinear and Linear EqualOrder-Interpolation Velocity-Pressure Elements,” Computer Methods in Applied Mechanics and Eng., vol. 95, no. 2, 1992, pp. 221–242. 5. J.N. Chadwick, S.S. An, and D.L. James, “Harmonic Shells: A Practical Nonlinear Sound Model for NearRigid Thin Shells,” ACM Trans. Graphics, vol. 28, no. 5, 2009, article 119. 6. T. Igarashi, T. Moscovich, and J.F. Hughes, “AsRigid-as-Possible Shape Manipulation,” ACM Trans. Graphics, vol. 24, no. 3, 2005, pp. 1134–1141. 7. F. Labelle and J.R. Shewchuk, “Isosurface Stuffing: Fast Tetrahedral Meshes with Good Dihedral Angles,” ACM Trans. Graphics, vol. 26, no. 3, 2007, article 57.
computer simulation and interactive design. Umetani has an MS in frontier science from the University of Tokyo. Contact him at [email protected]. Kenshi Takayama is a PhD student in computer science at the University of Tokyo. His research interests include computer graphics, volumetric modeling, user interfaces, texture synthesis, geometric modeling, and image editing. Takayama has an MS in computer science from the University of Tokyo. Contact him at [email protected]. Jun Mitani is an associate professor of computer science at the University of Tsukuba. His research interests include computer graphics, geometric modeling, and origami. He’s a member of the Japan Society of Graphics Science and the Society for Art and Science of Japan. Mitani has a PhD in precision engineering from the University of Tokyo. Contact him at [email protected]. Takeo Igarashi is an associate professor of computer science at the University of Tokyo. He’s also the director of the Japan Science and Technology Agency ERATO (Exploratory Research for Advanced Technology) Igarashi Design Interface Project. His research interests include user interfaces and interactive computer graphics. Igarashi has a PhD in information engineering from the University of Tokyo. Contact him at [email protected].
Call rticles for A IEEE seek s re
a
uting p m o sive C he latest peerPerva on t
le, u ccessib
d de viewe
u
ou s biquit
s e f ul
velop
pap e
m e nt
co m p
s
uting
rs
vasiv in per
. To
e, mo
clud pics in
ra s t r re inf
u c tur
bile, a
e har
nd
d war
e , rea
e
l-wor
ld
wa u te r , s of t co m p ology man u h , ding n rac tio s , inclu d inte ration n e a id g s s co n sensin . s y s te m rivac y , an d an d p n io t c urit y, c e s , intera y labilit t , sca ymen deplo es: or.htm de lin /auth r gui vasive r o e h /p t Au rg /mc uter.o .comp www
te chn
er Furth ils: det a p e r va co m p
sive@
uter.o
rg
. www ter. u c om p ive er vas or g / p
Nobuyuki Umetani is a PhD student in computer science at the University of Tokyo. His research interests include
IEEE Computer Graphics and Applications
53
