Method of Automatic Generation of a Computer for Numerical Simulation of Physical Phenomenon

The present invention relates to a method for generating a numerical calculation program in which a physical phenomenon represented by a partial differential equation is converted into a numerical expression and is schematically displayed (namely, simulated), and in particular, a method for generating a numerical calculation program suitable for simulating a distribution of physical quantity in a space, for example, for a distribution of an electromagnetic field, a thermal conduction analysis, and a fluidal analysis.

Since a numerical calculation program varies depending on the shape of a region to be simulated, an equation (a partial differential equation, in general) representing a physical rule dominating a physical phenomenon to be simulated, a condition on a boundary of the region, and so forth; a considerable amount of human power is required to generate such a program. It is therefore desirable to automatically generate the program based on a simple input data.

Conventionally, program languages and programs for generating a program have been proposed, for example, by Morris, S. M. et al, "SALEM-A Programming System for the Simulation of Systems Described by Partial Differential Equations", Proc. Fall Joint Computer Conf., vol. 33, pp. 353-357 (1968); Cardenas, A. F. et al, "A Language for Partial Differential Equations", Comm. ACM, vol. 13, No. 3, pp. 184-191 (1970.3); Rice, J. R. et al, "ELLPACK Progress and Plans in Elliptic Problem Solvers", pp. 135-162, Academic Press (1981); and IMSL Inc., "TWODEPEP; A Finite Element Program", 3rd ed. (1982).

In the prior art techniques, the form of partial differential equations and the shape of the region to be simulated are limited. To overcome this difficulty, there have been published a program language called Differential Equation Solver Language and a program for generating a numerical simulation program by use of the programming language. For example, Transactions of 26th Conference of Information Processing Society of Japan, pp. 427-428, 429-430, and 431-432 describe a method for generating a program to achieve a numerical calculation based on the finite difference calculus for partial differential equations to be simulated, and Transactions of 29th Conference of Information Processing Society of Japan, pp. 1519-1522 and 1523-1524 outline a method for generating a program to effect a numerical calculation based on the finite element method for equations to be simulated.

Details about the latter are disclosed by Japanese Patent Unexamined Publication No. 60-140433 filed by the present assignee.

In a calculation according to the finite element method, a two- or three-dimensional region to be subjected to the numerical calculation is divided into small triangular regions called elements and first-degree equations approximately equivalent to partial differential equations are determined for vertices (called nodes) of each element, thereby obtaining approximate solutions for unknown quantities at the nodes.

According to the technique disclosed by Japanese Patent Unexamined Publication No. 60-140433, when a numerical calculation program is generated from the results of the region division, data indicating the linkage information between the nodes is generated for each node and the program is generated by use of the data, thereafter the program is executed. As a result, the amount of the data increases in proportion to the number of the nodes, and the increase of data is particularly remarkable in a case of a region to be subjected to a three-dimensional simulation. This leads to a problem in that a memory having a considerably great capacity is necessary.

The conventional apparatus disclosed by the articles referenced above, however, have not disclosed a technique for automatically dividing a region into elements when the total numbers of division specified by the high level language for the opposing sides of the region are different.

Moreover, among the forms of the objective partial differential equations, only particular forms have been instantiated, namely, a technique for automatically generating a program to effect a numerical calculation for equations which vary according to the simulation model has not been disclosed.

Furthermore, the disclosed examples are associated with only particular boundary conditions, that is, a method for incorporating an arbitrary boundary condition has not been disclosed.

It is therefore an object of the present invention to provide a method for automatically generating a numerical simulation program which is capable of achieving a numerical simulation according to the finite element method based on a reduced number of data items related to the elements.

Another object of the present invention is to provide a method for automatically generating a program which performs a numerical simulation according to the finite element method when an arbitrary region is specified to be divided.

Still another object of the present invention is to provide a method for automatically generating a program which effects a numerical simulation according to the finite element method for arbitrary partial differential equations.

To this end, according to the present invention, the partial differential equations inputted with respect to an unknown quantity A are processed according to the finite element method in which simultaneous first-order equations approximately equivalent to the partial differential equations are obtained, and then a matrix equation equivalent to the simultaneous first-order equations, namely, (k.sub.ij) (a.sub.i)=(d.sub.i) is processed to derive a coefficient {a.sub.j } so as to generate a program which calculates the unknown quantity A. For each group of a plurality of elements, the program calculates a contribution determined by the positions of the nodes associated with the elements with respect to a portion of matrix element group such as k.sub.lm and k.sub.lk and a portion of constant group such as d.sub.i and d.sub.m determined by the numbers assigned to a plurality of nodes contained in the group of a plurality of elements and then generates the final values of the matrix element group {k.sub.ij } and the constant group {d.sub.i } from the contribution thus calculated for each group of a plurality of elements.

If data indicating the relationship between the elements and the nodes and the position of each node are known, this method enables to obtain the final values as described above, which does not require the conventially required data indicating for each node, other nodes which are connected to the node and a number of the element to which the above other nodes belong .

The present invention will be apparent from the following detailed description taken in conjunction with the accompanying drawings in which:

FIG. 1 is a schematic flowchart illustrating a program generation processing 3 of the present invention;

FIG. 2 is a schematic diagram illustrating an execution flow of a simulation accomplished by a program generated according to the present invention;

FIG. 3 is an explanatory diagram demonstrating an example of a region to be simulated;

FIG. 4 is an explanatory diagram illustrating an example of a shape division information 1 specifying the region of FIG. 3;

FIG. 5 is an explanatory diagram depicting an example of an equation information 2 used in the processing of FIG. 1;

FIG. 6 is an explanatory diagram depicting an example of the plurality of elements generated by the region division processing 4 of FIG. 1;

FIG. 7 is a schematic diagram demonstrating the node.element data table group 7 generated by the region division processing of FIG. 1;

FIG. 8 is a schematic diagram showing an example of the discretization table 171 generated by the discretization of equation 5 of FIG. 1;

FIGS. 9A-9D are explanatory diagrams collectively illustrating a main section of a calculation Fortran program generated by the calculation program generation processing 6 of FIG. 1;

