Computer aided design system

A computer-aided design system in which a mechanical quantity of a structure can be determined from stored geometrical representation of the structure, the material properties thereof and the loads imposed thereon, comprising: generating data describing a structure as a mesh; generating by application of the finite element method to the mesh, the elements of a stiffness matrix, a loading vector, and a vector including an associated degree of freedom index for each row and column of the stiffness matrix, generating from the stiffness matrix and the loading vector elements of matrix A and a right hand-side vector f, where A and f are related through the relation Ax=f, the vector x representing a mechanical quantity at a plurality of points of the mesh; and a linear solver of the preconditioned conjugate gradient type for generating, for subsequent processing or display, the elements of the vector x from A and f using a preconditioning matrix, the system being characterised in that the linear solver comprises generating a principal submatrix from the stiffness matrix by setting to zero those elements therein in positions such that the degree of freedom index associated with their row is not equal to the degree of freedom index associated with their column and generating the preconditioning matrix from the principal submatrix.

SUMMARY OF THE INVENTION 
Disclosed herein is an extension of the ordinary preconditioned conjugate 
gradient method which enables solutions for a plurality of loadings and 
which avoids the necessity of "starting from scratch" each time. The 
invention, in one aspect, is a technique for converting a stored 
geometrical representation of the structure, the material properties 
thereof and the loads imposed thereon into a visualization of a mechanical 
quantity of the structure. Initially, a mesh is generated from the stored 
geometrical representation of the structure. The technique then includes 
generating, by application of the finite element method to the mesh, the 
elements of the stiffness matrix and a loading vector; generating from the 
stiffness matrix and the loading vector a matrix A and a right hand-side 
vector f, where A and f are related through the relation Ax=f, the vector 
x representing the mechanical quantity at a plurality of points of the 
mesh; generating a preconditioner to accelerate the solution to Ax=f; 
approximating solutions of a plurality of problems Ax.sub.k =f.sub.k, the 
plurality of problems having a common matrix A and a plurality of 
different loadings/right hand-sides f.sub.k, including recursively 
constructing for each k an A-orthogonal basis W(k) using W(1), W(2) . . . 
W(k-1), and obtaining in approximation to x.sub.k using W(k); and 
providing a user the solutions x.sub.k. 
In another aspect, the invention relates to a technique for generating a 
numerically stable A-projection y.sub.q of x upon an A-orthogonal basis 
sequence W, which includes: partitioning the A-orthogonal basis sequence W 
into a plurality of subsequences; calculating an A-projection of x upon a 
given subsequence of the plurality of subsequences using a residual value 
r.sub.i-1, including calculating a subsequent residual vector r.sub.i, for 
i=1; accumulating the A-projection into a cumulative A-projection y.sub.i 
of x, for i=1; repeating the calculation of the A-projection for i=2 . . . 
q upon respective subsequent subsequences of the plurality of 
subsequences; and repeating the accumulating such that the numerically 
stable cumulative a-projection y.sub.q of x upon W results having 
associated therewith a numerically stable residual vector r.sub.q. 
By using the above techniques, a fast solution is provided when switching 
from one loading to another, while the stiffness matrix remains unaltered. 
This procedure eliminates the need to "start from scratch" each time the 
loading is replaced by a new loading, which is typical of many of the 
iterative solution techniques currently available. The disclosed technique 
therefore considerably shortens the processing time required for the 
entire design phase.

