A solid volumetric deformable muscle model for computer animation using weighted residual method

Comput. Methods Appl. Mech. Engrg. 190 (2000) 853±863

www.elsevier.com/locate/cma

A solid volumetric deformable muscle model for computer animation using weighted residual method L.H. You, J.J. Zhang *, P. Comninos National Centre for Computer Animation, Bournemouth University, Talbot Campus, Weymouth House, Fern Barrow Poole, Dorset BH12 5BB, UK Received 15 July 1999

Abstract One of the most important yet expensive tasks in realistic character animation is to represent the body deformation, i.e. the deformation of the muscles. Because of the complexity of the geometry and deformable behaviours of such muscle tissues, threedimensional (3D) ®nite element analysis (FEA) is usually employed. FEA however, requires enormous computing power, and often is impractical for computer animation where computational cost needs to be balanced with modelling accuracy. In this sense, it di€ers from engineering analysis where accuracy is often paramount for many applications, such as in aviation industry. Proposed in this paper is a 3D solid volumetric deformable muscle model which is almost as accurate as the FEA, but requires only a fraction of the computing cost. In this model, the muscle tissues are assumed to be elastic bodies with isotropic behaviour. The governing equations are derived and a weighted residual method is introduced to solve these equations. Suitable trial functions are proposed which satisfy the deformation compatibility and equilibrium equations. Using the least squares technique, the residual values on the boundaries are formulated and minimised to determine the unknown constants in the trial functions. In order to validate our model, the numerical results produced by our model are compared with the ®nite element computation and it is found that the two results closely agree with each other. Ó 2000 Elsevier Science S.A. All rights reserved. Keywords: Volumetric deformation; Muscle model; Weighted residual method; Computer simulation

1. Introduction Computer animation of deformable characters, such as humans and animals, has extensive applications in ®lm special e€ects, digital entertainment, commercials and surgical simulation. Primarily, there are two classes of muscle models: geometric models and physically based models. 1.1. Geometric models The geometric models do not consider the physical properties of the muscles. Therefore, real-time deformation in computer animation can be achieved more easily and cheaply than the physically based models. The major disadvantage however, is the diculty of achieving convincing realism. Most early models are geometric ones, as the computation cost of a physically based model was not practicable with the hardware performance.

*

Corresponding author. Tel.: +44-1202-595055; fax: +44-1202-595530. E-mail address: [email protected] (J.J. Zhang).

0045-7825/00/$ - see front matter Ó 2000 Elsevier Science S.A. All rights reserved. PII: S 0 0 4 5 - 7 8 2 5 ( 9 9 ) 0 0 4 4 8 - X

854

L.H. You et al. / Comput. Methods Appl. Mech. Engrg. 190 (2000) 853±863