FIG. 10 is a schematic diagram showing an example of a basis function .phi..sub.i used in the program generation processing 3 of FIG. 1;

FIGS. 11A-11C are schematic diagrams illustrating a general flow of the equation discretization processing 5 of FIG. 1;

FIG. 12 is an explanatory diagram showing an equation table 111 generated by the equation discretization processing 5 of FIG. 1;

FIG. 13 is a schematic diagram depicting a boundary condition table 112 generated by the equation discretization processing 5 of FIG. 1;

FIG. 14 is an explanatory diagram showing a boundary side table 113 generated by the equation discretization processing 5 of FIG. 1;

FIG. 15 is a schematic diagram illustrating an integration form equation table 141 generated by the equation discretization processing 5 of FIG. 1;

FIG. 16 is an explanatory diagram illustrating an integration evaluation table 161 beforehand generated to be used by the equation discretization processing 5 cf FIG. 1;

FIGS. 17A-17C are schematic diagrams collectively depicting a general flowchart of the calculation program generation processing 6 of FIG. 1;

FIG. 18A is a schematic diagram illustrating elements related to the calculation of matrix elements k.sub.ij and constant terms d.sub.i associated with a node N.sub.i ;

FIG. 18B is a schematic diagram illustrating elements for which integration is necessary to obtain a plurality of matrix elements related to three nodes N.sub.i, N.sub.j, and N.sub.k is an element;

FIG. 19 is a schematic diagram illustrating tables 190-191 generated by the calculation Fortran program 8 generated by the calculation program generation processing 6 of FIG. 1 to store the matrix elements {K.sub.ij } and constants {d.sub.i };

FIG. 20 is an explanatory diagram showing a general flowchart of the region division processing 4 of FIG. 1; and

FIG. 21 is a schematic diagram for explaining the region division processing 4 of FIG. 1.

FIG. 1 is a flowchart of the numeric calculation program generation processing 3 according to the present invention. First, the operator inputs a shape.division information 1 indicating positions of boundaries or a shape of a region to be simulated and a procedure for dividing the region and an equation information representing partial differential equations dominating a physical phenomenon and initial values of variables in the equations or values on boundaries of the region (boundary conditions) to a computer (not shown) in which the program executing above-mentioned processing is installed in advance.

FIG. 3 is a schematic diagram showing an example of a region 41 to be simulated. The region 41 comprises subregions S1 and S2, which include sides L.sub.1 -L.sub.4 and sides L.sub.3 and L.sub.5 -L.sub.7, respectively.

FIG. 4 is an explanatory diagram showing an example of the shape.division information 1 for the region 41. FIG. 3 will be described in detail later.

FIG. 5 is a schematic diagram depicting an example of the equation information 2 comprising an equation 53 and boundary conditions 54-55 associated with sides L.sub.1, L.sub.6, etc. and sides L.sub.2, L.sub.4, etc., respectively.

In FIG. 1, the program generation processing 3 executes the following three kinds of processing. A region division processing 4 divides the simulation objective region, based on the shape.division information 1 inputted by the operator, into triangular regions called elements and generates a node.element data table group 7 related to the vertices or nodes thereof and the elements.

The region division is effected for each subregion when the region 41 includes a plurality of subregions as shown in FIG. 3, however, the node.element data table group 7 is formed with respect to the entire region 41.

FIG. 6 is a schematic diagram depicting a result of a division conducted on the subregion S1 in which E.sub.i (i=1, 2, . . . ) and N.sub.j (j=1, 2, . . . ) indicate the element numbers and the node numbers, respectively.

FIG. 7 is a schematic diagram depicting a concrete example of the node.element data table group 7, which will be described in detail later.

Referring to FIG. 1, the equation discretization processing 5 accomplishes processing based on the inputted equation 53 and boundary conditions 54-55 to obtain a function comprising algebraic equations defining coefficients K.sub.ij and constants d.sub.i of simultaneous first-order equations related to values a.sub.i of the unknown quantity A at each node N.sub.i and then stores the function in the discretization table 171.

Based on the boundary condition information included in the equation information 2, the boundary condition table 112 and the boundary side table 113 are generated.

Using the generated node.element data table group 7, boundary condition table 112, boundary side table 113, and discretization table 171, the calculation program generation processing 6 generates a Fortran program code to determine the concrete values of {K.sub.ij } and {d.sub.j } and a Fortran code to call a subroutine program which solves the matrix equation related to {K.sub.ij }, {a.sub.i }, and {d.sub.i } based on the obtained values of {K.sub.ij } and {d.sub.i } (Fortran program 8 of FIGS. 9A-9D).

When a simulation is conducted by the Fortran program thus generated, the calculation program 8 coded in Fortran is translated into a machine language program 10 by a known compiler 9 and the resultant calculation program is loaded on the computer 11. On the other hand, the calculated node.element data table group 7, and boundary side table 113 are inputted as data to the computer 11, thereby obtaining the simulation result 13.

Prior to the detailed description of the program generation processing 3 as a central item of the present invention, a description will be given of the calculation method according to the finite element method which is the basis of the program generation processing 3.

The inputted equation 53 and boundary conditions 54-55 (FIG. 5) will be described before the procedure of the calculation method.

In this equation 53, V indicates a gradient differential operator (.differential./.differential.x, .differential./.differential.y) and V is a vector function (V.sub.x, V.sub.y), where T, V.sub.x, V.sub.y, .beta., and F are functions of space coordinates (x, y) in general. Moreover, "." denotes a scalar product between vectors and "*" designates a product between scalars. Reference numeral 54 is a constant boundary condition indicating that the value of the unknown quantity A must be a constant or fixed on the boundary sides L.sub.1, L.sub.6, etc. of a region .OMEGA.. Reference numeral 55 is an unfixed boundary condition indicating that the expression 55 must hold on the boundary sides L.sub.2, L.sub.4, etc., where n is a unitary vector along the direction of the normal line on the boundary sides of the region and .lambda. and .mu. indicate functions of space coordinates (x, y) in general.