DETAILED DESCRIPTION OF THE INVENTION 
FIG. 1 is a schematic diagram of a computer aided design system according 
to the invention. The geometry of a structure, finite element mesh 
information, elements information, material properties, Boundary 
conditions and loading information are stored in storage device 10 under 
the control of Data Base Management System (DBMS) 20. This data is input 
via a front-end graphical interface 30. A display system 40, which 
includes a visual display unit and a printer, is used for displaying 
images of the structure and the stress & strain tensor, as well as the 
displacements (deformation), at prescribed portions of the structure using 
appropriate data visualization techniques. The display system also 
includes software for relevant report generation for display or printing. 
The information stored in the database is used to control numerical control 
module 50 which generates numerical control code for controlling machine 
tools used in manufacturing the structure being designed. Tool design 
module 60 handles the design of the tools required to make the structure. 
Inspection and testing module 70 generates data for use in inspection and 
testing of the manufactured structure. The general operation of these 
modules is well understood in the art and will not be further described 
herein. 
The CAD system also includes a finite element method (FEM) module 80 which 
enables a stored structure to be tested and various of its properties to 
be determined using finite element analysis. 
FIG. 2 is a schematic diagram showing the Finite element analysis module 80 
and the interface with DBMS 20. FE module 20 comprises a graphical mesh 
generator 90 to generate the division of the structure into elements, 
equation generator 100 for generating the system of algebraic equations 
which describe the behavior of the structure, and solver 110 for solving 
the generated equations to calculate the desired property. Discretization 
of the problem by finite element method (FEM) leads to a linear system of 
equations, generated by the equation generator, Ax=f, where A=(A.sub.ij) 
is a stiffness matrix, f=(f.sub.i)--right hand side (r.h.s) vector and 
x=(x.sub.i)--vector of unknowns i, j=1,2, . . . ,n. 
In structural analysis, it is generally necessary to solve the following 
types of elastic problems, where the displacements of parts of the 
structure are unknown parameters: 
a) Linear Equilibrium problems; 
b) Non-Linear Equilibrium problems; 
c) Dynamical problems which include forces of inertia; and 
d) Eigenvalue/eigenvector problems--Normal mode dynamics. 
Various theories of elasticity can be used to formulate these problems, all 
of them resulting in a system of elliptic partial differential equations. 
For dynamical problems, the space derivative part is elliptic. Various 
types of finite elements can be used, including Serendipic and P-elements. 
When the finite element method is applied to the relevant system of partial 
differential equations which corresponds to a linear equilibrium problem, 
one obtains, a system of linear algebraic equations, in which a stiffness 
matrix multiplies the vector of unknown displacements, while a certain 
right-hand side vector prescribes the loadings information. 
In non-linear equilibrium problems as well as in various dynamical 
problems, the FEM discretization leads to a system of algebraic equations, 
and the application of various iterative techniques devised for their 
solution results, again, within each iteration, in a similar system of 
linear algebraic equations. 
There exist a maximum of 6 degrees of freedom at each point of an elastic 
body: 
u.sub.1, u.sub.2, u.sub.3 --translations; and 
u.sub.4, u.sub.5, u.sub.6 --rotations. 
Each entry x.sub.i of the solution vector x corresponds to a node j and a 
degree of freedom m .epsilon. (1,2, . . . , 6). A vector L of degrees of 
freedom is generated which has, as i-th component, that value of m which 
corresponds to the component x.sub.i. 
Solver module 110 is shown schematically in FIG. 3 and comprises two 
parts--an iteration mechanism 120, for handling non-linear and dynamic 
problems, and a linear solver 130 for solving both linear equilibrium 
problems as well as the linear systems arising within each iteration in 
solving non-linear and dynamic problems. 
Linear solver 130 is the critical element in the operation of the whole CAD 
system. It must be capable of handling Single & Multi-Point Constraints, 
be efficient and sufficiently accurate even handling large condition 
numbers and maintain its efficiency in solving for cases with variable 
loadings, or arbitrarily specified multi-loading systems, by avoiding, as 
far as possible, the need to start from scratch whenever a new 
loading-system is being treated. 
FIG. 4 is a schematic block diagram illustrating the linear solver. It 
takes as input the stiffness matrix and the right hand side vector 
generated by equation generator 100. From these a preconditioning matrix 
is generated and, using this preconditioning matrix, the preconditioned 
conjugate gradients method is applied using a suitable computer subroutine 
to solve the equations generated by the equation generator. The 
preconditioned conjugate gradients method was originated by R. Hestenes, 
E. Steifel, `Methods of Conjugate Gradients for Solving Linear Equations,` 
Journal of Research of the National Bureau of Standards, 49, 409-435 
(1952). The present form of PCG is described in the book by Gene H. Golub, 
Charles F. Van Loan, `Matrix Computations,` John Hopkins University Press 
(1984). 
In order to solve the equation Ax=f, the PCG method with the preconditioner 
K proceeds as follows. A trial solution x.sub.0 is selected. Vectors 
r.sub.0 and P.sub.0 are defined as r.sub.0 =f-Ax.sub.0 and P.sub.0 
=K.sup.-1 r.sub.0. 
The following steps are then performed iteratively until convergence is 
achieved: 
EQU .alpha..sub.i =(K.sup.-1 r.sub.i-1, r.sub.i-1)/(p.sub.i-1, Ap.sub.i-1)(1) 
EQU x.sub.i =x.sub.i-1 +.alpha..sub.i p.sub.i-1 (2) 
EQU r.sub.i =r.sub.i-1 -.alpha..sub.i Ap.sub.i-1 (3) 
EQU .beta..sub.i =(K.sup.-1 r.sub.i, r.sub.i)/(K.sup.-1 r.sub.i-1, r.sub.i-1)(4 
) 
EQU p.sub.i =K.sup.-1 r.sub.i +.beta..sub.i p.sub.i-1 (5) 
where the expression in the form of (Vector.sub.-- x, Vector.sub.-- y) 
represents the scalar product of the vectors Vector.sub.-- x and 
Vector.sub.-- y. 
A suitable efficient matrix-vector multiplication, which can be adapted for 
the particular computer being used, is used for the multiplication 
operations. 
FIG. 5 shows the steps performed in the generation of the preconditioning 
matrix. A `principal` submatrix B=(b.sub.ij) is extracted from the 
stiffness matrix A=(a.sub.ij) by setting to zero those elements in the 
stiffness matrix for which the degree of freedom index associated with 
their row is not equal to the degree of freedom index associated with 
their column, using a method described in the following pseudo-code: 
do i=1, n 
do j=1, n 
if l.sub.j =l.sub.i then b.sub.ij =a.sub.ij 
else b.sub.ij =0 
enddo 
enddo 
Where l.sub.i are components of vector L. The resulting matrix is a matrix 
in which the coupling between certain degrees of freedom is eliminated. 
Following both a theoretical and experimental verification, the inventors 
have found that by using this process either the matrix B defined above or 
the matrix B' can be used as preconditioning matrices for the stiffness 
matrix, where B' is the result of applying the incomplete LU factorization 
to the matrix B. 
The incomplete L-U factorization is obtained as follows. If it is desired 
to solve the matrix equation Ax=f, where A is a large sparse nxn matrix, 
choosing a Gaussian elimination would lead us to a factorization A=LU. 
However this method has the problem that the number of non-zero elements 
is much greater than the number of non-zero elements in the sparse matrix 
A, and since only the non-zero elements are stored this results in a high 
cost in terms of storage and complexity. 
To avoid these difficulties an incomplete factorization technique is used 
to obtain a splitting A=LU--R, where the fill-in has been controlled by 
zeroing previously determined coefficients in L and U during the 
elimination. Such a technique is described in D. S. Kershaw, `The 
Incomplete Cholesky-Conjugate Gradient Method for the Iterative Solution 
of Systems of Linear Equations,` Journal of Computational Physics 26, 
43-65 (1978). 
In the case of a boundary value problem for an elastic body with rigid 
inclusions, the problem is usually formulated as a boundary value problem 
for an elastic body with holes instead of rigid inclusions, but with 
multipoint constraints generated by the rigid inclusions. It is known that 
a rigid element has only 6 independent degrees of freedom in 3D problems 
(3 in 2D). For example it may consist of the degrees of freedom of a 
single node. 
Other degrees of freedom for 3D rigid elements can be found by the 
kinematic relations, 
EQU u.sub.j =u.sub.0 r.sub.j0 x.omega..sub.0 ; .omega..sub.j =.omega..sub.0, 
where u(u1,u2,u3) is the vector of displacements and .omega. (u4,u5,u6) is 
the vector of rotations, j is any node with dependent degrees of freedom 
and r.sub.j0 is a pointer vector from point 0 to point j. 
If the total number of degrees of freedom g consists of n independent 
degrees of freedom and m dependent degrees of freedom, then the relation 
between the dependent degrees of freedom and the independent ones can be 
written in matrix form as u.sub.m =G.sub.mn u.sub.n, where G.sub.mn is an 
m.times.n rectangular matrix. 
The saddle point formulation and finite element method are used for 
discretization of the problem, giving rise to the following algebraic 
equation, 
##EQU1## 
where K.sub.gg denotes the FEM matrix for the elastic body with holes 
replacing the rigid inclusions. This will be referred to as the `hole 
stiffness matrix.` G.sub.mg is the matrix of multipoint constraints, 
P.sub.g is the loading vector and .lambda. the vector of Lagrange 
multipliers. I.sub.mm is the m.times.m unit matrix, and G.sup.t.sub.mg 
denotes the transpose of G.sub.mg. 
In this case a principal submatrix K.sup.d.sub.gg is extracted from the 
hole stiffness matrix K.sub.gg using the method described above. It has 
been found by the inventors that a reduced resolving matrix 
EQU R.sub.nn =G.sup.t.sub.gn K.sup.d.sub.gg G.sub.gn 
where 
##EQU2## 
or the result of applying the incomplete factorization to R.sub.nn as 
above can be used as a preconditioner for the resolving matrix K.sub.nn in 
the equation, 
EQU K.sub.nn u.sub.n =G.sup.t.sub.gn P.sub.g, K.sub.nn =G.sup.t.sub.gn K.sub.gg 
G.sub.gn, u.sub.g =G.sub.gn u.sub.n 
This equation can then be solved in the same way using the preconditioned 
conjugate gradients method. 
In practice, many applications of structural analysis require a system of 
algebraic equations to be solved for several loading vectors/right-hand 
sides. For example, this situation arises in linear problems with several 
loadings, in the inverse power method for normal mode dynamics and in 
iterative methods for non-linear problems. 
One way of achieving this is by treating each system separately, for each 
i, i=1 . . . k, via PCG, using the above preconditioner. However, this 
would be expensive in terms of computer time. 
These situations are handled in this embodiment of the invention in the 
following manner, which is both efficient & numerically reliable. 
If the system to be solved is Ax.sub.k =f.sub.k, for k=1, 2, 3 . . . N, the 
linear problem for the first right hand side vector f.sub.1 is solved 
using the method described above. The sequence of coordinate vectors to 
the base sequence W=P(1)={p.sup.1.sub.1, p.sup.1.sub.2 . . . , 
p.sup.1.sub.m1 } obtained during the PCG procedure described above are 
stored. For the subsequent right hand side vectors f.sub.k the solutions 
x.sub.k are found in two steps. In the first step the A--orthogonal 
projection y.sub.k of x.sub.k on the space generated by the basic sequence 
W, i.e. the span W, as well as the corresponding residual r.sub.k are 
calculated. In the second step, the linear problem Az.sub.k =r.sub.k is 
solved using the PCGM algorithm described below. The sequence of 
coordinate vectors P(k)={p.sub.1.sup.k, p.sub.2.sup.k, . . . , 
p.sub.mk.sup.k } is now added to the basic sequence W to obtain an updated 
extended basic sequence. The approximate solution for x.sub.k is given by 
x.sub.k =y.sub.k +z.sub.k. 
A straightforward calculation of the A-orthogonal projection y.sub.k of 
x.sub.k on the current span W, i.e., 
##EQU3## 
suffers from loss of precision, since small inner products of relatively 
large vectors are hard to resolve accurately. Thus in the preferred 
embodiment a cycling scheme is used for computing the projection y of x on 
the W, as follows. 
A parameter s is chosen between 1 and m. A quantity q is defined as 
q=m/s!. Variables y, y and y are initialized to zero. r.sub.0 is set to 
f.sub.k. For i=1, 2, . . . , q, the following steps are performed. 
EQU j.sub.1 =1+(i-1)s, j.sub.2 =is (1) 
EQU r.sub.i =r.sub.i-1 -y.sub.i-1 (3) 
##EQU4## 
EQU y.sub.i =y.sub.i-1 +y.sub.i (5) 
Then set y=Y.sub.q and r=r.sub.q. 
If the parameter s is increased, the accuracy improves but more 
calculations are required. A typical good order of magnitude for s is 10. 
The purpose of PCGM is to obtain--in an efficient & reliable way--the 
solutions x(1), . . . , x(k) to the systems of equations Ax(i)=f(i), (i=1 
. . . k), corresponding to the loading vectors f(1), . . . f(k). 
The PCGM consists of k subsequent STAGES, each depending on the results of 
the previous stages. 
Each stage (i) produces the following: 
1) A set of vectors P(i), which is added to a previously-constructed set 
W(i-1), so as to form the new set of vectors W(i). The way by which the 
vectors in P(i) are constructed is provided below, but it is important to 
emphasize that as a result of this construction, ALL the vectors in W(i) 
come out to be MUTUALLY A-ORTHOGONAL. By this is meant that if u,v are any 
two vectors in W(i), then u.sup.t Av=0, where here "t" stands for the 
"transpose" operation for vectors. 
2) An approximate solution x(i) to the system Ax(i)=f(i). 
The details of this construction are as follows: The set of vectors W(1) 
and x(1) are constructed via PCG (Preconditioned Conjugate Gradients) 
applied to the matrix A and f(1), and using the above-described 
preconditioner. 
Having completed the i-th stage, the i+1 stage starts by construction of a 
vector y(i+1)--which is the numerically stable A-projection of x(i+1) on 
SPAN W(i)--and is done via the "PROJECTION" cyclic method, as described 
above. 
Then, a corresponding numerically stable "residual" r(i+1) is formed, which 
is in fact provided by: 
EQU r(i+1)=f(i+1)-Ay(i+1). 
The various vectors in P(i+1) are now constructed in the same way this is 
done in ordinary PCG, using the above described preconditioner,--with the 
following exception: 
Let v(1), . . . , v(i) each denote the LAST vector in P(1), . . . P(i) 
respectively. Then, once a vector Q' is constructed as one of the 
`conjugate-directions` obtained during the PCG process--it is REPLACED by 
the vector 
##EQU5## 
(In order to orthogonalyze Q' on W(i)) 
This new vector Q is now added to P(i+1), all the vectors of which are 
stored. It is THIS operation which assures the mutual A-orthogonality of 
all the vectors in W(i+1). 
The construction of an approximate solution x(i+1) then follows the 
procedure used in PCG. All the vectors in P(i+1) are now added to W(i), so 
as to form W(i+1), which is stored. 
The difference from the repeated, 'from scratch'approach is that the 
"conjugate directions" constructed at each stage serve also in obtaining 
the solutions for the next loading/stage. Furthermore, within each stage, 
the formation of "A-PROJECTIONS" and "CONJUGATE DIRECTIONS" are modified 
considerably, but these are necessary in order to obtain numerical 
stability and efficiency in computational speed. 
The PCGM algorithm ensures that the basic sequence W is orthogonalised by 
imposing that in the conjugate gradients procedure all the vectors of 
P(i)={p.sup.i.sub.1, p.sup.i.sub.2 . . . p.sup.i.sub.mi } are orthogonal 
only to the last vectors of P(k), k=1, 2, . . . , i-1. The procedure is as 
follows. 
x.sup.0 is set to y, r.sup.0 is set to r and p.sup.0 is set to K.sup.-1 
r.sup.0, where K is the preconditioning matrix. An index j is initialized 
to 1. 
For k=i the following steps are performed. 
(1) p.sup.j-1 is set as Q above 
(2) A value a.sup.j is set to (K.sup.-1 r.sup.j-1 r.sup.j-1)/(p.sup.j-1, 
Ap.sup.j-1), and x.sup.j is set to x.sup.j +a.sup.j p.sup.j-1 and r.sup.j 
to r.sup.j-1 -a.sup.j Ap.sup.j-1. 
(3) A value .beta..sup.j is set to (K.sup.-1 r.sup.j, r.sup.j)/(K.sup.-1 
r.sup.j-1, r.sup.j-1) and p.sup.j is set to Kr.sup.j +.beta..sup.j 
p.sup.j-1 and then to 
##EQU6## 
(4) If a stopping criterion is satisfied, then m.sub.k is set to j, and 
P(k) is set to , and, {p.sub.1.sup.k, p.sub.2.sup.k, . . . , 
P.sub.mk.sup.k } 
##EQU7## 
otherwise j is set to j+1 and the process repeated. 
The invention can be implemented as a set of tools for a general digital 
computer, to serve in structural-analysis packages, which yield fast & 
robust iterative solutions, especially for medium and large 3-D problems, 
and enable the Multi-Point Constraints and the several loadings cases to 
be handled efficiently when the number of Loading-Systems being solved 
does not exceed certain limits, which, with present-day machines seems to 
be around several tenths. These tools enable performance advantage to be 
taken of low-accuracy specs, and their storage requirements are very 
modest. 
The inner structure of most structural analysis packages makes them ideal 
candidates for incorporating a solver which operates according to the 
invention. The packages are large, and are constructed from several 
modules that communicate information through a common interface. More so, 
because typical problems take considerable machine resources, most 
packages would have restart facilities, and additional tools for enhanced 
end-user input. 
Even though the Solver module is only one among the several modules that 
makes the structural analysis package, it is the limiting one in terms of 
computer resources like memory, disk space and machine time. Incorporating 
a solver in accordance with this invention into the package will enable 
the user to execute larger problems using considerably less resources than 
those required by Direct Solvers.