There are a number of geometric models in the existing literature. Komatsu [1] parameterised the control points of spline patches to simulate a contracting biceps when an elbow is bent. Dow and Semwal [2] developed a method for modelling the muscles and bones in the human body. In their method, each muscle and bone is modelled by a generalised cylinder. For each frame in an animation, three steps are taken: (1) the shape of a separate muscle is determined, (2) intersections among the muscles are detected and corrected, and (3) the composite shape is wrapped with a skin surface. Shen et al. [3] proposed a method to create and animate a human body based on the human anatomy. Their model uses ellipsoids to model the deformed muscles and automatically determines the skin using crosssectional contours. 1.2. Physically based models To improve the realism of computer animation of human ®gures, an increasing amount of research has been carried out in physically based models. Roughly speaking, these models can be categorised as surfacevolumetric or solid-volumetric models. Surface-volumetric models represent human muscle tissue as a set Hookean springs attached to the skin surface. Most physically based models reported in the literature fall in this category. Terzopoulos and Waters [4] proposed a deformable model that approximates the mechanical properties of real human facial tissue. This deformable model is a lattice of nonlinear springs connecting point masses and resolved by numerically integrating a system of di€erential equations of motion. Celniker and Gossard [5] developed an alternative deformable model that combines a parametrically described geometry with an energy minimisation algorithm. The ®nite element method was used to obtain the overall energy of the deformable curve or surface. This approach automatically adjusts the deformable shape by minimising its energy considering the geometric constraints and loads. Subsequently, Koch et al. [6] adopted this method to simulate facial surgery. In their approach, a ®nite element surface model was used to represent the facial surface. The deformation of di€erent types of facial surface was computed by minimising the overall energy, subject to external and internal loading forces. Following the modelling approach proposed by Terzopoulos and Waters [4], Lee et al. [7] further developed a system which generates facial expressions by constructing synthetic muscles embedded in an anatomically motivated skin model composed of three spring-mass layers. Wu and Thalmann [8,9] developed a plastic-visco-elastic skin surface model that can animate facial wrinkles and skin ageing. In their model, the skin can slide over a supporting layer which constrains the skin surface by a set of spring forces that represent the connective fat tissue between the skin and the muscle. Rather than using a spring system, solid-volumetric models represent an object as a solid together with its physical properties. Comparatively speaking, one obvious advantage of the surface-volumetric models is their computational eciency, as the resolution of a spring-mass system is less computationally intensive than the ®nite element calculation. This property is crucially important for computer animation. Human muscular tissues however are very complicated. Their structures and deformation properties are very dif®cult to be represented accurately by a spring-mass system. In addition, surface-volumetric models have a detectable non-solid behaviour, which can cause problems in certain applications (for instance, in the simulation of surgical procedures). Clearly, a surgeon cannot cut a virtual organ that is constructed using a surface-volumetric model, as there is nothing inside the surface. Thus, solid-volumetric models, despite being more computationally expensive, become increasingly popular. Sagar et al. [10] created an anatomically detailed 3D computer graphic model of human eyes and the surrounding tissues and implemented a virtual environment for their surgical simulation. This simulation is implemented using a ®nite element model based on the orthotropic elasticity theory of large deformation. Chen and Zeltzer [11] developed a novel computational method of skeletal muscle. The geometry and material properties were captured using a 3D elastic ®nite element method. A biomechanical model of muscle was used to apply non-linear forces to the ®nite element mesh nodes. In real-time simulation, solid-volumetric models encounter a major problem, namely the computing speed. This is due to the need of a large number of 3D ®nite elements for these models that inevitably make the computation very expensive and time consuming. So far, the ®nite element analysis (FEA) models used in computer animation are still con®ned to simple elastic behaviours. Nielsen and Cotin [12] applied a 3D

L.H. You et al. / Comput. Methods Appl. Mech. Engrg. 190 (2000) 853±863

855

elastic solid ®nite element model to surgery simulation and studied the problem of real-time performance. The simulation speed of solid-volumetric models remains an unsolved problem. In this paper, we propose a 3D solid-volumetric model for the elastic deformation of human muscle. Instead of employing a FEA model, which unavoidably necessitates high computational cost, our model adopts a cheaper numerical approach aiming to achieve fast performance without losing noticeable accuracy.

2. Governing equations In a 3D problem, if we assume the displacement components along the x, y and z directions to be u, v and w, respectively, then the equilibrium equation with respect to displacements can be derived from the strain±displacement relation, the stress±strain relation and the equilibrium equations. Ignoring volumetric forces, the equilibrium equation with respect to displacements can be written as: G oh r2 u ‡ ˆ0 k‡G ox G oh r2 v ‡ ˆ0 k‡G oy G oh r2 w ‡ ˆ0 k‡G oz

…1†

where G is the shear modulus, k the LameÕs constant, h the volume expansion coef®cient and r2 is the Laplacian operator. They are de®ned as: E 2…1 ‡ m† mE kˆ …1 ‡ m†…1 ÿ 2m† ou ov ow ‡ ‡ hˆ ox oy oz o2 o2 o2 r2 ˆ 2 ‡ 2 ‡ 2 ox oy oz