In the finite element method, the calculation is accomplished by use of a basis function defined on each node. FIG. 10 is a schematic diagram illustrating an example of a basis function .phi..sub.i defined on the i-th node N.sub.i. The function .phi..sub.i is 1 on the i-th node N.sub.i and is .0. on a line linking the other adjacent nodes N.sub.j, N.sub.k etc. Namely, the basis function .phi..sub.i is represented in a form of a tent which takes the value of height of a hexagonal cone established by connecting the i-th node N.sub.i to the adjacent nodes N.sub.j, N.sub.k etc.

In the finite element method, the objective unknown quantity A is approximated by a linear combination of basis functions .phi..sub.j (j=1-N.sub.T) as represented by the following expression (1). In place of the unknown quantity A, coefficients {a.sub.j } are derived as the values of the unknown quantity A at the nodes {N.sub.j }. In the expression (1), NODE denotes the number of the nodes in the region .OMEGA. to be simulated. ##EQU1##

To determine the coefficients {a.sub.j }, the expression 53 (FIG. 5) is multiplied by .phi..sub.i and the resultant function is integrated in the simulation region .OMEGA. to be reduced to the following equation (2).

Here, the following autonomous relational expression (3) is used.

where .differential..OMEGA. indicates the boundary of the region .OMEGA.. On the sides L.sub.2, L.sub.4, etc. (to be collectively referred to as a boundary .differential..OMEGA..sub.2 hereinbelow), expression (4) holds as follows from the relationships of the equation 53 of FIG. 5.

By substituting the equations (4) and (1), the following equation (5) results. The equation (5) holds when i indicates a node number of a node on the boundary .differential..OMEGA..sub.2 or a node inside the region .OMEGA.. ##EQU2##

On the other hand, expression (6) holds on the boundary sides L.sub.1, L.sub.6, etc. (to be collectively referred to as a boundary .differential..OMEGA..sub.1 hereinbelow) based on the boundary condition 54 and the equation (1) as follows.

Using the equations (5) and (6), simultaneous first-order equations are obtained with respect to coefficients a.sub.j. The equations can be represented as a matrix equation (7) as follows.

where;

when i denotes a node number of a node on the boundary .differential..OMEGA..sub.2 or of a node inside the region .OMEGA.; and ##EQU3## when i denotes a node number of a node on the boundary .differential..OMEGA..sub.1.

Consequently, solving the simultaneous first-order equations (7), the coefficients {a.sub.j } (j=1, 2, . . . NODE) can be derived. The operation to form the simultaneous first-order equations (7) from the original equation 53 (FIG. 5) is called discretization.

According to the present invention, from the inputted equation 53 and boundary conditions 54-55 (FIG. 5), a Fortran program code is generated to calculate concrete values of {K.sub.ij } and {d.sub.i }, which is characterized in that the code is so generated as to minimize the amount of data required to execute the generated code. Moreover, it is also characteristic that even if the inputted terms of the partial differential equations are not ordered according to the order of differentiation, the Fortran program code can also be generated. The present invention is further characterized in that the Fortran program code is generated according to an inputted specification of an arbitrary region division.

To calculate the values of K.sub.ij and d.sub.i, the values of various integrations contained in the equation (8) must be calculated. Usually, this is accomplished by use of the following various relational expressions.

First, the first term on the right side of the equation (8) is approximated from the following equation (1) or (14). ##EQU4## where, E represents an element including the nodes N.sub.i and N.sub.j and k is a node number of a node other than those of the nodes N.sub.i and N.sub.j and belonging to the element E. T.sub.i, T.sub.j, and T.sub.k denote values of the function T at the nodes N.sub.i, N.sub.j, and N.sub.k respectively. In the following description, the same notation will be applied to the functions .beta., .lambda., and so forth. The integral of the equation (11) can be obtained from the following expression. ##EQU5## holds, where E is an element including the node N.sub.i and k and l represent node numbers other than that of the node N.sub.i in the element. Also in this case, the expressions (12), (13), etc. hold with respect to the integration of the expression (14) if i is substituted by j.

The second term on the right side of the equation (8) is derived from the following expression. ##EQU6## where,

The expressions (16)-(17) can be considered to hold in either cases of i=j and i.noteq.j.

The third term on the right side of the expression (8) can be obtained as follows. ##EQU7##

The fourth term on the right side of the equation (8) is reduced to be as: ##EQU8## where .differential.E represents a side linking the nodes N.sub.i and N.sub.j. Moreover, ##EQU9## where, .vertline.L.sub.ij .vertline. indicates distance between the nodes L.sub.i and L.sub.j.

The first term on the right side of the equation (9) is expressed as: ##EQU10## where,

The second term on the right side of the equation (9) is represented as: ##EQU11## where,

As can be seen from the foregoing description, the integral terms in the expressions (8) or (9) can be represented, based on the expressions (11)-(25), as algebraic expressions by use of the positions of a plurality of related nodes and the function values at the nodes such as T.sub.i, V.sub.xi, V.sub.yi, .lambda..sub.i, f.sub.i, and .mu..sub.i etc.

In this processing, the integral expressions (8)-(9) defining K.sub.ij and d.sub.i of the matrix equation associated with the equation information 2 inputted are first determined, and then algebraic expressions defining K.sub.ij and d.sub.i are derived from the expressions (11)-(25).

In this case, since the terms are not necessarily arranged in the sequence of the order of differentiation, the order of differentiation is determined in this processing, which is characteristic to the present invention. This will be described with reference to FIGS. 11A-11C which is a flowchart according to the Program Analysis Diagram (PAD). The flowcharts subsequent to FIGS. 11A-11C are also generated according to the PAD.

