Patent Publication Number: US-8527220-B2

Title: Vortex extraction apparatus, method, program storage medium, and display system

Description:
CROSS-REFERENCE TO RELATED APPLICATION 
     This application is based upon and claims the benefit of priority of the prior Japanese Patent Application No. 2009-285535, filed on Dec. 16, 2009, the entire contents of which are incorporated herein by reference. 
     FIELD 
     The embodiment discussed herein is related to a vortex extraction apparatus and method in which vortices in a fluid space are extracted. Furthermore, the embodiment is related to a non-transitory computer-readable storage medium that stores a vortex extraction program executed by an arithmetic apparatus that executes the program and causing the arithmetic apparatus to operate as a vortex extraction apparatus. Moreover, the embodiment is related to a vortex extraction display system including a vortex extraction apparatus and a vortex display apparatus that displays vortices extracted by the vortex extraction apparatus. 
     BACKGROUND 
     A technique to extract and visualize vortices in a fluid space is preferable for determining the behavior of a fluid in various fields including, for example, the behavior of wind around wings of an airplane, water flows around screws of a ship, and air flowing through the human trachea. 
     Various techniques to extract vortices have been proposed in order to realize visualization of the vortices. A major challenge to the visualization technique is to handle a very large amount of data in order to extract vortices accurately. Thus, a high-speed large-sized arithmetic machine that handles a large amount of data at a high speed is preferably used. The amount of calculation for vortex extraction is preferably reduced. However, a greater challenge is the level of the visualization. In association with a large amount of data used for the calculation for the vertex extraction, a large amount of data is used to express extracted vortices. Thus, the visualization process is difficult to achieve unless a computer is used which is specified so as to demonstrate the same level of performance as that of the computer used for the vortex extraction. That is, visualization with an insufficiently specified computer leads to slow responses to operations such as rotation, movement, enlargement, and reduction of displayed vortices. Thus, obtaining preferable display is difficult. Furthermore, the vortices vary momentarily. Hence, temporal variations in vortices are preferably determined in order to know the behavior of the vortices. However, owing to the large amount of data used for the vortex extraction, visualizing temporal variations in vortices is difficult. 
     When the data is decimated for visualization in order to achieve a high-speed visualization process, important phenomena may be lost by the decimation and fail to be displayed despite having spent the calculation cost to accurately extract the vortex. Thus, such phenomena may be missed. 
     The followings are documents related to the vortex extraction, documents disclosing examples of a region splitting method for allowing a large amount of data to be handled, and documents to be refer to in the description provided later.
     (1) Japanese Laid-Open Patent Publication No. 2001-132700   (2) Japanese Laid-Open Patent Publication No. 2005-309999   (3) Japanese Laid-Open Patent Publication No. 2000-339297   (4) D. Sujudi, R. Haimes, IDENTIFICATION OF SWIRLING FLOW IN 3-D VECTOR FIELDS, Paper of American Institute of Aeronautics and Astronautics, 1995   (5) T. Weinkauf, Cores of Swirling Particle Motion in Unsteady Flows, Paper of IEEE, 2007   (6) On the identification of a vortex, Jinhee Jeong and Fazle Hussain, Jaurnal of fluid mechanics, vol. 285, pp 64-94, 1995   (7) Marching cubes: A high resolution 3 D surface construction algorithm, William E. Lorensen, Harvey E. Cline, Proceedings of the 14th annual conference on Computer graphics and interactive techniques, Pages: 163-169, 1987   

     SUMMARY 
     According to an aspect of the invention, a vortex extraction apparatus includes: 
     a data acquisition section that acquires flow velocity data at points distributed in a fluid space; 
     a first numerical data calculation section that calculates, for each vertex of each of voxels into which the inside of the fluid space is split, first numerical data indicating whether the vertex is located inside or outside a vortex region in the fluid, based on the flow velocity data acquired by the data acquisition section; 
     a text data creation section that determines a center line of the vortex in the fluid based on the flow velocity data acquired by the data acquisition section and determines the size of the vortex region around the center line based on the first numerical data on each vertex calculated by the first numerical data calculation section, the text data creation section creating text data indicative of the center line and the size; 
     a polygon count reduction processing section that counts the number of polygons defining a surface of the vortex region and each including an isosurface of the first numerical data calculated by the first data calculation section and integrates the voxels with new voxels into which the fluid space is roughly split until the number of polygons is equal to or smaller than a set value, to decimate the first numerical data and thereby obtain first numerical data corresponding to each vertex of each of the new voxels; and 
     a data output section that outputs the text data created by the text data creation section and the first numerical data obtained by the polygon count reduction processing section and corresponding to each vertex of each of the new voxels. 
     The object and advantages of the invention will be realized and attained by means of the elements and combinations particularly pointed out in the claims. 
     It is to be understood that both the foregoing general description and the following detailed description are exemplary and explanatory and are not restrictive of the invention, as claimed. 
    
    
     
       BRIEF DESCRIPTION OF DRAWINGS 
         FIG. 1  is a schematic diagram illustrating an example of a system including a large-scale computer according to an embodiment of a vortex extraction apparatus according to the present invention; 
         FIG. 2  is a diagram schematically illustrating a vortex extraction program according to an embodiment of the present invention; 
         FIG. 3  is a diagram schematically illustrating an embodiment of the vortex extraction apparatus according to the present invention; 
         FIG. 4  is a flowchart schematically illustrating an embodiment of a vortex extraction method according to the present invention; 
         FIG. 5  is a diagram schematically illustrating the internal structure of the large-scale computer; 
         FIG. 6  is a flowchart illustrating steps of the vortex extraction program executed by the large-scale computer; 
         FIG. 7  is a diagram illustrating an interpolation process; 
         FIG. 8  is a diagram illustrating an example of an algorithm for determining the size of a pixel (voxel); 
         FIG. 9  is a diagram illustrating a region splitting process; 
         FIG. 10  is a diagram illustrating a vorticity calculation method; 
         FIG. 11  is a diagram illustrating a vortex region extraction process; 
         FIG. 12  is a flowchart illustrating the detailed flow of a vortex center segment creation process (S 28 ) in  FIG. 6 ; 
         FIG. 13  is a schematic diagram of a voxel; 
         FIG. 14  is a diagram illustrating a rotation determination algorithm; 
         FIG. 15  is a diagram illustrating an additive rotation determination algorithm; 
         FIG. 16  is a diagram illustrating a polygon count reduction process; 
         FIG. 17  is a diagram illustrating an example of vortex display; 
         FIG. 18  is a diagram translucently illustrating vortices as seen from a direction in which the center line extends in a lateral direction; 
         FIG. 19  is a diagram translucently illustrating vortices as seen from a direction in which the center line extends perpendicularly to the sheet of  FIG. 19 ; 
         FIG. 20  is a diagram illustrating the display screen in  FIG. 18  on which the center line and lines illustrating the vortex region are further superimposed; 
         FIG. 21  is a diagram illustrating the display screen in  FIG. 19  on which the center line and lines illustrating the vortex region are further superimposed; and 
         FIG. 22  is a diagram illustrating a modification of the present embodiment. 
     
    
    
     DESCRIPTION OF EMBODIMENTS 
     Embodiments of the present invention will be described below. 
       FIG. 1  is a schematic diagram illustrating an example of a system including a large-scale computer according to an embodiment of a vortex extraction apparatus according to the present invention. 
     Here, a large-scale computer  10 , a file server  20 , a visualization computer  30  are connected together via a communication line network  40  such as the Internet. 
     The large-scale computer  10  corresponds to an embodiment of the vortex extraction apparatus according to the present invention, and arithmetically processes a large amount of data at a high speed. 
     The file server  20  is an apparatus in which there are filed: flow velocity data on which vortex extraction calculations carried out by the large-scale computer  10  are based and data indicative of the results of the vortex extraction calculations carried out by the large-scale computer  10 . 
     The visualization computer  30  is an example of a vortex display apparatus and receives data indicative of the results of vortex extraction calculations carried out by the large-scale computer  10  to display the results on a display screen  31 . The visualization computer  30  may be a commercially available personal computer (hereinafter referred to as the “PC”). 
       FIG. 1  illustrates only one visualization computer  30 . However, this represents all such computers, and multiple or many visualization computers  30  may be provided. 
     Here, it is assumed that data on the flow velocities measured at multiple points distributed in a fluid space is stored in the file server  20 , the data being indicative of the distribution of the flow velocities in the fluid space. 
     The flow velocity data may be actually measured by carrying out, for example, wind tunnel experiments or may be obtained through simulation performed by a simulation computer (not illustrated). Alternatively, the large-scale computer  10  illustrated in  FIG. 1  may also serve as a simulation computer. 
     The flow velocity data stored in the file server  20  and indicating the flow velocity distribution in the fluid space is input to the large-scale computer  10  (arrow A). The large-scale computer  10  then carries out a vortex extraction calculation. Data indicative of the result of the vortex extraction calculation performed by the large-scale computer  10  is transmitted to the file server  20  again. The data is temporarily stored in the file server  20  (arrow B). Thereafter, the data stored in the file server  20  and indicating the vortex extraction calculation result is transmitted to the visualization computer  30  in response to a request from the visualization computer  30  (arrow C). The visualization computer  30  receives the data indicative of the vortex extraction calculation result and displays vortices based on the data. 
     In the file server  20 , data on the flow velocity distribution, which varies momentarily, is accumulated in a time series manner. The large-scale computer  10  may extract vortices at each point in time. Thus, the visualization computer  30  may display momentary variations in vortices (occurrence, extinction, confluence, branching, and the like). Furthermore, the vortex at a certain point in time may be displayed such that, for example, the vortex is redirected, enlarged, or reduced. 
       FIG. 2  is a diagram schematically illustrating a vortex extraction program according to an embodiment of the present invention. 
     A vortex extraction program  100  illustrated in  FIG. 2  is installed in and executed by the large-scale computer  10 . The vortex extraction program  100  allows the large-scale computer  10  to operate as an example of the vortex extraction apparatus according to the present invention. The vortex extraction program  100  includes a data acquisition section  110 , an interpolation operation section  120 , a first numerical data calculation section  130 , a second numerical data calculation section  140 , a text data creation section  150 , a polygon count reduction processing section  160 , and a data output section  170 . The sections  110  to  170  are program sections included in the vortex extraction program  100 . The operations of the sections  110  to  170  will be described later. 
       FIG. 3  is a diagram schematically illustrating an embodiment of the vortex extraction apparatus according to the present invention. 
     The vortex extraction apparatus illustrated in  FIG. 3  includes a data acquisition section  210 , an interpolation operation section  220 , a first numerical data calculation section  230 , a second numerical data calculation section  240 , a text data creation section  250 , a polygon count reduction processing section  260 , and a data output section  270 . The sections  210  to  260  illustrated in  FIG. 3  indicate functions implemented in the large-scale computer  10  when the program sections  110  to  170  included in the vortex extraction program  100  illustrated in  FIG. 2  and having the same names as those in  FIG. 3  are executed in the large-scale computer  10 . Although the same names are used in both  FIGS. 2 and 3 , the sections  110  to  170  in  FIG. 2  refer only to software, whereas the sections  210  to  270  in  FIG. 3  refer to functions implemented in the large-scale computer  10  by a combination of the program sections illustrated in  FIG. 2  and having the same names with hardware, OS (Operating System), and the like. Thus, here, the description of the operations of the sections  210  to  270  in  FIG. 3  also serves as the description of the sections  110  to  170  in  FIG. 2 . 
     The data acquisition section  210  included in the vortex extraction apparatus  200  illustrated in  FIG. 3  provides a function to acquire data on flow velocities measured at multiple points distributed in a fluid space from which vortices are to be extracted, from the file server  20  illustrated in  FIG. 1  via the communication line network  40 . 
     Furthermore, the interpolation operation section  220  executes an interpolation operation based on the flow velocity data acquired by the data acquisition section  210  to calculate flow velocity data on each vertex of each of multiple first voxels resulting from division of the fluid space into structured grids. 
     Additionally, the first numerical data calculation section  230  calculates first numerical data indicating whether each vertex of each of the multiple first voxels into which the inside of the fluid space is split is located inside or outside a fluid vortex region. The first numerical data calculation section  230  calculates the first numerical data based on the flow velocity data acquired by the data acquisition section  210 . 
     Moreover, the second numerical data calculation section  240  calculates second numerical data indicative of the intensity of vortices in each first voxel included in a vortex region determined based on the first numerical data calculated by the first numerical data calculation section  230 . The second numerical data calculation section  240  calculates the second numerical data based on the flow velocity data acquired by the data acquisition section  210 . The second numerical data calculation section  240  may calculate the second numerical data on each voxel, either only for the inside of the vortex region or for the whole fluid space without determination of whether each vertex of the voxel is located inside or outside the vortex region. 
     Furthermore, the text data creation section  250  determines the center line of vortices in the fluid based on the flow velocity data acquired by the data acquisition section  210 . Additionally, the text data creation section  250  determines the size of the vortex region around the center line based on the first numerical data on each vertex calculated by the first numerical data calculation section  230 . Then, the text data calculation section  250  creates text data indicative of the center line and size. 
     Moreover, the polygon count reduction processing section  260  counts the number of polygons defining the surface of the vortex region, such polygons including isosurfaces of the first numerical data calculated by the first data calculation section  230 . Then, the polygon count reduction processing section  260  integrates the multiple first voxels into second voxels into which the fluid space is roughly split until the number of polygons is equal to or smaller than a set value (a “set value”  261  illustrated in  FIG. 3 ). That is, the first numerical data is decimated such that only the first numerical data corresponding to the vertices of the second voxels remains. Based on the second numerical data calculated by the second numerical data calculation section  240 , the polygon count reduction processing section  260  further creates second numerical data corresponding to each of the second voxels resulting from the integration. Here, the “set value”  261  may be pre-specified for each visualization computer  30  (see  FIG. 1 ) from the visualization computer  30 . Alternatively, the “set value”  261  may be a uniform value set via the large-scale computer  10  regardless of the type of the visualization computer  30  or the like, so as to allow an average PC to carry out a visualization process at a high speed. Alternatively, the “set value”  261  may be contained in the vortex extraction program  100  illustrated in  FIG. 2 , as a constant. 
     Moreover, the data output section  270  outputs the text data created by the text data creation section  250 . Furthermore, the data output section  270  outputs the first numerical data corresponding to each vertex of each of the second voxels obtained by the polygon count reduction processing section  260 . The data output section  270  outputs not only the text data and first numerical data but also the second numerical data corresponding to the second voxels. 
     As described above, the data output by the vortex extraction apparatus  200  is temporarily stored in the file server  20  and then transmitted to the visualization computer  30 . The visualization computer  30  displays vortices on the display screen  31  based on the data. 
       FIG. 4  is a flowchart schematically illustrating an embodiment of a vortex extraction method according to the present invention. 
     The vortex extraction method illustrated in  FIG. 4  is implemented by executing the vortex extraction program  100  illustrated in  FIG. 2 , in the large-scale computer  10  illustrated in  FIG. 1 . The vortex extraction method includes a data acquisition step (S 11 ), an interpolation operation step (S 12 ), a first numerical data calculation step (S 13 ), and a second numerical data calculation step (S 14 ). The vortex extraction method further includes a text data creation step (S 15 ), a polygon count reduction processing step (S 16 ), and a data output step (S 17 ). 
     The execution contents of steps S 11  to S 17  are the same as the operations of the sections  210  to  270  of the vortex extraction apparatus  200  illustrated in  FIG. 3 . Here, duplicate descriptions are omitted. 
     Now, a further detailed embodiment will be described. Also in the embodiment described below, the general configuration of the system in  FIG. 1  is used without any change. 
       FIG. 5  is a diagram schematically illustrating the internal configuration of the large-scale computer illustrated in  FIG. 1 . 
     The large-scale computer  10  illustrated in  FIG. 5  includes multiple (in this case, four) CPUs (Central Processing Units)  11 A to  11 D that execute programs, and a memory  12  shared by the CPUs  11 A to  11 D. The multiple CPUs  11 A to  11 D serve to extract vortices in each of the regions into which the fluid space is split. In this case, the multiple CPUs carry out a vortex extraction process on the respective areas in the fluid space. This enables high-speed arithmetic processing. 
       FIG. 6  is a flowchart illustrating the steps of the vortex extraction program executed by the large-scale computer  10  illustrated in  FIGS. 1 and 5 . 
     In this case, processing starting with a data input process (S 21 ) and ending with a data output process (S 30 ) as illustrated in  FIG. 6  is executed. 
     Here, by way of example, the data input process (S 21 ), an interpolation process (S 22 ), a region splitting process (S 23 ), and the data output process (S 30 ) are carried out by the CPU  11 A, one of the multiple CPUs  11 A to  11 D illustrated in  FIG. 5 . For the other processes (S 24  to S 29 ), the other multiple CPUs  11 B to  11 D execute a calculation for the respective regions into which the fluid space is split. Flow velocity data that is momentarily varying time series data is input to the large-scale computer  10  illustrated in  FIGS. 1 and 5 . When the CPU  11 A finishes the region splitting process (S 23 ) on flow velocity data on a certain point in time, the CPUs  11 B to  11 D carry out a vortex extraction process on the respective regions for the same point in time. In the meantime, the CPU  11 A executes the data input process (S 21 ) and the like on flow velocity data on the next point in time while the CPUs  11 B to  11 D are carrying out the vortex extraction process on the respective regions. When the CPUs  11 B to  11 D finish the vortex extraction process (S 23  to S 29 ) for a certain point in time, the CPU  11 A outputs the results of the extraction process thereof (S 30 ). As described above, in the large-scale computer  10 , the multiple CPUs  11 A to  11 D play the respective rolls to process the respective momentarily varying flow velocity data in parallel. 
     In the data input process (S 21 ) in the flowchart illustrated in  FIG. 6 , processing corresponding to the data acquisition section  210  according to the embodiment (represented by the vortex extraction apparatus  200  illustrated in  FIG. 3 ) is executed. That is, in the data input process (S 21 ), the large-scale computer  10  sequentially receives the flow velocity data stored in the file server  20  as time series data and loads the received flow velocity data into the vortex extraction program illustrated in  FIG. 6 . 
     Furthermore, in the interpolation process (S 22 ), the large-scale computer  10  carries out processing corresponding to the interpolation operation section  220  of the vortex extraction apparatus  200  illustrated in  FIG. 3 . 
       FIG. 7  is a diagram illustrating the interpolation process. 
     Part (A) of  FIG. 7  is a diagram illustrating an example of an unstructured grid obtained before an interpolation process. Part (B) of  FIG. 7  is a diagram illustrating the unstructured grid in Part (A) of  FIG. 7  on which a structured grid is superimposed. Part (C) of  FIG. 7  is a structured grid resulting from the interpolation process.  FIGS. 7(A) to 7(C)  each illustrate, in the center, a cross section of wings of an airplane obtained when the wings are simulated. The fluid space (in this case, the space around the periphery of the wings) is three-dimensional, but is two-dimensionally illustrated here for simplification. 
     Part (A) of  FIG. 7  illustrates a mosaic pattern including irregular polygons. The polygons included in the mosaic pattern have a larger average area (a larger average side length) toward the periphery of the pattern and have a smaller average area (a smaller average side length) toward the wings, located in the center of the pattern. Each vertex of each of the polygons included in the mosaic pattern is associated with flow velocity data indicative of the flow velocity of a fluid (for example, air) flowing through the fluid space; the flow velocity is measured at the position of each vertex. The flow velocity data is indicative of a flow velocity vector including the three-dimensional flow velocity direction and size of the vortices measured at the associated position. 
     The flow velocity data is indicative of the flow velocities measured at the multiple points distributed in the fluid space but is inconvenient to handle in this state. Thus, as illustrated in Part (B) of  FIG. 7 , in which two patterns are superimposed on each other, quadrangular (hexahedral) pixels (voxels) orderly partitioned in two directions, that is, x and y directions (in actuality, three directions, that is, x, y, and z directions) are assumed. Then, an interpolation process is carried out based on the flow velocity data associated with each vertex of each of the irregular polygons illustrated in Part (A) of  FIG. 7  to determine the flow velocity data associated with each vertex of each pixel (each voxel). Thus, as illustrated in Part (C) of  FIG. 7 , flow velocity data is determined on the points arranged at equal intervals in the fluid space. This facilitates the subsequent arithmetic process. 
     The length (the interval between the flow velocity data in the fluid space) of each side of each pixel (each voxel) is not particularly limited. However, excessively large data intervals may prevent phenomena from being extracted. On the other hand, when the data intervals are excessively small, the calculation takes a long time even with the high throughput of the large-scale computer  10 . Furthermore, however small the data intervals are, a resolution equal to or higher than that of the flow velocity of the original unstructured grid (Part (A) of  FIG. 7 ) is obtained. Thus, the size of the pixel (voxel) is preferably determined such that the pixel has an area (volume) equivalent to the average area (in three dimensions, the average volume) of a large number of polygons included in the unstructured grid (Part (A) of  FIG. 7 ). Alternatively, the size of the pixel (voxel) may be determined such that the pixel has a side length equivalent to the average side length of a large number of polygons. 
       FIG. 8  is a diagram illustrating an example of an algorithm for determining the size of the pixel (voxel). 
     The axis of abscissas in  FIG. 8  corresponds to ten portions on a logarithmic axis into which the value between the minimum and maximum values min and max of side length of the polygon appearing in the unstructured grid in Part (A) of  FIG. 7  is divided. Furthermore, the axis of ordinate in  FIG. 8  indicates the number of times a side appears in each portion when the axis of abscissas (side length) is divided into the ten pieces on the logarithmic axis. 
     The length of each side of each of the polygons included in the unstructured grid illustrated in Part (A) of  FIG. 7  is determined. In connection with the number of times each side length appears in the corresponding portion on the axis of abscissas in  FIG. 8 , the numbers in portions A and B with the minimum and maximum side lengths, respectively, are counted. The number of times the minimum side length appears in the portion A is defined as “a”. The number of times the maximum side length appears in the portion B is defined as “b”. The representative value of the minimum side length in the portion A is defined as L min . The representative value of the maximum side length in the portion B is defined as L max . The length of each side of the structured grid (voxel) illustrated in Part (C) of  FIG. 7  is defined as L. Then, the length L of each side of the structured grid (voxel) may be calculated in accordance with: 
     
       
         
           
             
               
                 
                   L 
                   = 
                   
                     
                       
                         
                           a 
                           · 
                           L 
                         
                         ⁢ 
                         
                             
                         
                         ⁢ 
                         min 
                       
                       + 
                       
                         
                           b 
                           · 
                           L 
                         
                         ⁢ 
                         
                             
                         
                         ⁢ 
                         max 
                       
                     
                     
                       ( 
                       
                         a 
                         + 
                         b 
                       
                       ) 
                     
                   
                 
               
               
                 
                   ( 
                   1 
                   ) 
                 
               
             
           
         
       
     
     For example, the flow velocity data on each vertex of the voxel with the thus calculated side length L is determined by the interpolation process. However, the way of determining the length L mentioned here is an example, and there may be adopted any of various methods of calculating, for example, the simple average of the side lengths in the unstructured grid and the average length of sides in a sample region in a part of the fluid space. For example, Part (A) of  FIG. 7  illustrates the unstructured grid in which the average area is smaller (the average side length is shorter) toward the wings, located in the center. Thus, when the vicinity of the wings is of interest, the average side length of a part of the unstructured grid which is located close to the wings may be determined and then, an interpolation process may be carried out so as to form a structured grid including sides with the determined average length. 
     The interpolation process itself is a well-known technique and will not be described below. Furthermore, in the above description, the flow velocity data on the unstructured grid is input. However, when the original flow velocity data corresponds to such a structured grid as illustrated in Part (C) of  FIG. 7 , the interpolation process (S 22 ) is undesired and thus skipped. 
     When the interpolation process (S 22 ) is finished, the region splitting process (S 23 ) and then a flow velocity gradient calculation process (S 24 ) are carried out. As described above, the data input and interpolation processes (S 21 ) and (S 22 ) in  FIG. 6  are executed by the CPU  11 A, one of the multiple CPUs  11 A to  11 D illustrated in  FIG. 5 . In contrast, the flow velocity gradient calculation process (S 24 ) is shared by the other multiple CPUs  11 B to  11 D so that the CPUs  11 A to  11 D carry out the flow velocity gradient calculation process on the respective multiple regions into which the fluid space is split. Thus, here, the region splitting process (S 23 ) is carried out. 
       FIG. 9  is a diagram illustrating region splitting. Also in this case, a two-dimensional structured grid is illustrated for simplification. 
     Here, the fluid space D is split into three overlapping regions D 1  to D 3 . A region D 12 , one of two shaded regions D 12  and D 23 , is included in both regions D 1  and D 2 . The other region D 23  is included in both regions D 2  and D 3 . A region thus included in multiple regions is called a sleeve region. 
     The multiple (three) CPUs  11 B to  11 D illustrated in  FIG. 5  carries out, on each of the regions D 1  to D 3 , a part of the vortex extraction process which includes the local eigenvalue calculation process (S 25 ) and the subsequent processes in the flowchart in  FIG. 6 . Here, since the sleeve regions D 12  and D 23  are set as illustrated in  FIG. 9 , for the sleeve regions D 12  and D 23 , the multiple CPUs sharing the processing of the regions D 1  to D 3  including the sleeve regions D 12  and D 23  provide almost the same processing results. Thus, after the processing, the processing results for all the regions of the fluid space D may be smoothly integrated together. 
     In the description below, the region resulting from the splitting and the whole fluid space are not distinguished from each other unless otherwise specified. 
     In the flowchart in  FIG. 6 , a combination of the flow velocity gradient calculation process (S 24 ), the local eigenvalue calculation process (S 25 ), and a vortex region extraction process (S 27 ) corresponds to the processing in the first numerical data calculation section  230  in  FIG. 3 . However, the vortex region extraction process (S 27 ) is intended to improve the accuracy of vortex region extraction, and the first numerical data calculation process without the vortex region extraction process (S 27 ) may be sufficient. 
     In the flow velocity gradient calculation process (S 24 ), the flow velocity gradient of each vertex of each voxel is calculated based on the following Expression (2). 
     Here, the velocities in the x, y, and z directions in the fluid space are denoted as u, v, and w, respectively. u, v, and w are represented by “q”, and x, y, and z are represented by “r”. Furthermore, the coordinates of each vertex are representatively denoted as (i, j, k). Moreover, i, j, and k are representatively denoted as “n”, and i, j, k, and n may be omitted. 
     In this case, the flow velocity gradient of the vertex (i, j, k) is expressed by: 
     
       
         
           
             
               
                 
                   
                     q 
                     r 
                   
                   = 
                   
                     
                       
                         ∂ 
                         q 
                       
                       
                         ∂ 
                         r 
                       
                     
                     = 
                     
                       
                         
                           q 
                           
                             n 
                             + 
                             1 
                           
                         
                         - 
                         
                           q 
                           
                             n 
                             - 
                             1 
                           
                         
                       
                       
                         2 
                         ⁢ 
                         Δ 
                         ⁢ 
                         
                             
                         
                         ⁢ 
                         r 
                       
                     
                   
                 
               
               
                 
                   ( 
                   2 
                   ) 
                 
               
             
           
         
       
     
     However, Δr denotes the length of each side of each voxel in an r direction. 
     Here, processing is preferably carried out such that q r =0 on boundaries with quiescent objects such as boundaries with wings illustrated in the center of  FIG. 7  and wall surfaces defining the fluid space. 
     In this case, q r  is expressed by Expression (2) that is a single expression. However, Expression (2) denotes nine elements of a flow velocity matrix J in the following Expression (3). 
     Once the flow velocity gradient is calculated for each vertex (i, j, k) of each voxel in accordance with Expression (2), the local eigenvalue calculation process (S 25 ) is carried out. 
     In the local eigenvalue calculation process (S 25 ), a flow velocity matrix J containing, as elements, flow velocity gradients determined in accordance with Expression (2) is given as follows. 
     
       
         
           
             
               
                 
                   J 
                   = 
                   
                     ( 
                     
                       
                         
                           
                             u 
                             x 
                           
                         
                         
                           
                             u 
                             y 
                           
                         
                         
                           
                             u 
                             z 
                           
                         
                       
                       
                         
                           
                             v 
                             x 
                           
                         
                         
                           
                             v 
                             y 
                           
                         
                         
                           
                             v 
                             z 
                           
                         
                       
                       
                         
                           
                             w 
                             x 
                           
                         
                         
                           
                             w 
                             y 
                           
                         
                         
                           
                             w 
                             z 
                           
                         
                       
                     
                     ) 
                   
                 
               
               
                 
                   ( 
                   3 
                   ) 
                 
               
             
           
         
       
     
     Furthermore, S and Q are defined by: 
     
       
         
           
             
               
                 
                   S 
                   = 
                   
                     
                       1 
                       2 
                     
                     ⁢ 
                     
                       ( 
                       
                         J 
                         + 
                         
                           J 
                           T 
                         
                       
                       ) 
                     
                   
                 
               
               
                 
                   ( 
                   4 
                   ) 
                 
               
             
             
               
                 
                   Q 
                   = 
                   
                     
                       1 
                       2 
                     
                     ⁢ 
                     
                       ( 
                       
                         J 
                         - 
                         
                           J 
                           T 
                         
                       
                       ) 
                     
                   
                 
               
               
                 
                   ( 
                   5 
                   ) 
                 
               
             
           
         
       
     
     Here, J T  is a transposed matrix of J. Then, three eigenvalues λ a , λ b , and λ c  for a matrix A expressed as follows are calculated.
 
 A=S   2   +Q   2   (6)
 
     The calculated three eigenvalues λ a , λ b , and λ c  are arranged in numerical order to be λ 1 &lt;λ 2 &lt;λ 3 . Each of λ 1 , λ 2 , and λ 3  is any of λ a , λ b , and λ c . 
     Here, the eigenvalue λ 2 , corresponding to the median of the three eigenvalues λ 1 , λ 2 , and λ 3  arranged in numerical order, is extracted. λ 2  corresponds to an example of the first numerical data calculated by the first numerical data calculation section  230  in  FIG. 3 . 
     The eigenvalue λ 2  is calculated for each vertex (i, j, k) of each voxel. When the eigenvalue λ 2  calculated for a certain vertex (i, j, k) is λ 2 &lt;0, the vertex (i, j, k) is determined to be inside the vortex region. When λ 2 ≧0, the vertex (i, j, k) is determined to be outside the vortex region. 
     This vortex region determination algorithm focuses on the fact that when a broad view of a vortex is taken, the vortex makes two-dimensional motion by revolving around the center line, which extends one-dimensionally. The algorithm itself is well-known and is called a λ 2  method proposed by Jeong et al (see “On the identification of a vortex” listed above). 
     Here, when λ 2 &lt;0 is calculated for a certain vertex (i, j, k), λ 2  is associated with the vertex (i, j, k) with the value unchanged. On the other hand, when λ 2  is calculated to be λ≧0, λ 2  is substituted with λ 2 =0, which is then associated with the vertex. 
     Thus, the vortex region is roughly extracted such that a region with λ 2 &lt;0 is determined to be a vortex region, whereas a region with λ 2 =0 is located outside the vortex. 
     Then, an absolute vorticity value calculation process (S 26 ) is carried out. In the absolute vorticity value calculation process (S 26 ), an absolute value |ω| of a vorticity ω is calculated. The absolute vorticity value calculation process (S 26 ) corresponds to the processing carried out by the second numerical data calculation section  240  in  FIG. 3 . 
     Here, the vorticity ω is calculated for each voxel. 
     When the flow velocity vector of a certain voxel is defined as U, the vorticity ω is given as follows.
 