Gˆ

where E is the YoungÕs modulus and m is the PoissonÕs ratio. The deformation compatibility equations are given as: o2 ey o2 ez o2 cyz ‡ 2 ˆ oz2 oy oy oz o2 ez o2 ex o2 czx ‡ 2 ˆ ox2 oz oz ox 2 2 o ex o ey o2 cxy ‡ 2 ˆ oy 2 ox ox oy   ocyz oczx ocxy o o2 e x ‡ ‡ ˆ2 ÿ ox oy oz oyoz ox   ocxy ocyz o oc o2 e y ‡ ˆ2 ÿ zx ‡ oy oz ox ozox oy   ocxy ocyz oczx o o2 e z ‡ ‡ ˆ2 ÿ oz ox oy oxoy oz where ex ; ey ; ez are the normal strains and cyz ; czx ; cxy the shear strains.

…2†

856

L.H. You et al. / Comput. Methods Appl. Mech. Engrg. 190 (2000) 853±863

The resolution of an elastic problem requires to ®nd suitable displacement functions, which satisfy the equilibrium equation (1), the deformation compatibility equation (2), and the following displacement boundary conditions: u ˆ us ;

v ˆ vs ;

w ˆ ws

…3†

as well as the external force boundary conditions: lrx ‡ msxy ‡ nsxz ˆ Xs ;

lsxy ‡ mry ‡ nsyz ˆ Ys ;

lsxz ‡ msyz ‡ nrz ˆ Zs

…4†

where rx ; ry ; rz are the normal stresses and syz , szx ; sxy the shear stresses; us ; vs ; ws are the speci®ed displacement values on the displacement boundary S; l, m, n are the direction cosines of the boundary S; and Xs ; Ys ; Zs are the speci®ed stress values on the external force boundary S. The equilibrium equation (1) and the deformation compatibility equation (2) can be satis®ed if the displacement functions take the following form: uˆH

oX ‡ r2 / 1 ; ox

vˆH

1 ; 2…1 ÿ m†

Xˆ

oX ‡ r2 /2 ; oy

wˆH

oX ‡ r2 /3 oz

…5†

where H ˆÿ

o/1 o/2 o/3 ‡ ‡ ox oy oz

and /1 ; /2 and /3 are known as trial functions, which satisfy the following biharmonic equations: r4 /1 ˆ 0;

r4 /2 ˆ 0;

r4 / 3 ˆ 0

…6†

By introducing the above displacement functions, the resolution of a 3D elastic problem is now transformed into ®nding the suitable trial functions /1 ; /2 and /3 which satisfy the biharmonic equations (6), the external force boundary conditions (4) and the displacement boundary conditions (3).

3. Weighted residual solution of governing equations For many practical problems, it is very dicult to ®nd an exact analytical solution satisfying the governing equations and their boundary conditions. This is not to say, however, that such problems are unsolvable. As a matter of fact, apart from the FEA method, various approximately analytical methods and semi-analytical methods are promising alternatives [13,14,17]. These alternative approaches are in general much less expensive than FEA. The weighted residual method is such a semi-analytical method. The key to its success is the choice of the trial functions. Similar to our previous work [15,16], in this paper, the following trial functions are used, which automatically satisfy the biharmonic equations (6): /1 ˆ

k X

sin nx cos ny…ci1 egz ‡ ci2 eÿgz ‡ ci3 zegz ‡ ci4 zeÿgz †

iˆ1

‡

k X

sin ny cos nz…ci5 egx ‡ ci6 eÿgx ‡ ci7 xegx ‡ ci8 xeÿgx †

iˆ1

‡

k X iˆ1

sin nz cos nx…ci9 egy ‡ ci10 eÿgy ‡ ci11 yegy ‡ ci12 yeÿgy †

L.H. You et al. / Comput. Methods Appl. Mech. Engrg. 190 (2000) 853±863