First, the inputted partial differential equation 53 and boundary conditions 54-55 (FIG. 5) are registered to the equation table 111 (FIG. 12) and the boundary condition table 112 (FIG. 13), respectively, and then the boundary side table 113 (FIG. 14) is generated (step 101). The equation table 111 comprises entries 117, 116 and 114 keeping the sign, content, and order of differentiation for each term of the differential equation 53. Separation of the terms can be easily determined by counting the number of addition arithmetic expressions and the number of parentheses contained in the equation. The order of differentiation of each term is stored in a processing step 102 (FIG. 11A) to be described later. The boundary condition table 112 keeps the inputted boundary conditions 54-55 according to the content of each boundary condition, namely, each boundary condition expression and the name of the boundary side on which the condition of the boundary condition expression is satisfied are stored. The boundary side table 113 (FIG. 14) keeps, for each node constituting each boundary side contained in the boundary condition expressions, the node number and the name assigned to the side. The node numbers are obtained from the L table 97 included in the node element data table group 7 (FIG. 7), which will be described in detail later.

Next, in step 102, the equations stored in the equation table 111 (FIG. 12) are checked to identify the order of differentiation for each term thereof. That is, in step 121, the equation table 111 is checked to retrieve an unknown variable A contained in each term. Step 122 then identifies the number and kinds of differential operators related to the unknown variable A. Step 123 sets the number of the differential operators as the order of differentiation to an area 114 of the equation table 111. For the term .beta.*A proportional to the variable A and the term F not including the variable A, the order of differentiation is set to 0 and -1, respectively. That is, if the area 114 is set to a value equal to or greater than .0., the pertinent term contains the variable A; otherwise, the term does not contain the variable A. Step 103 converts each term stored in the equation table 111 into the integral form. That is, the equation 53 (FIG. 5) is converted into the equation (4). First, step 131 multiplies each term by the basis function .phi..sub.i and converts the respective terms of the entry 116 into the integral form 132-135 according to the order of differentiation of each term set to the entry 114. In this case, for a term with the second order of differentiation such as .gradient..(T.gradient.A), the integral form is translated, by use of the partial integral formula (3), into the integral form 132 constituted with the terms having the first order of differentiation. In this procedure, the boundary integral term 136 is derived. In step 137, n.(T.gradient.A) of the boundary integral term 136 is replaced with the item on the right side of the boundary condition equation stored in the boundary condition table 112 (FIG. 13), namely, the translation of the equation into the equation (4) is achieved. This operation is conducted for each different boundary condition and for each boundary edge satisfying the condition.

The translation result is stored in the integral form table 141 (FIG. 15), which contains the sign, content, and the related boundary side name for each term of the integral form equation. Step 139 derives a function comprising algebraic expressions defining K.sub.ij and d.sub.i from the integral terms stored in the integral form equation table 141. That is, step 151 substitutes the unknown variable A of each term in the integral form equation for a linear combination (1) of the basis function and removes the known quantity T from the integral. In step 152, the integral is replaced with the total .SIGMA. with respect to the index j. Steps 151-152 effect a conversion of the equation (4) into the equation (5).

According to the present embodiment, algebraic expressions representing the results of execution of various integral terms such as .intg..gradient..phi..sub.j..gradient..phi..sub.i, .intg..gradient..phi..sub.i, and .intg..phi..sub.i are beforehand stored in the integration evaluation table 161 (FIG. 16). In step 153, as instantiated by the expressions (11) and (14), the known quantity (T) is replaced with an average ##EQU12## of the values at the nodes N.sub.i, N.sub.j, and N.sub.k constituting the element E. Moreover, the number of the basis functions included in the result of the step 152 and the condition indicating whether or not the differential operator .gradient. is applied are identified. Based on the identification result, the integration evaluation table 161 is accessed to retrieve an execution result of the matched integration, and then the integral in the result of the step 152 is replaced with the algebraic expression representing the execution result.

However, for the integrals of which the results vary between i=j and i.noteq.j in the integration evaluation table, namely, for .intg..sub.E .phi..sub.j .phi..sub.i dx and .intg..sub..differential..OMEGA..andgate.E .phi..sub.j .phi..sub.i dL, the values are replaced with the evaluation equations for i=j in the integration evaluation table 161. In any case of the integration, the value associated with i.noteq.j is half the value associated with i=j. If it is desired to utilize the discretization table 171 (FIG. 8) or i.noteq.j, the values of the terms related to .intg..sub.E .phi..sub.j .phi..sub.i dx and .intg..sub.7.OMEGA..andgate.E .phi..sub.j .phi..sub.i dL need only be simply multiplied by 0.5. When the processing described above is effected for all integral terms in the integration equation table 141 (FIG. 15), the discretization table 171 (FIG. 8) is generated. Entry 172 keeps the sign of each term received from the integral form equation table 141, entry 173 stores an algebraic expression as a result of the execution of an integration, and entry 174 keeps an information indicating whether the discretized equation holds inside the region or on the boundary and indicating a boundary side of the boundary where equation holds on the boundary (this can be easily obtained from the integration range of each integral term). Entry 175 remains at 1 if the original integral form equation contains .phi..sub.i, otherwise, the entry 175 remains .0.. Judgement to determine whether or not the equation contains .phi..sub.i can be easily made from the identification result of the step 153. That is, based on the equations (8) and (9), if a term is found to include .phi..sub.i, the term is contained in the expression of k.sub.ij ; otherwise, the term is contained in the expression of d.sub.i. As a result, from the information of the entry 175, it can be determined whether the discretized term contributes to {k.sub.ij } or {d.sub.i }.

Entry 176 is used to indicate, for i.noteq.j, whether or not the equation in the entry 173 is to be multiplied by 0.5 before the equation is further processed.

The node.element table group 7 of FIG. 7 will be described, whereas the region division processing 4 which generates the node.element table group 7 will be described later.

In FIG. 7, the dimension table 910 keeps the number of the nodes NODE, the number of the elements or meshes MESH, the maximum value of the adjacent nodes adjacent to a node BAND, the number of the elements constituting each element 3, and the number of the elements each having two nodes on a side L.sub.i (i=1, 2, . . . ), namely, L1BNUM, L2BNUM, etc.

X table 91 keeps the x coordinate values X.sub.1, X.sub.2, etc. for all nodes, while Y table 92 keeps the y coordinate values Y.sub.1, Y.sub.2, etc. for all nodes. The x and y coordinate values are stored in the sequence of the node numbers.

