- Email: [email protected]
BEM higher-order surfaces using Java 3D textures
Advances in Engineering Software 34 (2003) 469–476 www.elsevier.com/locate/advengsoft
Generating contours on FEM/BEM higher-order surfaces using Java 3D textures G.P. Nikishkov* The University of Aizu, Aizu-Wakamatsu City, Fukushima 965-8580, Japan Received 10 December 2002; revised 24 February 2003; accepted 20 March 2003
Abstract An efficient technique to visualize primary and secondary results for combined finite element method/boundary element method models as contours is presented. The technique is based on dividing higher-order surfaces into triangles and on using texture interpolation to produce contour plots. Since results of high accuracy with significant gradients can be obtained using sparse meshes of boundary elements and finite elements, special attention is devoted to element face subdivision. Subdivision density is defined on the basis of both face edge curvature and ranges of result fields over element faces. Java 3D API is employed for code development. q 2003 Elsevier Science Ltd. All rights reserved. Keywords: Finite element method; Boundary element method; Visualization; Contours; Java 3D
1. Introduction The finite element method (FEM) and the boundary element method (BEM) are used for computational modeling for several decades. In the beginning, simple elements with linear approximation of sought functions were employed. Later, simple elements have been replaced by higher-order elements with quadratic (or even cubic) approximation due to their computational efficiency. At present, the most used elements in three-dimensional analysis are hexahedral finite element with 20-nodes and quadrilateral boundary element with 8-nodes. Both elements are built on quadratic interpolating functions. For some specific problems, it is advantageous to apply a combination of FEM and BEM. In Ref. [1], the combined FEM –BEM technique is used to analyze non-planar cracks in three-dimensional components. The FEM is used for stress analysis of the uncracked structural component and the symmetric Galerkin boundary element method (SGBEM) is employed to model the crack. According to the superposition principle, the solution is sought in an iteration procedure with alternating between the FEM and SGBEM model. Using quadratic elements for both FEM and BEM models and taking into account the ability of the BEM * Tel.: þ81-242-37-2649; fax: þ 81-242-37-2706. E-mail address: [email protected] (G.P. Nikishkov).
to model singular crack fields efficiently, it is possible to obtain displacement and stress fields with high precision even for sparse finite element and boundary element meshes. Visualization techniques for FEM and BEM models and analysis results have been developed during last decades. Commercial finite element codes such as ANSYS [2] usually have postprocessors with various options for result visualization. The Visualization Toolkit VTK [3] and the AVS system [4] provide many opportunities for field visualization. From this point of view, a researcher working with the FEM or BEM has some choice of tools for visualization of FEM/BEM results. On the other hand, available visualization tools cannot satisfy all researcher’s needs especially in cases when some non-standard methods used (such as SGBEM– FEM) or when the visualization should accompany calculations in real time. In these cases, specialized lightweight software may be useful to develop. This software can take into account specific features of the computational procedure and the results. The development of the customized visualization software can be attractive if the development cost is reasonable. Previously, when graphics libraries supported just simple graphical primitives, the development of visualization software was a complicated task. Recently developed Java 3D application programming interface (API) [5] makes the development of visualization software simpler. The Java 3D
0965-9978/03/$ - see front matter q 2003 Elsevier Science Ltd. All rights reserved. doi:10.1016/S0965-9978(03)00052-8
470
G.P. Nikishkov / Advances in Engineering Software 34 (2003) 469–476
API contains Java classes for a sophisticated threedimensional graphics rendering. The programmer works with high-level interface for creating and manipulating 3D geometric objects. Despite all its wide functionality, the API can be easily learned and used. The Java 3D employs linear triangular and quadrilateral filled polygons for surface representation. Because of this, the visualization of FEM/ BEM models consisting of simplest elements is almost straightforward. However, for higher-order elements, the transformation of element faces into triangular polygons should be done carefully taking into account both geometry and result fields. In this paper, a simple and at the same time efficient technique to visualize primary and secondary FEM/BEM result fields as color contours is presented. The technique is based on dividing higher-order surfaces into triangles and on using texture interpolation to produce contour plots. Since results of high accuracy with significant gradients can be obtained using sparse meshes of boundary or finite elements, special attention should be devoted to element face subdivision. We define subdivision density on the basis of both face edge curvature and ranges of result fields over element faces. After determining number of subdivisions along element face sides, control points are seeded inside the element face and then Delaunay triangulation is used to produce subdivision. Contours of results are not generated explicitly but obtained as a result of texture interpolation. The developed Java 3D code is illustrated with examples of visualization of FEM results for a valve and FEM/BEM results for a cube with an embedded circular crack.
2. Overview of the algorithm A computational procedure of the combined SGBEM– FEM method for fracture mechanics analysis is based on the superposition principle. According to this principle, a problem for a finite body with a crack is replaced by a superposition of two problems. The first problem for a finite body with a crack is solved by the FEM. The second problem for a crack in an infinite medium is modeled using SGBEM. For a correct superposition corresponding to the solution for a finite body with a crack, fictitious forces on the boundary of the finite element model should be found in order to compensate for the stresses caused by the presence of a crack in an infinite body. The alternating method [1] provides for an efficient solution, without assembling the joint SGBEM– FEM matrix. The SGBEM – FEM alternating method alternates between the finite element solution for an uncracked body and the boundary element solution for a crack in an infinite body. Usually, several iterations are sufficient for convergence. It should be noted that because of ability of the SGBEM to efficiently model singularities, displacements and stresses of high precision might be obtained using relatively sparse meshes of finite elements
Fig. 1. Texture interpolation inside a triangle.
and boundary elements. It the vicinity of the crack, stress fields are characterized by high gradients. The input data for the visualization consists of a set of nodes defined by spatial coordinates, a set of elements that is specified by nodal connectivities, and a set of result values. Primary results are obtained at nodes after solution of the global equation system. Secondary results, which are expressed through derivatives of the primary results, usually have the best precision at some points inside elements. In elasticity problems, the primary results are displacements, and the secondary results are strains and stresses. The visualization algorithm consists of the following main steps. 1. Obtain continuous field of results by extrapolation from integration points inside elements to element nodes with subsequent averaging (omit this step for primary results defined at nodal points). 2. Create the surface of the finite element model (omit this step for boundary element meshes). 3. Subdivide curved faces into flat triangles on the basis of face curvature and gradient of results. 4. Create contour pictures inside triangles by specifying texture coordinates of one-dimensional color pattern at triangle vertices. The last step of the visualization algorithm is performed with the use of texture interpolation available in Java 3D API. Instead of explicit contour drawing, we prepare onedimensional color texture and specify texture coordinates at triangle vertices as shown in Fig. 1.
3. Geometry and results data It is supposed that the finite element model consists of brick-type 20-node elements shown in Fig. 2a and the boundary element model is composed of 8-node quadrilateral
G.P. Nikishkov / Advances in Engineering Software 34 (2003) 469–476
471
node i; j; h; z are local coordinates ð21 # j; h; z # 1Þ; ji ; hi ; zi are nodal values of the local coordinates. Faces of 20-node finite elements have the same geometry as 8-node boundary elements (Fig. 2b). Quadratic two-dimensional shape functions are employed for coordinate and result interpolation inside 8-node quadrilateral element: Ni ¼ ð1 þ j0 Þð1 þ h0 Þðj0 þ h0 2 1Þ=4 Fig. 2. (a) Brick-type 20-node finite element and (b) 8-node boundary element.
for corner nodes Ni ¼ ð1 2 j2 Þð1 þ h0 Þ=2
for ji ¼ 0
2
elements (Fig. 2b). It can be easily noticed that the finite element faces are geometrically same as the boundary elements. Therefore, the visualization of results for the surface of the finite element model and for the boundary element model can be performed in a similar way. The finite element model is defined by two arrays: nodal coordinate array and connectivity array. Primary results (such as displacements in elasticity problems) are related to nodes. Secondary results (such as stresses depending on derivatives of primary results) have the most precise values at reduced integration Gaussian points 2 £ 2 £ 2 inside finite elements. For model and results visualization with Java 3D the following geometry information should be prepared: 1. element faces located at the surface of the finite element model and divided into triangles of appropriate size; 2. surface normal vectors for triangle vertices; 3. surface element edges divided into line segments of appropriate size; 4. texture coordinates for the one-dimensional color pattern interpolation (defined by results values at vertices of triangles).
ð3Þ
Ni ¼ ð1 2 h Þð1 þ j0 Þ=2
j0 ¼ jji ;
for hi ¼ 0
h0 ¼ hhi
A normal to the surface at any point with the specified local coordinates is calculated as a vector product of two tangent vectors along j and h directions: n~ ¼ ~e1 £ ~e2 ›x 1 ~e1 ¼ ›j ›x 1 ~e2 ¼ ›h
›x 2 ›j ›x 2 ›h
›x 3 ›j ›x 3 ›h
ð4Þ
Derivatives of the global coordinates xi in respect to the local coordinates j; h; z are determined with the use of the shape functions: X ›Nm m ›x i ¼ x ›j ›j i X ›Nm m ›x i ¼ x ›h ›h i
ð5Þ
3.1. Geometry relations for 20-node element
3.2. Surface of the finite element model
Quadratic shape functions are used to interpolate coordinates and results inside 20-node finite element: X X xk ¼ Ni xik ; f ¼ Ni fi ð1Þ
To visualize the finite element model or results at the surface of the model, it is necessary to select finite element faces that belong to the surface of the finite element model (or part of the finite element model). A topology of the finite element model is described by element connectivities. Connectivities for each finite element are global (model) node numbers listed in order of local (element) node numbers. The 20-node brick-type element has six element faces. Connectivity numbers for each face can be easily extracted from element connectivity array. The surface of the finite element model is created from outer element faces. Outer faces are mentioned in the model connectivity array only once while inner faces are mentioned exactly two times. A fast procedure to select outer element faces is to characterize each element face by two connectivity numbers instead of eight. To avoid
Ni ¼ ð1 þ j0 Þð1 þ h0 Þð1 þ z0 Þðj0 þ h0 þ z0 2 2Þ=8
ð2Þ
for vertex nodes Ni ¼ ð1 2 j2 Þð1 þ h0 Þð1 þ z0 Þ=4
for ji ¼ 0
Ni ¼ ð1 2 h2 Þð1 þ j0 Þð1 þ z0 Þ=4
for hi ¼ 0
2
Ni ¼ ð1 2 z Þð1 þ j0 Þð1 þ h0 Þ=4
j0 ¼ jji ;
h0 ¼ hhi ;
for ji ¼ 0
z0 ¼ zzi
where xk are coordinates of a point inside the element; xik is the kth coordinate of the ith element node; f is the result function value at a point; fi is the result function value at
472
G.P. Nikishkov / Advances in Engineering Software 34 (2003) 469–476
ambiguity, we selected the following pair of connectivity numbers: the first connectivity number is the minimum global node number for the face and the second connectivity number is the global number of node that is diagonal to the first node. 3.3. Continuous field of secondary results Secondary results (stresses in elasticity problems) can be determined at any point inside the finite element. Since stresses are combinations of displacement derivatives, then stress field is discontinuous across element boundaries. Stresses have best precision at 2 £ 2p£ffiffi 2 integration points with local coordinates j; h; z ¼ ^1= 3 as shown in Fig. 3. In order to build a continuous field of secondary results, first, it is necessary to extrapolate result values from 2 £ 2 £ 2 integration points to vertices of 20-node element. Results are calculated at eight integration points, and a trilinear extrapolation in the local coordinate system j; h; z is used fi ¼ LiðmÞ fðmÞ
ð6Þ
where fðmÞ are known function values at reduces integration points; fi are function values at vertex nodes and LiðmÞ is the symmetric extrapolation matrix: 2
LiðmÞ
6 6 6 6 6 6 6 6 6 6 ¼6 6 6 6 6 6 6 6 6 4
A
B
C
B
B
C
D
A
B
C
C
B
C
A
B
D
C
B
A
C
D
C
A
B
C
A
B A
C
3
7 D7 7 7 C7 7 7 7 B7 7 7 B7 7 7 7 C7 7 7 B7 5
pffiffi A ¼ ð5 þ 3Þ=4; pffiffi C ¼ ð 3 2 1Þ=4;
pffiffi B ¼ 2ð 3 þ 1Þ=4; pffiffi D ¼ ð5 2 3Þ=4
Secondary results are extrapolated from integration points to all nodes of elements. Values for midside nodes can be calculated as an average between values for two vertex nodal values. Then averaging of contributions from the neighboring finite elements is performed for all nodes of the finite element model. Averaging produces a continuous field of secondary results specified at nodes of the model with quadratic variation inside finite elements. Later, the results can be interpolated to any point inside element or on its surface using relation (1).
4. Subdivision of quadratic faces Subdivision of quadratic element faces depends on two factors: curvature of the face and range of result function over the face. We want to create compatible mesh of triangular elements while it is desirable to generate mesh locally, i.e. considering one face at a time. The latter means that number of subdivisions along the face side should be same in elements sharing this side. Thus number of side subdivisions should be based only on nodal values of coordinates and results for this side. 4.1. Curvature-based subdivision
ð7Þ
A
Fig. 3. Reduced 2 £ 2 £ 2 integration points (1)–(8) and vertices 1–8 of the 20-node finite element.
Curvature-based subdivision depends on the spatial locations of three nodes, which define the side. Curvature radius can be calculated using one-dimensional shape functions and coordinates of three nodes on the side of a finite element face or the side of a boundary element. Since we do not need precise value of the curvature radius, its estimate can be simplified as shown in Fig. 4. Consider a triangle with vertices at nodes 1, 2 and 3 and sides a1 ; a2 and a3 : In the ordinary finite and boundary elements, the midside node 2 is equally distant from corner nodes 1 and 3. The curvature radius R can be
Fig. 4. Approximate calculation of a curvature of the FEM/BEM side.
G.P. Nikishkov / Advances in Engineering Software 34 (2003) 469–476
estimated as: a1 R¼ 2 sin a
473
ð8Þ
Sine of the angle a can be found using dual representation of the area of triangle 123 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi s ¼ pðp 2 a1 Þðp 2 a2 Þðp 2 a3 Þ; ð9Þ 1 p ¼ ða1 þ a2 þ a3 Þ 2 and s¼
1 a a sin a 2 1 3
ð10Þ
It is adopted that the number of geometry subdivisions ng for the side is proportional to the curvature parameter r
2s 1 1 ng ¼ kg r; r ¼ a1 =R ¼ þ ð11Þ a1 a2 a3 where kg is the empirical geometry coefficient. 4.2. Results-based subdivision We are going to use texture interpolation for creating contours of scalar results on FEM/BEM surfaces. The range of results for the finite element model is defined by the minimum and maximum result values fmin and fmax : The number of subdivisions on the side should be proportional to the range of results on this side. If result values at side nodes of the element face are denoted as f1 ; f2 and f3 (Fig. 4), then the number of side subdivisions nc due to color change is determined by the relation nc ¼ k c
max{absðf2 2 f1 Þ; absðf3 2 f2 Þ} fmax 2 fmin
ð12Þ
where kc is the empirical coefficient. The final number of subdivisions on the particular face side is selected as maximum of numbers of side subdivision due to geometry curvature ng and due to result range nc : 4.3. Delaunay triangulation Delaunay triangulation [6] is used to subdivide element faces into triangles. The triangulation is performed in the local coordinate system j; h (Fig. 2b). Then local coordinates of triangle vertices are transformed to the global coordinate system using relations (1) with the shape functions (3). To apply Delaunay triangulation algorithm, it is necessary to have locations of both boundary and interior points. Locations of boundary points are defined by dividing element face sides into equal parts using number of subdivisions obtained with Eqs. (11) and (12). To seed interior points, the following procedure is used. Boundary points located at the first side (except points at vertices) are copied to interior with interval, which is equal to the distance between points at the side. The copying is
Fig. 5. Delaunay triangulation of the element face in the local coordinate system.
stopped at the face center filling half of the face. The same procedure is performed sequentially for all sides. But starting from the second side, the actual placement of a point takes place only if this point has a distance from any existing point, which is larger than the specified distance. The specified distance is computed depending on the number of side subdivisions. After completion of the seeding process, the Delaunay triangulation is performed. Example of Delaunay triangulation in the local coordinate system is given in Fig. 5 for the side subdivisions 8, 3, 2, 4.
5. Using Java 3D API for visualization The Java 3D API [5] is a set of several Java packages designed for easier development of applications and applets with three-dimensional graphics capabilities. Java 3D is an extension to the Java 2 SDK. A graphical code contains Java 3D objects, which compose a scene graph. Java 3D provides high-level programming interface by hiding implementation details. The developer specifies geometry of visual objects, their appearance and behavior and light sources as Java 3D objects. After compiling, the scene is rendered automatically with ‘quasi’-photographic quality. The latter means that the effects of light source shading are shown, however, visual object shades and reflections are ignored. The highest in hierarchy object is called a virtual universe in Java 3D. The structure of the virtual universe is described by a scene graph. A Java 3D program creates instances of Java 3D objects and places them into a scene graph data structure. The scene graph is an arrangement of 3D objects in a tree structure that completely specifies the content of a virtual universe, and how it is to be rendered. Example of a scene graph that was used in producing graphics images in this paper is shown in Fig. 6. The scene object is an instance of the BranchGroup class. It serves as a root of scene graph. The following
474
G.P. Nikishkov / Advances in Engineering Software 34 (2003) 469–476
Fig. 6. Scene graph.
objects are attached to the root: background, scale and lights. The scale object is an instance of the TransformGroup class. The scale object is used for the initial scaling of models. With the background object, it is easy to specify background color of the scene. Using Background constructor or correspondent methods, it is also possible to set geometry or image as the scene background. Both an ambient light and directional lights are used to illuminate objects. The classes MouseRotate, MouseZoom and MouseTranslate are employed for rotating, zooming and translating the scene visual objects on the screen. The capability allowing writing the object’s transform information should be specified so that the behavior code can modify this information at runtime. It is worth noting that with Java 3D, the implementation of basic behavior features is very simple. The geometry and appearance objects contain the whole information about visual objects. The TriangleArray class constructor creates an empty geometry object. The number of triangle vertices and the format of information for each vertex are specified. Both vertex coordinates and vertex normals should be supplied. The Appearance object allows
one to specify various attributes, which define rendering of the visual objects. We use the PolygonAttributes class to define that back surfaces should not be rendered. This is suitable for the FEM model. For the BEM model of the crack, both sides of surfaces should be rendered. This can be done by specifying that vertex normals of backfacing polygons are flipped and that neither front nor back surfaces are discarded. The Material class is used to set material colors (ambient, diffuse, specular and emissive) and material shininess. The material object is also used to enable lighting. The View object contains all parameters needed in rendering a three-dimensional scene from one viewpoint. It is the main Java 3D object for managing the Java 3D viewing model. Numerous view policies can be specified. For example, it is possible to set parallel or perspective projection, scene antialiasing, etc. In simple cases, it is easy to set up a user environment with the class SimpleUniverse. This utility class creates all the objects, which are necessary to view the scene graph, such as a locale, a ViewingPlatform, and a Viewer object. For many Java 3D programs, the SimpleUniverse provides all necessary functionality with just several lines of the source code. When visualizing results, the coordinates for onedimensional color pattern (texture) should be specified at vertices in addition to vertex coordinates and vertex normals. It is more convenient to produce one-dimensional color patterns inside a visualization code instead of using previously prepared images. RGB piecewise linear functions shown in Fig. 7a are used for a mapping of result values to their color representation. It can be seen that a color composed of RGB components has five intervals, each of which has the length 0:2ðfmax 2 fmin Þ=fmax : Here fmax and fmin are maximum and minimum result values for the whole FEM/BEM model. Using RGB piecewise linear functions and specifying the number of color intervals, it is easy to produce one-dimensional patterns (textures with the width
Fig. 7. (a) Functions for mapping results to color representation. Each color has a range from 0 to 1. Values fmax and fmin are maximum and minimum result values for the whole FEM/BEM model. (b) One-dimensional color patterns.
G.P. Nikishkov / Advances in Engineering Software 34 (2003) 469–476
475
Fig. 8. Mises equivalent stress for an FEM model of a valve using textures with (a) 16 colors and (b) 256 colors; (c) element face subdivision.
of one pixel) shown in Fig. 7b. Continuous color interpolation pattern can be obtained by specifying number of color intervals equal to the pattern size in pixels. Besides the geometry object for faces, we use the class LineArray for element boundaries and the class PointArray for nodes of the model. At most, six geometry objects are used: faces, element boundaries and nodes for both FEM model and BEM model.
fields are obtained. Fig. 8a and b show graphical representation of the Mises equivalent stress using textures with 16 and 256 colors. It can be seen that even for such sparse mesh, secondary results can be visualized with good quality. Triangular subdivision of finite element faces based both on curvature and on function ranges in elements is presented in Fig. 8c. 6.2. FEM/BEM model of cube containing circular crack
6. Examples In order to illustrate the developed algorithms and Java 3D code several examples of visualizing FEM/BEM models and result fields are presented. 6.1. Finite element model of a valve Finite element model of a valve is composed of 454 brick-type 20-node elements. For the considered geometry, the mesh is very sparse. The elasticity problem for the valve model is solved and displacements at nodes and stresses at reduced integration points are computed. Using element extrapolation and node averaging, the continuous stress
Analysis of a cube containing embedded circular crack at its center is performed using SGBEM/FEM alternating method. Continuous fields of finite element stresses are obtained using element extrapolation and node averaging. Boundary element stresses are computed directly at finite element nodes by numerical integration over boundary element crack model. Thus, precise enough stress values are specified on a sparse mesh. Fig. 9a and b show Mises equivalent stress for a section of the finite element model and the triangular subdivision used for creating contours using patterns with 16 and 256 colors. Four elements located near the crack front demonstrate the possibility of catching a complicated field inside one element face with quadratic function
Fig. 9. Stress field for a section of a finite element model of a cube with a circular crack modeled by boundary elements: (a) using texture with 16 colors; (b) using texture with 256 colors; (c) element face subdivision.
476
G.P. Nikishkov / Advances in Engineering Software 34 (2003) 469–476
interpolation. Triangular subdivision of the finite element faces is depicted in Fig. 9c.
7. Conclusion In this paper, a simple and efficient approach to the visualization of FEM/BEM results as color contours is presented. The technique is based on irregular subdivision of higher-order element faces into triangular elements and on using texture interpolation to produce contours. The subdivision density is determined using both face edge curvatures and ranges of result fields over element edges. After determining number of subdivisions along element edges, control points are seeded inside the element face and the face is subdivided into triangles using Delaunay triangulation procedure. Contours of results are generated by texture interpolation. The Java 3D API is used for implementation of the proposed algorithm. Examples demonstrate that contour plots of high quality can be produced even for sparse meshes of finite or boundary elements. Using the proposed algorithm and Java 3D, a customized lightweight visualization software can be easily developed.
Acknowledgements The author would like to acknowledge the helpful discussion with Prof. N.N. Mirenkov at the University of Aizu.
References [1] Nikishkov GP, Park JH, Atluri SN. SGBEM–FEM alternating method for analyzing 3D non-planar cracks and their growth in structural components. Comput Model Engng Sci 2001;2:401–22. [2] Moaveni S. Finite element analysis: theory and application with ANSYS. Englewood Cliffs, NJ: Prentice-Hall; 1999. [3] Schroeder W, Martin K, Lorensen B. The visualization toolkit: an object-oriented approach to 3D graphics. Englewood Cliffs, NJ: Prentice-Hall; 1998. [4] Upson C, Faulhaber Jr. T, Kamins D, Laidlaw D, Schlegel D, Vroom J, Gurwitz R, Dam van A. The application visualization system: a computational environment for scientific visualization. IEEE Comput Graph Appl 1989;9:30–42. [5] Sowizral H, Rushforth K, Deering M. The Java 3D API specification. Reading, MA: Addison-Wesley; 2000. [6] Du C. An algorithm for automatic Delaunay triangulation of arbitrary planar domains. Adv Engng Software 1996;27:21–6.
Generating contours on FEM/BEM higher-order surfaces using Java 3D textures G.P. Nikishkov* The University of Aizu, Aizu-Wakamatsu City, Fukushima 965-8580, Japan Received 10 December 2002; revised 24 February 2003; accepted 20 March 2003
Abstract An efficient technique to visualize primary and secondary results for combined finite element method/boundary element method models as contours is presented. The technique is based on dividing higher-order surfaces into triangles and on using texture interpolation to produce contour plots. Since results of high accuracy with significant gradients can be obtained using sparse meshes of boundary elements and finite elements, special attention is devoted to element face subdivision. Subdivision density is defined on the basis of both face edge curvature and ranges of result fields over element faces. Java 3D API is employed for code development. q 2003 Elsevier Science Ltd. All rights reserved. Keywords: Finite element method; Boundary element method; Visualization; Contours; Java 3D
1. Introduction The finite element method (FEM) and the boundary element method (BEM) are used for computational modeling for several decades. In the beginning, simple elements with linear approximation of sought functions were employed. Later, simple elements have been replaced by higher-order elements with quadratic (or even cubic) approximation due to their computational efficiency. At present, the most used elements in three-dimensional analysis are hexahedral finite element with 20-nodes and quadrilateral boundary element with 8-nodes. Both elements are built on quadratic interpolating functions. For some specific problems, it is advantageous to apply a combination of FEM and BEM. In Ref. [1], the combined FEM –BEM technique is used to analyze non-planar cracks in three-dimensional components. The FEM is used for stress analysis of the uncracked structural component and the symmetric Galerkin boundary element method (SGBEM) is employed to model the crack. According to the superposition principle, the solution is sought in an iteration procedure with alternating between the FEM and SGBEM model. Using quadratic elements for both FEM and BEM models and taking into account the ability of the BEM * Tel.: þ81-242-37-2649; fax: þ 81-242-37-2706. E-mail address: [email protected] (G.P. Nikishkov).
to model singular crack fields efficiently, it is possible to obtain displacement and stress fields with high precision even for sparse finite element and boundary element meshes. Visualization techniques for FEM and BEM models and analysis results have been developed during last decades. Commercial finite element codes such as ANSYS [2] usually have postprocessors with various options for result visualization. The Visualization Toolkit VTK [3] and the AVS system [4] provide many opportunities for field visualization. From this point of view, a researcher working with the FEM or BEM has some choice of tools for visualization of FEM/BEM results. On the other hand, available visualization tools cannot satisfy all researcher’s needs especially in cases when some non-standard methods used (such as SGBEM– FEM) or when the visualization should accompany calculations in real time. In these cases, specialized lightweight software may be useful to develop. This software can take into account specific features of the computational procedure and the results. The development of the customized visualization software can be attractive if the development cost is reasonable. Previously, when graphics libraries supported just simple graphical primitives, the development of visualization software was a complicated task. Recently developed Java 3D application programming interface (API) [5] makes the development of visualization software simpler. The Java 3D
0965-9978/03/$ - see front matter q 2003 Elsevier Science Ltd. All rights reserved. doi:10.1016/S0965-9978(03)00052-8
470
G.P. Nikishkov / Advances in Engineering Software 34 (2003) 469–476
API contains Java classes for a sophisticated threedimensional graphics rendering. The programmer works with high-level interface for creating and manipulating 3D geometric objects. Despite all its wide functionality, the API can be easily learned and used. The Java 3D employs linear triangular and quadrilateral filled polygons for surface representation. Because of this, the visualization of FEM/ BEM models consisting of simplest elements is almost straightforward. However, for higher-order elements, the transformation of element faces into triangular polygons should be done carefully taking into account both geometry and result fields. In this paper, a simple and at the same time efficient technique to visualize primary and secondary FEM/BEM result fields as color contours is presented. The technique is based on dividing higher-order surfaces into triangles and on using texture interpolation to produce contour plots. Since results of high accuracy with significant gradients can be obtained using sparse meshes of boundary or finite elements, special attention should be devoted to element face subdivision. We define subdivision density on the basis of both face edge curvature and ranges of result fields over element faces. After determining number of subdivisions along element face sides, control points are seeded inside the element face and then Delaunay triangulation is used to produce subdivision. Contours of results are not generated explicitly but obtained as a result of texture interpolation. The developed Java 3D code is illustrated with examples of visualization of FEM results for a valve and FEM/BEM results for a cube with an embedded circular crack.
2. Overview of the algorithm A computational procedure of the combined SGBEM– FEM method for fracture mechanics analysis is based on the superposition principle. According to this principle, a problem for a finite body with a crack is replaced by a superposition of two problems. The first problem for a finite body with a crack is solved by the FEM. The second problem for a crack in an infinite medium is modeled using SGBEM. For a correct superposition corresponding to the solution for a finite body with a crack, fictitious forces on the boundary of the finite element model should be found in order to compensate for the stresses caused by the presence of a crack in an infinite body. The alternating method [1] provides for an efficient solution, without assembling the joint SGBEM– FEM matrix. The SGBEM – FEM alternating method alternates between the finite element solution for an uncracked body and the boundary element solution for a crack in an infinite body. Usually, several iterations are sufficient for convergence. It should be noted that because of ability of the SGBEM to efficiently model singularities, displacements and stresses of high precision might be obtained using relatively sparse meshes of finite elements
Fig. 1. Texture interpolation inside a triangle.
and boundary elements. It the vicinity of the crack, stress fields are characterized by high gradients. The input data for the visualization consists of a set of nodes defined by spatial coordinates, a set of elements that is specified by nodal connectivities, and a set of result values. Primary results are obtained at nodes after solution of the global equation system. Secondary results, which are expressed through derivatives of the primary results, usually have the best precision at some points inside elements. In elasticity problems, the primary results are displacements, and the secondary results are strains and stresses. The visualization algorithm consists of the following main steps. 1. Obtain continuous field of results by extrapolation from integration points inside elements to element nodes with subsequent averaging (omit this step for primary results defined at nodal points). 2. Create the surface of the finite element model (omit this step for boundary element meshes). 3. Subdivide curved faces into flat triangles on the basis of face curvature and gradient of results. 4. Create contour pictures inside triangles by specifying texture coordinates of one-dimensional color pattern at triangle vertices. The last step of the visualization algorithm is performed with the use of texture interpolation available in Java 3D API. Instead of explicit contour drawing, we prepare onedimensional color texture and specify texture coordinates at triangle vertices as shown in Fig. 1.
3. Geometry and results data It is supposed that the finite element model consists of brick-type 20-node elements shown in Fig. 2a and the boundary element model is composed of 8-node quadrilateral
G.P. Nikishkov / Advances in Engineering Software 34 (2003) 469–476
471
node i; j; h; z are local coordinates ð21 # j; h; z # 1Þ; ji ; hi ; zi are nodal values of the local coordinates. Faces of 20-node finite elements have the same geometry as 8-node boundary elements (Fig. 2b). Quadratic two-dimensional shape functions are employed for coordinate and result interpolation inside 8-node quadrilateral element: Ni ¼ ð1 þ j0 Þð1 þ h0 Þðj0 þ h0 2 1Þ=4 Fig. 2. (a) Brick-type 20-node finite element and (b) 8-node boundary element.
for corner nodes Ni ¼ ð1 2 j2 Þð1 þ h0 Þ=2
for ji ¼ 0
2
elements (Fig. 2b). It can be easily noticed that the finite element faces are geometrically same as the boundary elements. Therefore, the visualization of results for the surface of the finite element model and for the boundary element model can be performed in a similar way. The finite element model is defined by two arrays: nodal coordinate array and connectivity array. Primary results (such as displacements in elasticity problems) are related to nodes. Secondary results (such as stresses depending on derivatives of primary results) have the most precise values at reduced integration Gaussian points 2 £ 2 £ 2 inside finite elements. For model and results visualization with Java 3D the following geometry information should be prepared: 1. element faces located at the surface of the finite element model and divided into triangles of appropriate size; 2. surface normal vectors for triangle vertices; 3. surface element edges divided into line segments of appropriate size; 4. texture coordinates for the one-dimensional color pattern interpolation (defined by results values at vertices of triangles).
ð3Þ
Ni ¼ ð1 2 h Þð1 þ j0 Þ=2
j0 ¼ jji ;
for hi ¼ 0
h0 ¼ hhi
A normal to the surface at any point with the specified local coordinates is calculated as a vector product of two tangent vectors along j and h directions: n~ ¼ ~e1 £ ~e2 ›x 1 ~e1 ¼ ›j ›x 1 ~e2 ¼ ›h
›x 2 ›j ›x 2 ›h
›x 3 ›j ›x 3 ›h
ð4Þ
Derivatives of the global coordinates xi in respect to the local coordinates j; h; z are determined with the use of the shape functions: X ›Nm m ›x i ¼ x ›j ›j i X ›Nm m ›x i ¼ x ›h ›h i
ð5Þ
3.1. Geometry relations for 20-node element
3.2. Surface of the finite element model
Quadratic shape functions are used to interpolate coordinates and results inside 20-node finite element: X X xk ¼ Ni xik ; f ¼ Ni fi ð1Þ
To visualize the finite element model or results at the surface of the model, it is necessary to select finite element faces that belong to the surface of the finite element model (or part of the finite element model). A topology of the finite element model is described by element connectivities. Connectivities for each finite element are global (model) node numbers listed in order of local (element) node numbers. The 20-node brick-type element has six element faces. Connectivity numbers for each face can be easily extracted from element connectivity array. The surface of the finite element model is created from outer element faces. Outer faces are mentioned in the model connectivity array only once while inner faces are mentioned exactly two times. A fast procedure to select outer element faces is to characterize each element face by two connectivity numbers instead of eight. To avoid
Ni ¼ ð1 þ j0 Þð1 þ h0 Þð1 þ z0 Þðj0 þ h0 þ z0 2 2Þ=8
ð2Þ
for vertex nodes Ni ¼ ð1 2 j2 Þð1 þ h0 Þð1 þ z0 Þ=4
for ji ¼ 0
Ni ¼ ð1 2 h2 Þð1 þ j0 Þð1 þ z0 Þ=4
for hi ¼ 0
2
Ni ¼ ð1 2 z Þð1 þ j0 Þð1 þ h0 Þ=4
j0 ¼ jji ;
h0 ¼ hhi ;
for ji ¼ 0
z0 ¼ zzi
where xk are coordinates of a point inside the element; xik is the kth coordinate of the ith element node; f is the result function value at a point; fi is the result function value at
472
G.P. Nikishkov / Advances in Engineering Software 34 (2003) 469–476
ambiguity, we selected the following pair of connectivity numbers: the first connectivity number is the minimum global node number for the face and the second connectivity number is the global number of node that is diagonal to the first node. 3.3. Continuous field of secondary results Secondary results (stresses in elasticity problems) can be determined at any point inside the finite element. Since stresses are combinations of displacement derivatives, then stress field is discontinuous across element boundaries. Stresses have best precision at 2 £ 2p£ffiffi 2 integration points with local coordinates j; h; z ¼ ^1= 3 as shown in Fig. 3. In order to build a continuous field of secondary results, first, it is necessary to extrapolate result values from 2 £ 2 £ 2 integration points to vertices of 20-node element. Results are calculated at eight integration points, and a trilinear extrapolation in the local coordinate system j; h; z is used fi ¼ LiðmÞ fðmÞ
ð6Þ
where fðmÞ are known function values at reduces integration points; fi are function values at vertex nodes and LiðmÞ is the symmetric extrapolation matrix: 2
LiðmÞ
6 6 6 6 6 6 6 6 6 6 ¼6 6 6 6 6 6 6 6 6 4
A
B
C
B
B
C
D
A
B
C
C
B
C
A
B
D
C
B
A
C
D
C
A
B
C
A
B A
C
3
7 D7 7 7 C7 7 7 7 B7 7 7 B7 7 7 7 C7 7 7 B7 5
pffiffi A ¼ ð5 þ 3Þ=4; pffiffi C ¼ ð 3 2 1Þ=4;
pffiffi B ¼ 2ð 3 þ 1Þ=4; pffiffi D ¼ ð5 2 3Þ=4
Secondary results are extrapolated from integration points to all nodes of elements. Values for midside nodes can be calculated as an average between values for two vertex nodal values. Then averaging of contributions from the neighboring finite elements is performed for all nodes of the finite element model. Averaging produces a continuous field of secondary results specified at nodes of the model with quadratic variation inside finite elements. Later, the results can be interpolated to any point inside element or on its surface using relation (1).
4. Subdivision of quadratic faces Subdivision of quadratic element faces depends on two factors: curvature of the face and range of result function over the face. We want to create compatible mesh of triangular elements while it is desirable to generate mesh locally, i.e. considering one face at a time. The latter means that number of subdivisions along the face side should be same in elements sharing this side. Thus number of side subdivisions should be based only on nodal values of coordinates and results for this side. 4.1. Curvature-based subdivision
ð7Þ
A
Fig. 3. Reduced 2 £ 2 £ 2 integration points (1)–(8) and vertices 1–8 of the 20-node finite element.
Curvature-based subdivision depends on the spatial locations of three nodes, which define the side. Curvature radius can be calculated using one-dimensional shape functions and coordinates of three nodes on the side of a finite element face or the side of a boundary element. Since we do not need precise value of the curvature radius, its estimate can be simplified as shown in Fig. 4. Consider a triangle with vertices at nodes 1, 2 and 3 and sides a1 ; a2 and a3 : In the ordinary finite and boundary elements, the midside node 2 is equally distant from corner nodes 1 and 3. The curvature radius R can be
Fig. 4. Approximate calculation of a curvature of the FEM/BEM side.
G.P. Nikishkov / Advances in Engineering Software 34 (2003) 469–476
estimated as: a1 R¼ 2 sin a
473
ð8Þ
Sine of the angle a can be found using dual representation of the area of triangle 123 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi s ¼ pðp 2 a1 Þðp 2 a2 Þðp 2 a3 Þ; ð9Þ 1 p ¼ ða1 þ a2 þ a3 Þ 2 and s¼
1 a a sin a 2 1 3
ð10Þ
It is adopted that the number of geometry subdivisions ng for the side is proportional to the curvature parameter r
2s 1 1 ng ¼ kg r; r ¼ a1 =R ¼ þ ð11Þ a1 a2 a3 where kg is the empirical geometry coefficient. 4.2. Results-based subdivision We are going to use texture interpolation for creating contours of scalar results on FEM/BEM surfaces. The range of results for the finite element model is defined by the minimum and maximum result values fmin and fmax : The number of subdivisions on the side should be proportional to the range of results on this side. If result values at side nodes of the element face are denoted as f1 ; f2 and f3 (Fig. 4), then the number of side subdivisions nc due to color change is determined by the relation nc ¼ k c
max{absðf2 2 f1 Þ; absðf3 2 f2 Þ} fmax 2 fmin
ð12Þ
where kc is the empirical coefficient. The final number of subdivisions on the particular face side is selected as maximum of numbers of side subdivision due to geometry curvature ng and due to result range nc : 4.3. Delaunay triangulation Delaunay triangulation [6] is used to subdivide element faces into triangles. The triangulation is performed in the local coordinate system j; h (Fig. 2b). Then local coordinates of triangle vertices are transformed to the global coordinate system using relations (1) with the shape functions (3). To apply Delaunay triangulation algorithm, it is necessary to have locations of both boundary and interior points. Locations of boundary points are defined by dividing element face sides into equal parts using number of subdivisions obtained with Eqs. (11) and (12). To seed interior points, the following procedure is used. Boundary points located at the first side (except points at vertices) are copied to interior with interval, which is equal to the distance between points at the side. The copying is
Fig. 5. Delaunay triangulation of the element face in the local coordinate system.
stopped at the face center filling half of the face. The same procedure is performed sequentially for all sides. But starting from the second side, the actual placement of a point takes place only if this point has a distance from any existing point, which is larger than the specified distance. The specified distance is computed depending on the number of side subdivisions. After completion of the seeding process, the Delaunay triangulation is performed. Example of Delaunay triangulation in the local coordinate system is given in Fig. 5 for the side subdivisions 8, 3, 2, 4.
5. Using Java 3D API for visualization The Java 3D API [5] is a set of several Java packages designed for easier development of applications and applets with three-dimensional graphics capabilities. Java 3D is an extension to the Java 2 SDK. A graphical code contains Java 3D objects, which compose a scene graph. Java 3D provides high-level programming interface by hiding implementation details. The developer specifies geometry of visual objects, their appearance and behavior and light sources as Java 3D objects. After compiling, the scene is rendered automatically with ‘quasi’-photographic quality. The latter means that the effects of light source shading are shown, however, visual object shades and reflections are ignored. The highest in hierarchy object is called a virtual universe in Java 3D. The structure of the virtual universe is described by a scene graph. A Java 3D program creates instances of Java 3D objects and places them into a scene graph data structure. The scene graph is an arrangement of 3D objects in a tree structure that completely specifies the content of a virtual universe, and how it is to be rendered. Example of a scene graph that was used in producing graphics images in this paper is shown in Fig. 6. The scene object is an instance of the BranchGroup class. It serves as a root of scene graph. The following
474
G.P. Nikishkov / Advances in Engineering Software 34 (2003) 469–476
Fig. 6. Scene graph.
objects are attached to the root: background, scale and lights. The scale object is an instance of the TransformGroup class. The scale object is used for the initial scaling of models. With the background object, it is easy to specify background color of the scene. Using Background constructor or correspondent methods, it is also possible to set geometry or image as the scene background. Both an ambient light and directional lights are used to illuminate objects. The classes MouseRotate, MouseZoom and MouseTranslate are employed for rotating, zooming and translating the scene visual objects on the screen. The capability allowing writing the object’s transform information should be specified so that the behavior code can modify this information at runtime. It is worth noting that with Java 3D, the implementation of basic behavior features is very simple. The geometry and appearance objects contain the whole information about visual objects. The TriangleArray class constructor creates an empty geometry object. The number of triangle vertices and the format of information for each vertex are specified. Both vertex coordinates and vertex normals should be supplied. The Appearance object allows
one to specify various attributes, which define rendering of the visual objects. We use the PolygonAttributes class to define that back surfaces should not be rendered. This is suitable for the FEM model. For the BEM model of the crack, both sides of surfaces should be rendered. This can be done by specifying that vertex normals of backfacing polygons are flipped and that neither front nor back surfaces are discarded. The Material class is used to set material colors (ambient, diffuse, specular and emissive) and material shininess. The material object is also used to enable lighting. The View object contains all parameters needed in rendering a three-dimensional scene from one viewpoint. It is the main Java 3D object for managing the Java 3D viewing model. Numerous view policies can be specified. For example, it is possible to set parallel or perspective projection, scene antialiasing, etc. In simple cases, it is easy to set up a user environment with the class SimpleUniverse. This utility class creates all the objects, which are necessary to view the scene graph, such as a locale, a ViewingPlatform, and a Viewer object. For many Java 3D programs, the SimpleUniverse provides all necessary functionality with just several lines of the source code. When visualizing results, the coordinates for onedimensional color pattern (texture) should be specified at vertices in addition to vertex coordinates and vertex normals. It is more convenient to produce one-dimensional color patterns inside a visualization code instead of using previously prepared images. RGB piecewise linear functions shown in Fig. 7a are used for a mapping of result values to their color representation. It can be seen that a color composed of RGB components has five intervals, each of which has the length 0:2ðfmax 2 fmin Þ=fmax : Here fmax and fmin are maximum and minimum result values for the whole FEM/BEM model. Using RGB piecewise linear functions and specifying the number of color intervals, it is easy to produce one-dimensional patterns (textures with the width
Fig. 7. (a) Functions for mapping results to color representation. Each color has a range from 0 to 1. Values fmax and fmin are maximum and minimum result values for the whole FEM/BEM model. (b) One-dimensional color patterns.
G.P. Nikishkov / Advances in Engineering Software 34 (2003) 469–476
475
Fig. 8. Mises equivalent stress for an FEM model of a valve using textures with (a) 16 colors and (b) 256 colors; (c) element face subdivision.
of one pixel) shown in Fig. 7b. Continuous color interpolation pattern can be obtained by specifying number of color intervals equal to the pattern size in pixels. Besides the geometry object for faces, we use the class LineArray for element boundaries and the class PointArray for nodes of the model. At most, six geometry objects are used: faces, element boundaries and nodes for both FEM model and BEM model.
fields are obtained. Fig. 8a and b show graphical representation of the Mises equivalent stress using textures with 16 and 256 colors. It can be seen that even for such sparse mesh, secondary results can be visualized with good quality. Triangular subdivision of finite element faces based both on curvature and on function ranges in elements is presented in Fig. 8c. 6.2. FEM/BEM model of cube containing circular crack
6. Examples In order to illustrate the developed algorithms and Java 3D code several examples of visualizing FEM/BEM models and result fields are presented. 6.1. Finite element model of a valve Finite element model of a valve is composed of 454 brick-type 20-node elements. For the considered geometry, the mesh is very sparse. The elasticity problem for the valve model is solved and displacements at nodes and stresses at reduced integration points are computed. Using element extrapolation and node averaging, the continuous stress
Analysis of a cube containing embedded circular crack at its center is performed using SGBEM/FEM alternating method. Continuous fields of finite element stresses are obtained using element extrapolation and node averaging. Boundary element stresses are computed directly at finite element nodes by numerical integration over boundary element crack model. Thus, precise enough stress values are specified on a sparse mesh. Fig. 9a and b show Mises equivalent stress for a section of the finite element model and the triangular subdivision used for creating contours using patterns with 16 and 256 colors. Four elements located near the crack front demonstrate the possibility of catching a complicated field inside one element face with quadratic function
Fig. 9. Stress field for a section of a finite element model of a cube with a circular crack modeled by boundary elements: (a) using texture with 16 colors; (b) using texture with 256 colors; (c) element face subdivision.
476
G.P. Nikishkov / Advances in Engineering Software 34 (2003) 469–476
interpolation. Triangular subdivision of the finite element faces is depicted in Fig. 9c.
7. Conclusion In this paper, a simple and efficient approach to the visualization of FEM/BEM results as color contours is presented. The technique is based on irregular subdivision of higher-order element faces into triangular elements and on using texture interpolation to produce contours. The subdivision density is determined using both face edge curvatures and ranges of result fields over element edges. After determining number of subdivisions along element edges, control points are seeded inside the element face and the face is subdivided into triangles using Delaunay triangulation procedure. Contours of results are generated by texture interpolation. The Java 3D API is used for implementation of the proposed algorithm. Examples demonstrate that contour plots of high quality can be produced even for sparse meshes of finite or boundary elements. Using the proposed algorithm and Java 3D, a customized lightweight visualization software can be easily developed.
Acknowledgements The author would like to acknowledge the helpful discussion with Prof. N.N. Mirenkov at the University of Aizu.
References [1] Nikishkov GP, Park JH, Atluri SN. SGBEM–FEM alternating method for analyzing 3D non-planar cracks and their growth in structural components. Comput Model Engng Sci 2001;2:401–22. [2] Moaveni S. Finite element analysis: theory and application with ANSYS. Englewood Cliffs, NJ: Prentice-Hall; 1999. [3] Schroeder W, Martin K, Lorensen B. The visualization toolkit: an object-oriented approach to 3D graphics. Englewood Cliffs, NJ: Prentice-Hall; 1998. [4] Upson C, Faulhaber Jr. T, Kamins D, Laidlaw D, Schlegel D, Vroom J, Gurwitz R, Dam van A. The application visualization system: a computational environment for scientific visualization. IEEE Comput Graph Appl 1989;9:30–42. [5] Sowizral H, Rushforth K, Deering M. The Java 3D API specification. Reading, MA: Addison-Wesley; 2000. [6] Du C. An algorithm for automatic Delaunay triangulation of arbitrary planar domains. Adv Engng Software 1996;27:21–6.