/2 ˆ

k X

857

sin nx cos ny…ci13 egz ‡ ci14 eÿgz ‡ ci15 zegz ‡ ci16 zeÿgz †

iˆ1

‡

k X

sin ny cos nz…ci17 egx ‡ ci18 eÿgx ‡ ci19 xegx ‡ ci20 xeÿgx †

iˆ1

‡

k X

sin nz cos nx…ci21 egy ‡ ci22 eÿgy ‡ ci23 yegy ‡ ci24 yeÿgy †

iˆ1

/3 ˆ

k X

…7† gz

sin nx cos ny…ci25 e ‡ ci26 e

ÿgz

gz

‡ ci27 ze ‡ ci28 ze

ÿgz

†

iˆ1

‡

k X

sin ny cos nz…ci29 egx ‡ ci30 eÿgx ‡ ci31 xegx ‡ ci32 xeÿgx †

iˆ1

‡

k X

sin nz cosnx…ci33 egy ‡ ci34 eÿgy ‡ ci35 yegy ‡ ci36 yeÿgy †

iˆ1

p where n ˆ ip=L; g ˆ 2n, e is the natural logarithm base constant, L the maximum size of the elastic object in question, /1 ; /2 ; /3 the trial functions and cij …i ˆ 1; 2; . . . ; k; j ˆ 1; 2; . . . ; 36† are the unknown constants. Since the biharmonic equations (6) are satis®ed, after substituting the trial functions into them, only the residual values of Eqs. (3) and (4) need to be minimised, which will be carried out using the discrete least squares method. Assuming that there are p collocation points on all displacement boundaries and q collocation points on all external force boundaries, we substitute the co-ordinate values of these points into Eq. (7). By making use of Eqs. (3) and (4), we obtain the following set of linear equations representing the residual values: Ruj …c; xj † ˆ H

 xj † oX…c; ‡ r2 /1 …c; xj † ÿ us ox

Rvj …c; xj † ˆ H

 xj † oX…c; ‡ r2 /2 …c; xj † ÿ vs oy

Rwj …c; xj † ˆ H

 xj † oX…c; ‡ r2 /3 …c; xj † ÿ ws oz

…j ˆ 1; 2; . . . ; p†

…8†

rx …c; xj † ‡ m sxy …c; xj † ‡ n sxz …c; xj † ÿ Xs RXj …c; xj † ˆ l RYj …c; xj † ˆ l sxy …c; xj † ‡ m ry …c; xj † ‡ n syz …c; xj † ÿ Ys RZj …c; xj † ˆ l sxz …c; xj † ‡ m syz …c; xj † ‡ n rz …c; xj † ÿ Zs

…j ˆ 1; 2; . . . ; q†

where c is the matrix of the unknown constants cij …i ˆ 1; 2; . . . ; k; j ˆ 1; 2; . . . ; 36†, and xj represents the jth collocation point, and  xj † ˆ o/1 …c; xj † ‡ o/2 …c; xj † ‡ o/3 …c; xj † X…c; ox oy oz Eqs. (8) can also be written in matrix form R ˆ AC ÿ B

…9†

858

L.H. You et al. / Comput. Methods Appl. Mech. Engrg. 190 (2000) 853±863

where R, C and B are 3…p ‡ q†  1 matrix consisting of the residual values, unknown constants and constant terms, and A is a 3…p ‡ q†  36k matrix consisting of the coecients of the unknown constants cij …i ˆ 1; 2; . . . ; k; j ˆ 1; 2; . . . ; 36†. The square sum of the residual values in Eq. (9) gives I ˆ RT R

…10†

The minimum of Eq. (10) is determined by oI ˆ0 oC

…11†

which gives the following linear equations: AT AC ˆ AT B

…12†