IADATA table 93 stores, for each element, the node numbers of the nodes contained therein in the ascending sequence of the node numbers.

ICOEF table 94 keeps, for each node, the node numbers of the nodes adjacent thereto in the ascending node number order. If the number of the nodes currently adjacent to a node is less than the maximum number BAND, a number of 0's which number is equal to the difference therebetween is stored for each of the nodes.

L table 97 comprises tables 97-1 to 97-n each keeping, for each side L.sub.i (i=1, 2, . . . ) constituting a region, the node numbers of the nodes contained in each element which have two nodes on the side. For example, for a side L.sub.1 of FIG. 6, the nodes belonging to L.sub.1 includes N.sub.1, N.sub.2, N.sub.6, and N.sub.10, and the elements containing two nodes thereon are E.sub.1, E.sub.4, and E.sub.8. In the table 97-1, the nodes N.sub.1 and N.sub.2 existing on the side L.sub.1 are stored in the first and second columns and the node N.sub.4 existing on other than L.sub.1 is stored in the third column. That is, from the contents of the first and second columns of the table 97-1, the nodes existing on the side L.sub.1 can be determined to include N.sub.1, N.sub.2, N.sub.6, and N.sub.10.

In this processing, a Fortran program code (FIG. 9) is generated from the node.element data table group 7 (FIG. 7) generated by the region division processing 4 and the discretization table 171 (FIG. 8) generated by the equation discretization processing 5 so as to calculate the values of the coefficient matrix {k.sub.ij } and the constant vector {d.sub.i } constituting the matrix equation (7).

The basic concept for generating the Fortran code is as follows. The coefficient matrix {k.sub.ij } and the constant vector {d.sub.i } are represented by the expressions (8) and (9) when the node N.sub.i does not exist on the boundary .differential..OMEGA..sub.1, and the expressions (8) and (9) are converted into algebraic expressions by evaluating the integral terms, and the resultant algebraic expressions are stored in the discretization table 171 in advance. The boundary integration terms .intg..sub..differential..OMEGA..sbsb.2 .lambda..phi..sub.2 .phi..sub.l and .intg..sub..differential..OMEGA..sbsb.2 .mu..phi..sub.i in the expressions (8) and (9) are to be integrated on the boundary .differential..OMEGA..sub.2 ; consequently, if the node N.sub.i does not exist on the boundary .differential..OMEGA..sub.2, the results become .0.. That is, to determine {k.sub.ij } and {d.sub.i } corresponding to a node existing on .differential..OMEGA..sub.1, the calculation need only be effected for the values of the discretized equations (associated with the regional integration terms) which hold inside .OMEGA..sub.2 and which are stored in the discretization table 171 if the node N.sub.i exists inside the region .OMEGA..sub.2. If the node N.sub.i exists on .differential..OMEGA..sub.2, in addition to the discretized equations which hold inside the region and which are stored in the discretization table 171, the discretized equations (associated with the boundary integration terms) which hold on the boundary side region including the N.sub.i need only be calculated.

Whether a discretized term is used in the calculation of {k.sub.ij } or {d.sub.i } can be determined by referring to the information kept in the entry 175 as described above.

On the other hand, when the N.sub.i exists on .differential..OMEGA..sub.1, the values of {k.sub.ij } and {d.sub.i } are determined from the equation (10). In this case therefore, the discretization table need not be used, and the determination of the values of {k.sub.ij } and {d.sub.i } can be effected independent of the node not existing on .differential..OMEGA..sub.1.

A description will be given of a case where the discretized equations which hold inside the region .OMEGA..sub.2 are evaluated.

To calculate the values of the discretized equations which hold inside the region is identical to achieving the regional integration equivalent thereto; consequently, the calculation of the values of the discretized equations can be conducted by calculating the value of the discretized equation for each element contained in the integration region and by adding the resultant values.

FIG. 18A shows a region in which the values of the discretized equations are calculated. To obtain the values of k.sub.ii and d.sub.i, the values of the discretized equations are calculated on the elements E.sub.s -E.sub.x including the nodes N.sub.i and the obtained values are added to each other; whereas to calculate the values of k.sub.ij (i.noteq.j), the values of the discretized equations are calculated on two elements E.sub.w and E.sub.x including the nodes N.sub.i and N.sub.j and the resultant values are added to each other.

Element E (comprising the nodes N.sub.i, N.sub.j, and N.sub.k) which is included in the integration region for obtaining the values of nine matrix elements k.sub.ii, k.sub.ij, k.sub.ik, k.sub.ji, k.sub.jj, k.sub.jk, k.sub.ki, k.sub.kj, and k.sub.kk and three constant vector components d.sub.i, d.sub.j, and d.sub.k as shown in FIG. 18B is selected; then, the calculation of the discretized equations associated with the contribution from the element E to k.sub.ii -k.sub.jj and d.sub.i -d.sub.k is effected at a time. Processing for all elements, can be considered to be identical to calculating the values of the discretized equations for obtaining the values of all coefficient matrices {k.sub.ij } and all constant vectors {d.sub.i }.

The calculation of the values of the discretized equations which hold on the boundary side region can also be accomplished in the same fashion.

Based on the method described above, the Fortran code for achieving the procedure to determine {k.sub.ij } and {d.sub.i } can be generated. In this processing, the calculation can be conducted for each element by use of the element table 93 indicating the nodes included in the element and the x and y coordinate tables 91-92.

Conventionally, the final values of k.sub.ii, k.sub.il -k.sub.iNODE, and d.sub.i are calculated for each node by use of the data indicating, for each node other nodes which are connected to the node and a number of the element to which the above other nodes belong, and hence the amount of the data is considerably increased, however, the data associated with this processing is unnecessitated by the present invention.

The following paragraphs will describe the flowchart of FIGS. 17A-17C. Step 182 generates a program code in which the entry 174 (see FIG. 8) is searched in the discretization table 171 to find out a plurality of discretized equations related to the regional integration, and then the value of each discretized equation thus obtained is calculated for all nodes of all elements, thereby calculating a portion of the values of {k.sub.ij } or {d.sub.i }, the portion being determined depending on the regional integration terms.