ω=∇× U   (7)
 
     Here, ∇ is an operator expressed as follows. 
     
       
         
           
             ∇ 
             
               = 
               
                 ( 
                 
                   
                     ∂ 
                     
                       ∂ 
                       x 
                     
                   
                   , 
                   
                     ∂ 
                     
                       ∂ 
                       y 
                     
                   
                   , 
                   
                     ∂ 
                     
                       ∂ 
                       z 
                     
                   
                 
                 ) 
               
             
           
         
       
     
     Here, x denotes an outer product (vector product). The use of Expression (8) results in Expression (9). 
     
       
         
           
             
               
                 
                   U 
                   = 
                   
                     ( 
                     
                       
                         u 
                         c 
                       
                       , 
                       
                         v 
                         c 
                       
                       , 
                       
                         w 
                         c 
                       
                     
                     ) 
                   
                 
               
               
                 
                   ( 
                   8 
                   ) 
                 
               
             
             
               
                 
                   
                     
                       
                         ω 
                         = 
                           
                         ⁢ 
                         
                           ∇ 
                           
                             × 
                             U 
                           
                         
                       
                     
                   
                   
                     
                       
                         = 
                           
                         ⁢ 
                         
                           ( 
                           
                             
                               
                                 
                                   ∂ 
                                   
                                     w 
                                     c 
                                   
                                 
                                 
                                   ∂ 
                                   y 
                                 
                               
                               - 
                               
                                 
                                   ∂ 
                                   
                                     v 
                                     c 
                                   
                                 
                                 
                                   ∂ 
                                   z 
                                 
                               
                             
                             , 
                             
                               
                                 
                                   ∂ 
                                   
                                     u 
                                     c 
                                   
                                 
                                 
                                   ∂ 
                                   z 
                                 
                               
                               - 
                               
                                 
                                   ∂ 
                                   
                                     w 
                                     c 
                                   
                                 
                                 
                                   ∂ 
                                   x 
                                 
                               
                             
                             , 
                             
                               
                                 
                                   ∂ 
                                   
                                     v 
                                     c 
                                   
                                 
                                 
                                   ∂ 
                                   x 
                                 
                               
                               - 
                               
                                 
                                   ∂ 
                                   
                                     u 
                                     c 
                                   
                                 
                                 
                                   ∂ 
                                   y 
                                 
                               
                             
                           
                           ) 
                         
                       
                     
                   
                 
               
               
                 
                   ( 
                   9 
                   ) 
                 
               
             
           
         
       
     
     Here, U denotes the representative flow velocity vector of the voxels, and the flow velocity data is related to each vertex of each voxel. As one technique, it is conceivable to calculate flow velocity data on the center of the voxel by the interpolation process based on the flow velocity data on each vertex and to define the calculated flow velocity data as U so that the calculation in Expression (9) is carried out. However, here, ω is calculated as follows without separating the interpolation operation from the calculation of ω based on Expression (9). Here, for easy understanding, the description will be provided based on two dimensions. 
     Expression (9) is rewritten for the two dimensions as follows. 
     
       
         
           
             
               
                 
                   
                     
                       
                         ω 
                         = 
                         
                           ∇ 
                           
                             × 
                             U 
                           
                         
                       
                     
                   
                   
                     
                       
                         = 
                         
                           
                             
                               ∂ 
                               
                                 v 
                                 c 
                               
                             
                             
                               ∂ 
                               x 
                             
                           
                           - 
                           
                             
                               ∂ 
                               
                                 u 
                                 c 
                               
                             
                             
                               ∂ 
                               y 
                             
                           
                         
                       
                     
                   
                 
               
               
                 
                   ( 
                   10 
                   ) 
                 
               
             
           
         
       
     
       FIG. 10  is a diagram illustrating a method for calculating the vorticity ω. Here, an x-y two-dimensional plane is assumed for simplification. 
     The velocities of each of the vertices (i, j, k), (i+1, j, k), (i, j+1, k), and (i+1, j+1, k) of a certain voxel (corresponding to a certain pixel herein because the two dimensions are assumed) in the x and y directions are expressed as u and v with the coordinates of the vertex (for example, i, j, and k) as illustrated in  FIG. 10 . 
     In this case, in accordance with Expression (10), the vorticity of a central point in the pixel is expressed as follows. 
     
       
         
           
             ω 
             
               
                 i 
                 + 
                 
                   1 
                   2 
                 
               
               , 
               
                 j 
                 + 
                 
                   1 
                   2 
                 
               
               , 
               k 
             
           
         
       
     
     Then, the following expression is given. 
     
       
         
           
             
               
                 
                   
                     ω 
                     
                       
                         i 
                         + 
                         
                           1 
                           2 
                         
                       
                       , 
                       
                         j 
                         + 
                         
                           1 
                           2 
                         
                       
                       , 
                       k 
                     
                   
                   = 
                   
                     
                       
                         1 
                         2 
                       
                       ⁡ 
                       
                         [ 
                         
                           
                             
                               
                                 v 
                                 
                                   
                                     i 
                                     + 
                                     1 
                                   
                                   , 
                                   
                                     j 
                                     + 
                                     1 
                                   
                                 
                               
                               - 
                               
                                 v 
                                 
                                   i 
                                   , 
                                   
                                     j 
                                     + 
                                     1 
                                   
                                 
                               
                             
                             
                               Δ 
                               ⁢ 
                               
                                   
                               
                               ⁢ 
                               x 
                             
                           
                           + 
                           
                             
                               
                                 v 
                                 
                                   
                                     i 
                                     + 
                                     1 
                                   
                                   , 
                                   j 
                                 
                               
                               - 
                               
                                 v 
                                 
                                   i 
                                   , 
                                   j 
                                 
                               
                             
                             
                               Δ 
                               ⁢ 
                               
                                   
                               
                               ⁢ 
                               x 
                             
                           
                         
                         ] 
                       
                     
                     + 
                     
                       
                         1 
                         2 
                       
                       ⁡ 
                       
                         [ 
                         
                           
                             
                               
                                 u 
                                 
                                   
                                     i 
                                     + 
                                     1 
                                   
                                   , 
                                   
                                     j 
                                     + 
                                     1 
                                   
                                 
                               
                               - 
                               
                                 u 
                                 
                                   
                                     i 
                                     + 
                                     1 
                                   
                                   , 
                                   j 
                                 
                               
                             
                             
                               Δ 
                               ⁢ 
                               
                                   
                               
                               ⁢ 
                               y 
                             
                           
                           + 
                           
                             
                               
                                 u 
                                 
                                   i 
                                   , 
                                   
                                     j 
                                     + 
                                     1 
                                   
                                 
                               
                               - 
                               
                                 u 
                                 
                                   i 
                                   , 
                                   j 
                                 
                               
                             
                             
                               Δ 
                               ⁢ 
                               
                                   
                               
                               ⁢ 
                               y 
                             
                           
                         
                         ] 
                       
                     
                   
                 
               
               
                 
                   ( 
                   11 
                   ) 
                 
               
             
           
         
       
     
     In actuality, a three-dimensional calculation in accordance with Expression (9) is executed instead of a two-dimensional calculation in accordance with Expressions (10) and (11). Thus, the vorticity ω of each voxel is calculated to determine the intensity of the vortices expressed as the absolute value |ω| thereof. The intensity |ω| of the vortices corresponds to an example of the second numerical data calculated by the second numerical data calculation section  240  in  FIG. 3 . 
     Then, in the present embodiment, the vortex region extraction process (S 27 ) in the flowchart in  FIG. 6  allows, if desired, modification of the eigenvalue λ 2  temporarily determined in the local eigenvalue calculation process (S 25 ) and indicating whether the vertex is located inside (λ 2 &lt;0) or outside (λ 2 =0) the vortex region. 
       FIG. 11  is a diagram illustrating the vortex region extraction process. 
     Part (A) of  FIG. 11  is a diagram schematically illustrating a region with λ 2 &lt;0. Part (B) of  FIG. 11  is a diagram illustrating separation of the vortex region based on the value of |ω|. Also in this case, a two-dimensional plane is assumed for simplification. 
     In Part (A) of  FIG. 11 , it is assumed that the shaded region corresponds to the region with λ 2 &lt;0 (inside the vortex region). The condition of vortices is unknown at this stage, but it is assumed that the illustrated two vortices are present. Here, for the inside of a region with an eigenvalue λ 2  calculated in the local eigenvalue calculation process (S 25 ) to be smaller than zero, the absolute vorticity value determined in the absolute vorticity value calculation process (S 26 ), that is, the vortex intensity |ω|, is checked. Then, |ω| is compared with a certain threshold. λ 2  of voxels with |ω| equal to or smaller than the threshold is replaced with λ 2 =0. For example, here, it is assumed that the voxel in the shaded region in Part (B) of  FIG. 11  has |ω| equal to or smaller than the threshold. As described above, the vortex region temporarily extracted using λ 2  is separated using |ω| if desired. This allows the vortex region to be more accurately extracted. 
     After the vortex region is extracted as described above, a process of creating a segment representing the vortex center (S 28 ) in  FIG. 6  is carried out. The process of creating a segment representing the vortex center (S 28 ) corresponds to the processing carried out by the text data creation section  250  in  FIG. 3 . 
       FIG. 12  is a flowchart illustrating the detailed flow of the process of creating a segment representing the vortex center (S 28 ). Furthermore,  FIG. 13  is a schematic diagram of a voxel. 
     As illustrated in  FIG. 13 , each vertex of the voxel is associated with flow velocity data expressed by such a vector as illustrated by an arrow in  FIG. 13 . Here, the flow velocity data on one vertex is representatively expressed as (u, v, w) (“u” denotes the flow velocity in the x direction, “v” denotes the flow velocity in the y direction, and “w” denotes the flow velocity in the z direction). 
     In the process of creating a segment representing the vortex center as illustrated in  FIG. 12 , processing described below is carried out on two surfaces (I and II surfaces) directed in the z direction, two surfaces directed in the y direction, and two surfaces directed in the x direction in this order. Here, the two surfaces (I and II surfaces) directed in the z direction will be described. The same processing is executed on the two surfaces directed in each of the y and x directions except for the direction. 
     First, for the I surface, one of the two surfaces (I and II surfaces) directed in the z direction, the flow velocity vectors of the four vertices corresponding to the four corners of the I surface are mapped to the I surface and processed. In other words, the program neglects a z direction component (w), one of the flow velocity vectors (u, v, w) of the four vertices. The program determines whether or not the fluid is rotating, based on the x and y components (u, v). A criterion for determination of whether or not the fluid is rotating will be described later. Upon determining that the fluid is rotating on the I surface, the program carries out similar determination on the II surface, the other surface directed in the same z direction. That is, the program neglects a (z) component  2 , one of the flow velocity vectors (u, v, w) of the four vertices corresponding to the four corners of the II surface. The program determines whether or not the fluid is rotating in the same direction as the rotation direction on the I surface, based on the x and y components (u, v). Moreover, the program determines whether or not outward and inward flux conditions are the same for the I surface and for the II surface. Upon determining that rotation is occurring in the same direction as that for the I surface and for the II surface and that the outward and inward flux conditions are the same for the I surface and for the II surface, the program assigns a positive constant Vc to the voxel. The constant Vc indicates that the voxel is located on the vortex center axis. The same processing is carried out on voxels in the fluid space for which at least one vertex has λ 2 &lt;0 (this is indicative of the inside of the vortex region). Moreover, the same processing is carried out on the two surfaces directed in the y direction. The same processing is carried out on the two surfaces directed in the x direction. 
       FIG. 12  is a flowchart basically illustrating the above-described processing. The contents of the processing will be described below in accordance with the flowchart. 
     In step S 31 , on which of the z, y, and x directions processing is to be carried out is determined. The above-described processing is executed on the z direction. 
     In step S 32 , whether or not the processing has been carried out on all of the three directions, namely the z, y, and x directions, is determined. When the processing has been carried out on all of the three directions, namely the z, y, and x directions, the program proceeds to step S 42 . The processing in step S 42  will be described later. In the following description, it is assumed that the z direction is set in step S 31 . In step S 33 , one of the unprocessed voxels is set to be the next processing target. In step S 34 , the program determines whether or not all the voxels have been processed with respect to the direction (here, the z direction) set in step S 31 . In step S 34 , upon determining that all the voxels have been processed with respect to the currently set direction (here, the z direction) set in step S 31 , the program returns to step S 31  to set the next direction (for example, the y direction follows the z direction). When an unprocessed voxel is set in step S 33 , the program proceeds to step S 35  via step S 34 . In step S 35 , the program determines whether or not at least one of the eight vertices of the voxel set in step S 33  meets λ 2 &lt;0. When the eight vertices all have λ2=0 (as described above, λ 2 ≧0 has been replaced with λ 2 =0), the voxel is located outside the vortex region. Thus, the processing is undesired for the voxel. The program returns to step S 33  to set the next voxel. 
     Upon determining, in step S 35 , that at least one of the eight vertices has λ 2 &lt;0, the program proceeds to step S 36  to determine whether or not for the I surface, one of the two surfaces (I and II surfaces) directed in the current set direction (here, the z direction), the fluid is rotating. In the determination, as described above, the flow velocity components, in the currently set direction (here, the z direction), of the flow velocity vectors of the four vertices corresponding to the four corners of the I surface are neglected, with the flow velocity components in the remaining two directions taken into account. The determination criterion for rotation will be described further later. 
     Upon determining that rotation is not occurring on the I surface, the program suspends the processing executed on that direction (here, the z direction) of the voxel, and returns to step S 33  (step S 37 ). Upon determining that rotation is occurring on the I surface, the program then carries out similar determination on the other surface (II surface) directed in the same direction (here, the z direction) (step S 38 ). Upon determining that rotation is not occurring on the II surface (step S 39 ) and upon determining that rotation is occurring on the II surface but that the rotation is not in the same direction as that on the I surface (the same direction also for the outward and inward fluxes described later) (step S 40 ), the program returns to step S 33 . Upon determining that the rotation on the II surface is in the same direction as that of the rotation on the I surface (the same direction also for the outward and inward fluxes described), the program proceeds to step S 41 . Then, the program sets, for the current processing target voxel, a value Vc indicating that the voxel is located on the center axis. The program then returns to step S 33 . In step S 33 , the next one of the unprocessed voxels is set to be the next processing target with respect to the direction corresponding to the current processing (here, the z direction). 
     When the above-described processing has been carried out on all of the three directions, the z, y, and x directions, the program proceeds to step S 42 . In step S 42 , in the voxels for which the value Vc, indicating that the voxel is located on the center axis, is set, a spline function is applied to a group of almost one-dimensionally arranged voxels to calculate the vortex center line. The condition for the group of almost one-dimensionally arranged voxels is used because besides the one-dimensionally arranged voxels, another secondary center line may be present away from the voxels. When the vortex center line is determined, the program then proceeds to step S 43 . Here, when a plane is assumed which crosses the center line at a right angle at a longitudinally central point on the center line, the minimum and maximum distances r min  and r max  from the center line (the intersection between the center line and the plane) to the surface of the vortex region (region with λ 2 &lt;0) on the plane are determined. Then, a region generally representing the vortex region is determined based on the minimum and maximum distances r min  and r max . A segment is then created which corresponds to extension of the region along the center line and which represents the vortex center. 
     In the process of creating a segment representing the vortex center as illustrated in  FIG. 6  (S 28 ), data on the segment is created which is based on the vortex center line and vortex region determined as described above and which represents the vortex center. The data on the segment representing the vortex center includes a function indicative of the center line and data including r min  and r max  and indicating the vortex region. The data is created as text data unassociated with the individual polygons. 
     Now, the process of determining whether or not rotation is occurring in steps S 36  and S 38  in  FIG. 12  will be described. 
       FIG. 14  is a diagram illustrating a rotation determination algorithm. Also in this case, for simplification, each voxel is illustrated as a pixel mapped onto a two-dimensional plane. 
     In both Part (A) and Part (B) of  FIG. 14 , a fluid is determined to be rotating in a shaded voxel and not to be rotating in a voxel on the immediate right of the shaded voxel. In Part (A) of  FIG. 14 , the flow velocity vectors of the four corners of the shaded pixel are directed in the same rotating direction (here, clockwise rotating direction) outside and along each side of the pixel with an angle of at most 90° to the side. In this case, the central point of an outward flux of a vortex is present in the pixel. 
     Furthermore, in Part (B) of  FIG. 14 , the flow velocity vectors of the four corners of the shaded pixel are directed inward of each side and in the same rotating direction along each side with an angle of at most 90° to the side. In this case, the central point of an inward flux of a vortex is present in the pixel. 
     In  FIGS. 14(A) and 14(B) , the pixel on the immediate right of the shaded pixel fails to meet both the outward and inward flux conditions described above. 
     Here, only one of the I and II surfaces directed in the same direction (for example, the z direction) has been described. However, if both outward and inward fluxes are rotating in the same direction in the same manner on both the I and II surfaces directed in the same direction, then for the pixel (voxel), the value Vc indicating that the voxel is located on the center line of the vortex is set. 
     However, the rotation determination algorithm described with reference to  FIGS. 14(A) and 14(B)  infrequently fails to determine the voxel actually located on the center line to be located on the center line. However, here, no problem occurs because the center line is specified through application of the spline function. That is, the spline function works continuously even though the value Vc infrequently fails to be set for the voxel on the center line. 
     However, the determination accuracy for the voxel on the center line may be increased by adding a determination algorithm described below to the determination algorithm described with reference to  FIGS. 14(A) and 14(B) . 
       FIG. 15  is a diagram illustrating an additive rotation determination algorithm. 
     In both Part (A) and Part (B) of  FIG. 15 , the fluid is determined to be rotating in the lower right pixel (voxel). However, Part (A) and Part (B) of  FIG. 15  illustrate only one of the I and II surfaces, and both the I and II surfaces preferably meet the rotation conditions as described above. 
     The lower right quadrangular pixels in both Part (A) and Part (B) of  FIG. 15  fail to meet the rotation determination conditions in both Part (A) and Part (B) of  FIG. 14 . Thus, the opposing corners of the pixel are connected together with lines to form triangles. Paired opposing corners are connected together with a line to form two triangles. Other paired opposing corners are connected together to form two other triangles. The same determination criteria as those described with reference to  FIG. 14  are applied to each of the four triangles. In other words, when the shaded triangles in Part (A) and Part (B) of  FIG. 15  are noted, the flow velocity vectors of the three vertices of each of the triangles are directed in the same rotating direction along the respective sides with an angle of at most 90° to the sides. Furthermore, the conditions for the outward and inward fluxes (both Part (A) and Part (B) of  FIG. 15  illustrate an outward flux) are met. The value Vc is also set for voxels in which such triangles are present. The determination accuracy for the voxel on the center line of the vortex may be increased by adding the determination algorithm in  FIG. 15  to the determination algorithm in  FIG. 14 . 
     Now, a polygon count reduction process (S 29 ) in the flowchart in  FIG. 6  will be described. 
       FIG. 16  is a diagram illustrating the polygon count reduction process. Also in this case, for simplification, polygons are illustrated as two-dimensional pixels. Furthermore, here, also for simplification, only 4×4 polygons (pixels), namely, a total of 16 polygons, are illustrated. 
     Each vertex of each polygon (pixel) is associated with λ 2 . Here, the subscript  2  of λ 2  is omitted, and i, j, and the like which are indicative of coordinates are illustrated as subscripts. However, the representative of all the values λ 2  of the coordinates carries the subscript  2  and is illustrated as λ 2 . As described above, each λ 2  is either λ 2 &lt;0 or λ 2 =0. Here, a negative number with a small absolute value, for example, −0.1, is set, and a Marching Cube method is applied. Thus, polygons each formed of an isosurface with λ 2 =−0.1 are generated. The polygons are expressed as segments on  FIG. 16 . The Marching Cube method itself is a well-known technique described in “Marching cubes: A high resolution 3D surface construction algorithm” listed above. The polygons each formed of an isosurface with λ 2 =−0.1 define the surface of the vortex region. In  FIG. 16 , oblique solid lines correspond to the respective polygons. In a fluid space in which three-dimensional voxels are arranged, each polygon forms a triangle. 
     The number of the polygons thus generated is counted. Here, as described with reference to  FIG. 5 , the multiple CPUs  11 B to  11 D share the region to perform the processing (see  FIG. 9 ). Thus, here, the CPU  11 A receives reports of the polygon counts in the regions D 1  to D 3 , from the CPUs  11 B to  11 D, assigned to the regions D 1  to D 3 , respectively. The CPU  11 A aggregates the polygon counts so as to avoid doubly counting the polygon counts in the sleeve regions D 12  and D 23 . If the aggregate polygon count exceeds a set value, the CPU  11 A instructs the CPUs  11 B to  11 D to roughly split the fluid space and integrate resultant new voxels with the above-described voxels until the polygon count is equal to or smaller than the set value. Specifically, the 4×4 voxels illustrated in  FIG. 16  are integrated with adjacent 2×2 voxels to obtain new voxels indicated by an alternate long and short dash line in  FIG. 16 . The integration with the new voxels means that the resolution is reduced by removing the values λ 2  of the points other than the vertices of the new voxels, with only the values λ 2  of the vertices of the new voxels left. For the new voxels thus created, isosurfaces (polygons) with λ 2 =−0.1 are created and the numbers of the polygons are aggregated, as described above. The numbers of the original voxels arranged in the x, y, and z directions are defined as m, n, and o, respectively. Then, with the numbers of the voxels sequentially reduced to, for example, m/2, n/2, and o/2; m/3, n/3, and o/3; . . . , that is, with the resolution gradually reduced, polygons are repeatedly created until the aggregate number of polygons is equal to or smaller than the set value. 
     Here, the set value means the number (for example, 1,000,000) of polygons that may be processed by the visualization computer  30  illustrated in  FIG. 1 , at an operation speed sufficient for visualization of the vortex. The set value may be pre-specified by the visualization computer  30  or may be a uniformly set value at which a standard PC preferably used as a visualization computer may process the polygons at a sufficient operation speed. 
     Once the polygon count becomes equal to or smaller than the set value, the intensity |ω| of the vorticity is set for each of the voxels leading to the polygon count. The intensity |ω| of the vorticity has been determined for each of the original, non-integrated voxels (step S 26  in  FIG. 6 ). For each of the integrated voxels, a representative value is set, such as the average value or median of the vorticity intensities |ω| of the multiple original voxels contained in the integrated voxel. 
     Such a polygon count reduction process (step S 29  in  FIG. 6 ) reduces the resolution for λ 2  and |ω| in the fluid space but enables the visualization computer  30  to provide smooth display. This allows temporal variations in vortices (occurrence, extinction, branching, confluence, and the like) to be easily observed, and improves responsiveness to rotation used to allow the vortices to be viewed at a different angle and other operations such as enlargement, reduction, and movement of the vortices. Furthermore, the display resolution decreases from the display resolution obtained during the vortex extraction operation. However, the process of creating text data indicative of the vortex center line and region is carried out before the resolution is reduced. Thus, such important transition phenomena as hindered from appearing by reducing the resolution may be appropriately sensed. 
     In the data output process (S 30 ) in  FIG. 6 , processing corresponding to the data output section  270  in  FIG. 3  is carried out. That is, here, the text data indicative of the segment representing the vortex center and created in the process of creating the segment (S 28 ) is output. Moreover, in the data output process (S 30 ), the following data is also output: λ 2  associated with each vertex of each of the new voxels resulting from the integration in the polygon count reduction process (S 29 ) and |ω| associated with each of the new voxels resulting from the integration. Furthermore, the positional coordinates of the voxels located on the vortex center line (voxels associated with the value V c ) may be output. However, in the present embodiment, data indicative of polygons created during the execution of the polygon count reduction process (S 29 ) is not output. This is because data indicative of polygons is large in amount and is disadvantageous in terms of communication time and memory capacity for storage and because polygons may be created by the visualization computer  30  at a sufficient operation speed if desired. 
     The data thus output by the large-scale computer  10  illustrated in  FIG. 1  is temporarily stored in the file server  20 . The data is then transmitted to the visualization computer  30  in response to a request from the visualization computer  30 . Based on the received data, the visualization computer  30  displays, on the display screen  31 , the segment representing the vortex center based on the text data and vortices themselves based on the voxel data. 
     Although the display mode of the visualization computer  30  is not particularly limited, typical display modes of the visualization computer  30  will be described below. 
     The minimum display includes the display of the segment representing the vortex center, based on the text data indicative of the segment. The minimum display allows the occurrence, extinction, confluence, branching, movement, and the like of vortices to be viewed. 
     Detailed display includes the display of the surface shape of the vortices on the display screen  31  based on polygon data created by the visualization computer  30  using the Marching Cube method. 
     Then, the detailed shape of the vortices, which is unclear from the display, based on the text data, of the segment representing the vortex center, may be determined. 
     Furthermore, in addition to the display of the surface shape of the vertices based on polygons, the display of |ω| associated with the voxel containing the surface polygons may be provided in a color corresponding to the magnitude of the value |ω|. For example, |ω| with a large value is displayed in red, whereas |ω| with a small value is displayed in blue. Then, the vortex intensity is displayed in terms of a color distribution, enabling the direction of the vortex center to be understood. 
     Furthermore, in addition to the display of the surface shape of the vertices based on polygons, the translucent display of the surface shape may be provided so that the vortex center line and region based on the text data may be superimposed on the surface shape. Moreover, in the above-described display modes, a cross section may be displayed so as to allow the inside of the vortices to be observed. 
     Now, display examples will be described with reference to the drawings. 
       FIG. 17  is a diagram illustrating an example of vortex display. 
     Here, a fluid is filled in a cubic. Vortices are illustrated which are extracted based on a flow velocity distribution created by simulating the case where the top surface of the cubic is moved at a constant speed. Furthermore, the magnitude of |ω| of the surface of the vortices may be expressed by shading. 
       FIG. 18  is a diagram translucently illustrating vortices as viewed from a direction in which the center line extends in the lateral direction.  FIG. 19  is a diagram translucently illustrating the vortices as viewed from a direction in which the center line extends perpendicularly to the sheet of  FIG. 19 . 
       FIGS. 18 and 19  illustrate voxels for which the value Vc corresponding to the center line is set. 
       FIGS. 20 and 21  are diagrams corresponding to the display screens in  FIGS. 18 and 19  on which the center line and lines representing the segment region representing the vortex center are further superimposed. 
     According to the present embodiment, the vortices may be displayed in such various modes at a sufficient drawing speed. Furthermore, according to the present embodiment, text data on the segment representing the vortex center is created by calculations with a high resolution obtained before the integration of the voxels. This allows possible missing of preferable information to be avoided. 
       FIG. 22  is a diagram illustrating a modification of the present embodiment. 
       FIG. 22  illustrates a large-scale computer that is the same as the large-scale computer  10  illustrated in  FIG. 5 . In the description of  FIG. 5 , data indicative of the flow velocity distribution in the fluid space is stored in the file server  20  (see  FIG. 1 ). In contrast, in  FIG. 22 , simulation for determining the flow velocity distribution in the fluid space is also performed by the large-scale computer  10 . The fluid space is split into regions assigned to the respective CPUs  11 A to  11 D. Each of the CPUs  11 A to  11 D first calculates the flow velocity distribution in the corresponding region by simulation, and then extracts vortices. When the vortex extraction is finished, the simulated flow velocity distribution calculation is carried out on the next point in time, and vortices at that point in time are then extracted. The CPUs  11 A to  11 D in  FIG. 22  alternately carry out the simulated flow velocity distribution calculation and the vortex extraction. Thus, in the present case, the large-scale computer  10  may avoid specializing in vortex extraction and may be used for the simulated flow velocity distribution calculation. 
     According to the present embodiment and the like, the text data indicative of the center line and size of the vortex region is created based on the flow velocity data with a relatively high resolution. These preferable pieces of information are created first. Therefore, even when voxel data, polygon data, and the like are reduced to perform display in an apparatus with specifications inferior to those of the vortex extraction apparatus according to the present embodiment and the like, the vortex display process may be carried out at a high speed without missing preferable information. 
     All examples and conditional language recited herein are intended for pedagogical purposes to aid the reader in understanding the invention and the concepts contributed by the inventor to furthering the art, and are to be construed as being without limitation to such specifically recited examples and conditions, nor does the organization of such examples in the specification relate to a showing of the superiority and inferiority of the invention. Although the embodiment of the present invention has been described in detail, it should be understood that the various changes, substitutions, and alterations could be made hereto without departing from the spirit and scope of the invention.