The resolution of (12) gives the unknown constants and hence the elastic deformation of the problem. 4. Numerical examples The above proposed model allows us to compute the deformation of an elastic object using a semianalytical form. Compared with the numerical approaches, notably 3D FEA, the obvious advantage is its computational eciency, which is imperative in the animation of deformable characters. Apart from computational eciency, the other important factor is the computational accuracy. In order to verify the accuracy of our model, we decided to make a quantitative comparison between our method and the FEA computation. Instead of using complicated models, we chose a simple geometric object, a block, as its deformation can be accurately obtained and hence the results are recognisable. The mechanical properties of the block were chosen the same as those of the muscles tissues used by Chen et al. [11] in their research. The YoungÕs modulus and PoissonÕs ratio were taken to be E ˆ 100 Pa and m ˆ 0:49, respectively. The FEA computation was performed using the commercial FEA package ADINA version 82. Case 1. A group of uniformly distributed end tensile force is applied to one end of the block. The other end …x ˆ 0† is fully ®xed as shown in Fig. 1(a). The dimension ratio of the three sides is 20:1:1 and the applied tensile force is 2 Pa. On each end face, 81 collocation points are uniformly distributed and the constant k is taken to be 7 in the trial functions (7). The resultant relative displacements in the x-direction at

Fig. 1. Point distributions for displacement list at the speci®ed section.

L.H. You et al. / Comput. Methods Appl. Mech. Engrg. 190 (2000) 853±863

859

points 1; . . . ; 25 (see Fig. 1(b)) of section x ˆ 20 are given in Table 1. In this table, NOP denotes the number of points, FEM represents the solution arrived at by the ®nite element method and WRM represents the solution achieved by our method. The initial and deformed ®gures are shown in Fig. 2(a), (b) and (c), respectively, where Fig. 2(b) was produced using the FEA method and Fig. 2(c) was produced by the proposed method. Case 2. The same block is now under an end tensile displacement, rather than force. The dimension ratio of the three sides is chosen to be 5:1:1 and the applied relative displacement (i.e. the ratio of the speci®ed end displacement to the unit length) is 0.5. The same boundary conditions and collocation points as those in the ®rst case were applied. Constant k was set to 5 in the trial functions (7). The resultant relative

Table 1 The displacements in x-direction at section x ˆ 20 NOP

FEM

WRM

Error (%)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25

3.43757 3.43757 3.43757 3.43757 3.43757 3.43757 3.43757 3.43757 3.43757 3.43757 3.43757 3.43757 3.43757 3.43757 3.43757 3.43757 3.43757 3.43757 3.43757 3.43757 3.43757 3.43757 3.43757 3.43757 3.43757

3.38108 3.38107 3.38097 3.38078 3.38052 3.40540 3.40545 3.40540 3.40524 3.40499 3.43025 3.43036 3.43035 3.43022 3.42999 3.45550 3.45566 3.45568 3.45558 3.45537 3.48105 3.48123 3.48127 3.48119 3.48102

1.64 1.64 1.65 1.65 1.66 0.94 0.93 0.94 0.94 0.95 0.21 0.21 0.21 0.21 0.22 0.52 0.53 0.53 0.52 0.52 1.26 1.27 1.27 1.27 1.26

Fig. 2. Deformation under end tensile force.

860

L.H. You et al. / Comput. Methods Appl. Mech. Engrg. 190 (2000) 853±863

displacements of points 1; . . . ; 25 at section x ˆ 4:375 in the x-direction are given in Table 2. The initial and deformed geometry are shown in Fig. 3(a), (b) and (c), respectively. Case 3. In the third case study we apply an end shear deformation to the same object with the same boundary conditions. The relative end displacement is chosen to be 0.6 and constant k to be 9 used in the trial function. The computed relative displacements in the x-direction at points 1; . . . ; 25 for section x ˆ 4:375 are listed in Table 3. The initial and deformed ®gures are shown in Fig. 4(a)±(c).

Table 2 The displacements in x-direction at section x ˆ 4:375 NOP