Consequently, when the program codes generated by the processing steps 181 and 182 are executed, the final values of the matrix elements {k.sub.ij } and the constant terms {d.sub.i } can be derived for the nodes inside the region to be simulated and the nodes on the boundary .differential..OMEGA..sub.2 with variable values (sides L.sub.2, L.sub.4, etc. in the case of FIG. 5).

Step 182 is a section where the matrix elements {k.sub.ij } and the content terms {d.sub.i } determined according to the node numbers of the nodes on the boundary .differential..OMEGA..sub.1 with fixed values (sides L.sub.1, L.sub.6, etc. in the case of FIG. 5) are obtained from the equation (10).

In the steps 181 and 182 of the embodiment, based on the positions of all nodes including the nodes on the boundary .differential..OMEGA..sub.1, the values of the discretized equations related to the regional integration or the values of the discretized equations related to the boundary integration are first calculated, and then the matrix elements k.sub.ij and d.sub.i related to the nodes on the boundary .differential..OMEGA..sub.1 with fixed values are replaced with the values obtained by the step 182. Finally, step 184 generates a code for calling a matrix solution Fortran code routine (not shown) beforehand provided which derives {a.sub.j } from the equation (1) based on the {k.sub.ij } and {d.sub.i } thus determined. A description will be given of the details about the respective processing in conjunction with the example of the code generated as shown in FIGS. 9A-9D.

The program to be generated here first reserves the RCOEF table 190 and CONS table 191 (FIG. 19) to store the calculated values of {k.sub.ij } and {d.sub.i } in the memory (not shown) of the computer (not shown) and then conducts the calculation of the values of {k.sub.ij } and {d.sub.i }.

In FIGS. 9A-9D, RCOEF (p, q) and CONS (l) are codes for referencing the entries for k.sub.pq and d.sub.l in the RCOEF table 190 and the CONS table 191, respectively. Similarly, X(p), Y(p), IADATA (p, q), and ICOEF (p, q) are codes for referencing the X table 91, the Y table 92, IADATA 93, and the ICOEF table 94 of FIG. 9.

According to this processing method, an element E (comprising the nodes N.sub.i, N.sub.j, and N.sub.k) is selected so as to effect the calculation at a time for determining the values of the discretized equations associated with nine coefficient matrix components k.sub.ii, k.sub.jk, k.sub.ik, k.sub.ji, k.sub.jj, k.sub.jk, k.sub.ki, k.sub.kj and k.sub.kk and three constant vector components d.sub.i, d.sub.j, and d.sub.k on the element E and the values are added to the respective coefficient matrix elements and the constant vector components. However, if all such codes generating the discretized equations of the 12 components are generated, the program becomes redundant. To overcome this difficulty, p, q, and r are introduced to assign i, j, and k in the different fashions, respectively; consequently, only the calculation codes associated with k.sub.pp, k.sub.pq, and d.sub.p are generated, and the calculation is conducted for the 12 components by sequentially changing i, j, and k to be assigned to p, q, and r. Namely, p indicates one of i, j, and k; q indicates one of i, j, and k and is different from the item assigned to p; and r indicates one of i, j, and k and is different from those items assigned to p and q. Assignment of p, q, and r and a concrete example thereof will be described in the following explanation about the code generation processing hereinbelow.

First of all, step 191 (FIG. 17A) generates a code 211 (FIG. 9A) for initializing the RCOEF table 190 and the CONS table 191. Steps 192-200 generate a code in which for each element the discretized equations corresponding to the regional integration of the coefficient matrix components k.sub.ii, k.sub.ij, k.sub.ik, k.sub.ji, k.sub.jj, k.sub.jk, k.sub.ki, k.sub.kj, and k.sub.kk and the constant vector components d.sub.i, d.sub.j, and d.sub.k of the element (to be represented as E comprising the nodes N.sub.i, N.sub.j, and N.sub.k) are calculated, and then the obtained values are added to the respective areas allocated to the respective components in the RCOEF table 190 and the CONS table 191 as shown in FIG. 19. Step 192 generates a code 212 as a DO loop which repeats the execution of the code generated by the steps 193-200 as many times as there are elements. Step 193 generates a code 213 as a DO loop which repeatedly executes the code generated by the steps 194-200 as many times as there are nodes of the element (i.e., three times). Step 194 generates a code 214 which obtains the nodes constituting each element from the IADATA table 93 and assigns the node numbers thereof to the p, q, and r described above. If the component nodes are N.sub.i, N.sub.j, and N.sub.k and the node numbers thereof are i, j, and k, then i, j, and k are assigned to p, q, and r, for example. The DO loop generated by the step 194 has a function which changes the node number to be assigned to p among the p, q, and r as i, j, and k. Step 195 generates a code which calculates, from the expression (13), b.sub.p, c.sub.p, and .DELTA.E based on the information stored in the x coordinate table 91 and the y coordinate table 92. Here, b.sub.p, c.sub.p, and .DELTA.E are obtained by replacing i with p, j with q, and k with r. Step 196 generates a code in which, based on the added result of the equation information of the entry 173 associated with the terms 1711, 1714, and 1715 for which the entry 174 of the discretization table 171 contains .OMEGA. (indicating the discretized equations which hold inside the region .OMEGA.) and for which the entry 175 contains 1 (indicating the discretized terms contributive to the matrix elements {k.sub.ij }), the discretized equations associated with the regional integration terms contributive to the matrix elements {k.sub.ij } in each element are calculated and the resultant values are added to the entry k.sub.pp of the RCOEF table 190. The step 196 further generates a code in which, based on the equation information of the entry 173 associated with the term 1716 for which the entry 174 of the discretization table 171 contains .OMEGA. and for which the entry 175 contains .0. (indicating the discretized equations contributive to the constant vectors {d.sub.i }), the discritized equations corresponding to the regional integration terms contributive to the constant vectors {d.sub.i } in each element are calculated and the resultant values are added to d.sub.p of the CONS table 191 (code 216). Since the node number assigned to p by the DO loop generated by the step 194 changes as i, j, and k in this order, the values of the discretized equations associated with the regional integration of k.sub.ii, k.sub.jj, k.sub.kk d.sub.i, d.sub.j, and d.sub.k in the element E are calculated by the step 196, and as a result, the values are added to the respective areas of the components as shown in FIG. 19. Step 197 generates a code 217 as a DO loop which repeals an operation as many times as indicated by <number of the nodes of the element>-1 (namely, two times). Step 198 generates a code 218 which, without changing the node number assigned to p, assigns the different node numbers to q and r. For example, if i has been assigned to p, j and k are assigned to q and r, respectively, or k and j are assigned thereto. The DO loop generated by the step 197 is used to change the node numbers assigned to q and r among the p, q, and r.