FEM

WRM

Error (%)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25

0.42718 0.42686 0.42666 0.42686 0.42718 0.42686 0.42653 0.42633 0.42653 0.42686 0.42666 0.42633 0.42614 0.42633 0.42666 0.42686 0.42653 0.42633 0.42653 0.42686 0.42718 0.42686 0.42666 0.42686 0.42718

0.37919 0.37905 0.38237 0.38914 0.40026 0.39106 0.38860 0.38925 0.39317 0.40012 0.40082 0.39659 0.39455 0.39530 0.39790 0.40795 0.40240 0.39819 0.39613 0.39478 0.41323 0.40664 0.40094 0.39669 0.39185

11.23 11.20 10.38 8.84 6.30 8.39 8.89 8.70 7.82 6.2 6.06 6.98 7.41 7.28 6.74 4.43 5.66 6.60 7.1 7.52 3.27 4.74 6.03 7.07 8.27

Fig. 3. Deformation under end tensile displacement.

L.H. You et al. / Comput. Methods Appl. Mech. Engrg. 190 (2000) 853±863

861

Table 3 The displacements in x-direction at section x ˆ 4:375 NOP

FEM

WRM

Error (%)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25

0.47576 0.47677 0.47710 0.47677 0.47576 0.47490 0.47601 0.47636 0.47601 0.47490 0.47459 0.47573 0.47609 0.47573 0.47459 0.47490 0.47601 0.47636 0.47601 0.47490 0.47576 0.47677 0.47710 0.47677 0.47576

0.48055 0.48956 0.49650 0.50354 0.51210 0.47955 0.48854 0.49497 0.50103 0.50822 0.48200 0.48890 0.49411 0.49934 0.50575 0.48599 0.48922 0.49290 0.49765 0.50348 0.48830 0.48727 0.48993 0.49576 0.50234

1.01 2.68 4.07 5.61 7.64 0.98 2.63 3.91 5.26 7.02 1.56 2.77 3.78 4.96 6.58 2.34 2.78 3.47 4.55 6.02 2.63 2.20 2.69 3.98 5.59

Fig. 4. Deformation under end shear displacement.

4.1. Observations From the case studies, it is clear that the proposed model o€ers satisfying accuracy. The maximum errors compared with the FEA method are only around 10% and on the average, the error is less than 5%. This is accurate enough for computer animation. Since the weighted residual method is a semi-analytical method, substantially fewer unknowns than those of the ®nite element method are introduced. Hence less data storage and CPU resource are required.

862

L.H. You et al. / Comput. Methods Appl. Mech. Engrg. 190 (2000) 853±863

Table 4 The number of the unknowns for both methods Case number

1

2

3

FEM WRM

2187 252

2187 180

2187 324

The unknown constants of the ®nite element calculation and the weighted residual analysis for the above three cases are listed in Table 4. It is clear that the proposed model is much cheaper. 5. Conclusion and discussion Two factors are crucial to the computer animation of deformable characters: visual realism and computing eciency. It is therefore not the same as those engineering applications where precision cannot be compromised. Based on this understanding, we have developed a solid volumetric deformable muscle model, which primarily aims at its eciency, while also satisfying the accuracy. This model was derived from the theory of elasticity. To resolve this model and obtain the deformation formula of an object, the weighted residual method was introduced. This method works by choosing suitable trial functions, which automatically eliminate the residual values (errors) inside the 3D elastic object. The formula coecients can then be determined by minimising the residual values only on the boundaries. It is worth pointing out that since the coecients of the formulae relate merely to the boundary conditions, it is especially advantageous to animation applications. As far as the animated characters are concerned, usually just a few boundary conditions may change during the course of animation and most of them will stay unchanged. Thus only those residual values with respect to the a€ected boundaries need recalculation. As a consequence, only a small number of coecients need to be recalculated. Unlike the FEA method, the results from the weighted residual method are stored in the form of coecients taken in the analytical formulae. As it is no longer necessary to store the displacements of the objects, much less memory and computation are required. Therefore our method is more attractive for realtime modelling applications. To con®rm the accuracy of the proposed model, the FEA method was used to perform the calculation of some numerical examples. Close agreement was found between our method and the FEA.