Step 199 generates, based on the equation (13), b.sub.p, c.sub.p, b.sub.q, c.sub.q, and .DELTA.E from the information stored in the X table 91 and the Y table 92. Like step 196, step 200 generates a code 220 in which, based on the discretization table 171, the calculation of the discretized equations associated with the regional integration contributive to the matrix elements {k.sub.ij } of each element is conducted and then the resultant values are added to the entry k.sub.pq in the RCOEF table 190. Since the DO loop generated by the step 193 changes the node numbers assigned to q and r; for example, if i has been assigned to p, the discretized equations associated with the regional integration of k.sub.ij and k.sub.ik of the element E are calculated and the obtained values are added to the respective entries of k.sub.ij and k.sub.ik shown in FIG. 19. If j has been assigned to p, the discretized equations related to k.sub.ji and k.sub.jk are calculated; whereas, if k has been assigned to p, the discretized equations associated with k.sub.ki and k.sub.kj are calculated.

For each boundary side having a boundary condition with a variable value, a code (not shown) is generated to execute the processing of the steps 201-207. Step 201 generates a code 221 as a DO loop for repeating an operation as many times as there exist elements each having two nodes on each side. Step 203 generates a code 223 as a DO loop which repeats an operation as many times as there exist nodes of the pertinent element (i.e., two times). Step 204 generates a code 224 which obtains the nodes of each element on the boundary by use of the boundary side table 97 and assigns the node numbers of the obtained nodes to p and q. If the nodes on the boundary are N.sub.i and N.sub.j and the node numbers thereof are i and j; i and j or j and i are assigned to p and q, respectively. The DO loop generated by the step 203 is used to change the node numbers assigned to p and q. Step 205 generates a code which calculates length of a line segment L.sub.pq between two nodes having the node numbers p and q. Step 206 generates a code in which, if the pertinent boundary side is L.sub.2, based on the equation information of the entry 173 related to the item 1712 for which the entry 174 is L.sub.2 and the entry 175 is 1 in the discretization table 171, the calculation of the discretized equations associated with the boundary integration terms contributive to the matrix elements {k.sub.ij } of a line segment (to be referred to as a boundary line segment herebelow) comprising two nodes on the boundary related to the element is effected and then the obtained values are added to the entry allocated to k.sub.pp in the RCOEF table 190; the step 206 further generates a code in which, based on the equation information of the entry 173 related to the item 1713 for which the entries 174 and 175 include L.sub.2 and .0., respectively in the discretization table 171, the calculation of the discretized equations associated with the boundary integration terms contributive to the constant vector d.sub.i in the boundary line segment is conducted and then the obtained values are added to the entry d.sub.p in the CONS table 191 (code 226). Like step 206, step 207 generates a code 227 in which, based on the discretization table 171, the calculation of the discretized equations associated with the boundary integration terms contributive to the matrix elements {k.sub.ij } of the boundary line segment is achieved and the obtained values are added to the entry of k.sub.pq in the RCOEF table 190. The steps 201-207 generate the codes only for the boundary regions for which the boundary conditions with variable values are defined as described above. Since the information indicating the boundary conditions and the regions associated therewith is kept in the field 174 of the discretization table 171, the codes are so generated as to reference the information.

Steps 208, 209, and 2010 in step 183 generate codes 228, 229, and 2210, respectively which set the boundary conditions represented by the equation (10) to the nodes on the boundary .differential..OMEGA..sub.1 with fixed values.

To determine whether a side is on a fixed-value boundary or not, the information stored in the boundary condition table 112 (FIG. 13) need only be referenced, and the number of the nodes and the node numbers on the boundary are kept in the side table 97 (FIG. 7). Finally, the step 184 generates a code 2211 in which, to obtain {a.sub.j } based on the RCOEF table 190 and CONS table 191, a matrix solution routine prepared by the system is called. The values of {a.sub.j } are thus obtained, which completes the generation of the codes.

Among the calculation program generation flow, the steps 196, 1910, 206, and 207 excepting the integral calculation code generation are constituted with codes of the fixed structure, namely, these codes are independent of the types of expressions and boundary conditions. Consequently, the code generation does not require any special idea. The integral calculation codes are so generated as to replace the equation information kept in the discretization table 171 (FIG. 8) with the corresponding Fortran designations and symbols.

The operation flow of the program generation processing 3 has been described.

Prior to this processing, the shape division information 1 will be described.

FIG. 5 shows a concrete example of the shape division information 1.

This is an example of the shape.division information 1 in a case where a region 41 comprising a subregion S1 defined by points P.sub.1, P.sub.2, P.sub.3, P.sub.4, P.sub.5 and a subregion S2 defined by points P.sub.5, P.sub.4, P.sub.3, P.sub.6, P.sub.9, P.sub.8, and P.sub.7 in FIG. 3 is set as an object of the simulation. The information includes a group of program statements defining the shapes of the areas and a group of program statements dividing the areas into elements. In the first program statement L.sub.1 =LINE (P.sub.1, P.sub.2), R(3, r.sub.1), LINE (P.sub.1, P.sub.2) defines a line segment L.sub.1 having two end points P.sub.1 and P.sub.2, whereas R (3, r.sub.1) indicates to divide the line segment into three partitions with a common ratio r.sub.1.L.sub.2 =LINE (P.sub.2, P.sub.3), D(2) means that a line segment L.sub.2 having two end points P.sub.2 and P.sub.3 is to be divided into two equal parts. ARC in L.sub.3 indicates an arc to be defined, whereas SPLINE in L.sub.6 denotes a spline curve. L.sub.1 to L.sub.7 are related to the boundary sides necessary to define the areas S1 and S2, which are to be automatically divided into elements. For example, a statement S1=QUAD (P.sub.1, P.sub.2, P.sub.3, P.sub.5), A indicates that an automatic region division is to be conducted on a region enclosed with the sides L.sub.1 -L.sub.4 defined by the points P.sub.1, P.sub.2, P.sub.3, and P.sub.5 in advance.

Next, the content of the region division processing 4 will be described with reference to FIG. 20.

First, division of the boundary sides L.sub.1 -L.sub.7 of the subregions S1 and S2 (FIG. 3) is accomplished according to the shape.division information 1 (FIG. 4; step 61). Namely, the line segments of the boundary sides are divided into equal parts or according to the equal ratio, thereby determining the nodes on the sides.

Subsequently, the inside areas of the subregions S1 and S2 are divided into elements according to the division of the boundary sides (step 62).

The processing flow of the automatic division will be described with reference to FIG. 21. In this diagram, reference numeral 71 schematically denotes a subregion obtained after the sides L.sub.1 -L.sub.4 are divided, and a small circle represents a node on the side thus determined by the division of the side in the step 60. Here, the number of divided parts is equal for the opposing sides L.sub.2 and L.sub.4, but the number of divided parts varies between the sides L.sub.1 and L.sub.3. In this case, virtual nodes are located between nodes on the side L.sub.1 having a smaller number of divided parts, where the number of virtual nodes is indicated by the difference between the number of divided parts of the side L.sub.3 having a larger number thereof and that of the side L.sub.3 having a smaller number thereof (step 67). In this case, since the virtual node cannot be a vertex of the element, as the number of virtual nodes increases, the density of elements after the region division lowers, consequently, it is essential that the virtual nodes should not be concentrated on a single location. In this processing, therefore, the virtual nodes are so located between the real nodes as to be dispersed therebetween. The virtual nodes to be located between the real nodes are selected according to a predetermined function, for example, a function linear with respect to the distance between the virtual nodes and an end of the side having the smaller number of divided parts. As a result, the number of real nodes becomes equal to the number of virtual nodes for a pair of the opposing sides. If a subregion is a triangular region, a side having the greatest number of real nodes is selected from three sides of the triangle. Assuming the side to be divided into two sides by an appropriate node thereon, namely, the triangle is regarded as a rectangular region, thereby applying completely the same processing thereto. Next, consider a normal coordinate region 72 having the same division number as that of the boundary sides of the subregion S1. The region 72 includes an orthogonal lattice in which the distance between the lattice points is fixed to 1. The system generates a function F which maps each lattice point on each boundary of the region 72 onto one of the real or virtual nodes on the boundary of the region 71. For example, a lattice point 77 is mapped onto a node 78 (a virtual node in this case) on the side L.sub.1 of the subregion S1. Several methods can be considered to generate the mapping function. For example, as is well known, the mapping from the region 72 onto the region 71 can be defined under conditions that transformation between the coordinates satisfies the Laplace equation and that the lattice points on the boundary of the lattice region 72 are mapped onto the real or virtual nodes on the sides of the region 71. Next, among the lattice points of the region 72, the lattice points (virtual points) corresponding to the virtual nodes of the region 71 are obtained. Such virtual lattice points are indicated by X on the boundary sides of the region 72. Furthermore, among the lattice points inside the region 77, the lattice points to be regarded as the virtual lattice points are determined (step 68). The determination of the virtual lattice points inside the region is stepwise achieved such that for each determination, two virtual lattice points on the ends are removed from the group of virtual lattice points adjacent to each other on the boundary. For example, the virtual lattice points 88 and 90 on the ends are removed from the virtual lattice points 88, 89, and 90, and the lattice point 91 existing at an upper location is assumed to be a virtual lattice point. From the virtual lattice points 72 and 92, any virtual lattice point cannot be generated. Next, the region 72 is divided into small rectangular regions each comprising real lattice points as the four corner points of the rectangle. Refernce numeral 79 indicates a diagram showing the result of division. Finally, in each small rectanguler region, the real lattice points thereof are appropriately connected to each other to generate a group of rectangular elements (step 69), which are indicated by reference numeral 75. The method for connecting the lattice points need only be determined so as to minimize the length of the sides resulting from the connection.

Next, according to the mapping function F generated beforehand, the region 5 is mapped onto the region 71 to obtain the element division diagram 76. The inside of the subregion S1 has been thus divided. The other subregions are also processed in the same fashion. Thereafter, according to the region division results, the node.element data table group (FIG. 7) is generated (step 63).

An embodiment has been described, but various changes or modifications thereof can be easily considered, for example, the shape.discretization information 1 and the mesh division processing 4 may be replaced with a CAD system, moreover, the node.element data table group 7 need not be necessarily inputted at an execution of the calculation but the table group 7 may be incorporated as a code in the calculation program 10.

Although the embodiment has been described in a case of a two-dimensional space region, the automatic generation of the programs can also be implemented for a one- or three-dimensional space region in which only the integral formula of the basis function is changed.

Since the calculation code is determined according to the equations, this method can be applied to various equations.

In addition to the kinds of basis functions of FIG. 10, basis functions of the higher order can also be easily introduced to the system by changing the integral formula and the formula for transforming equations.

The program language of the generated program is not limited to FORTRAN, for example, PL/I and PASCAL having the same level of description function as FORTRAN are applicable to the automatic program generation.

According to the present invention, the generated program code requires a reduced memory capacity for the program code execution as compared with the prior art technique.

The number of subregions to be obtained by the region division can be arbitrarily controlled, which enables one to generate a program accomplishing a numerical simulation with desired precision.

Moreover, the system enables generation of a program achieving a numerical simulation on arbitrary mathematical model without imposing any restrictions on the objective partial differential equations and boundary conditions.