References [1] K. Komatsu, Human skin model capable of natural shape variation, The Visual Computer 3 (1988) 265±271. [2] D.E. Dow, S.K. Semwal, A framework for modeling the human muscle and bone shapes, in: Proceedings of the Third International Conference on CAD and Computer Graphics, vol. 1, 1993, pp. 110±115. [3] J. Shen, N.M. Thalmann, D. Thalmann, Muscle based human body deformation, in: Proceedings of the Third International Conference on CAD and Computer Graphics, vol. 1, 1993, pp. 95±100. [4] D. Terzopoulos, K. Waters, Physically based facial modelling analysis and animation, The Journal of Visualization and Computer Animation 1 (1990) 73±80. [5] G. Celniker, D. Gossard, Deformable curve and surface ®nite-elements for free-form shape design, Computer Graphics 25 (4) (1991) 257±266. [6] R.M. Koch, M.H. Gross, F.R. Carls, D.F. Von B uren, G. Fankhauser, Y.I.H. Parish, Simulating facial surgery using ®nite element models, in: Computer Graphics Proceedings, Annual Conference Series, 1996, pp. 421±428. [7] Y. Lee, D. Terzopoulos, K. Waters, Realistic modelling for facial animation, in: Computer Graphics Proceedings, Annual Conference Series, 1995, pp. 55±62. [8] Y. Wu, D. Thalmann, A plastic±visco-elastic model for wrinkles in facial animation and skin ageing, in: Proceedings of the Second Paci®c Conference on Computer Graphics and Applications, Paci®c Graphics Õ94, Fundamentals of Computer Graphics, 1994, pp. 201±213. [9] Y. Wu, N.M. Thalmann, D. Thalmann, A dynamic wrinkle model in facial animation and skin ageing, The Journal of Visualization and Computer Animation 6 (1995) 195±205.

L.H. You et al. / Comput. Methods Appl. Mech. Engrg. 190 (2000) 853±863

863

[10] M.A. Sagar, D. Bullivant, G.D. Mallinson, P.J. Hunter, A virtual environment and model of the eye for surgical simulation, in: Computer Graphics Proceedings Annual Conference Series 1995, pp. 205±212. [11] D.T. Chen, D. Zeltzer, Pump it up: Computer animation of a biomechanically based model of muscle using the ®nite element method, Computer Graphics 26 (2) (1992) 89±98. [12] M.B. Nielsen, S. Cotin, Real-time volumetric deformable models for surgery simulation using ®nite elements and condensation, Eurographics Õ96 15 (3) (1996) 57±66. [13] L.H. You, J.J. Zhang, H.B. Wu, R.B. Sun, A simpli®ed calculating method of wall plates of dry gas holders under gas pressure, International Journal of Pressure Vessels and Piping 74 (1) (1997) 13±18. [14] L.H. You, S.Y. Long, E€ects of material properties of interfacial layer on stresses in ®brous composite subjected to thermal loading, Composite A 29 (1998) 1185±1192. [15] L.H. You, H.B. Wu, R.B. Sun, S.Y. Long, A series method calculating the body of dry gas holders under gas pressure, Journal of International Association for Shell and Spatial Structures 39 (1) (1998) 15±24. [16] L.H. You, J.J. Zhang, Elastic±plastic stresses in a rotating solid disk, International Journal of Mechanical Science 41 (3) (1999) 269±282. [17] L.H. You, S.Y. Long, L. Rohr, Elastic±plastic stress ®eld in a coated continuous ®brous composite subjected to thermomechanical loading, Journal of Applied Mechanics, Transactions of the ASME 66 (3) (1999) 750±757.