.__
fi!B
cm
Computer methods
in applied
mechanics and
engineering
e
__
ELSEVIER
Comput.
Methods
Appl. Mech. Engrg.
150 (1997)
133-153
PZ: An object oriented environment for scientific programming
Philippe Remy Bernard Devloo
Faculdade de Engenharia Civil UNICAMP, C.P. 6021, Campinas SP 13081-970, Brasil
Paper dedicated to Prof. J.T. Oden
Abstract
An object oriented environment for scientific programming
is presented. The environment includes classes for matrix computations
(Tmatrix) and a set of classes for the implementation
of finite element algorithms (PZ). The PZ environment implements one and two
dimensional finite elements with arbitrary orders of interpolation and applicable to a variety of systems of differential equations. The
TMatrix environment includes an abstract matrix class interface and a variety of storage patterns and solution methods.
Introduction
The object oriented programming philosophy emerged in the decade of 1990 as a new style/philosophy
of
developing software. Its main benefits are (would be) the reusability of code, the modularization
of software
components
and therefore the reduction of the program development
cycle. While most of the software
developed today uses the object oriented programming philosophy, the dominant programming style in scientific
computing remains the structured programming style. This can be attributed to the following factors:
l Most scientists
use programming merely to verify the algorithms which they develop on a theoretical basis
l Switching
programming
language requires an investment
which takes several years. Object oriented
programming may be modular, but developing the (reusable) modules takes a lot of time.
Object oriented programming is more complex than structured programming. Therefore, it takes the students
more time to get accustomed and become productive.
l C +
+ compil.ers are not necessarily as available on computer systems as fortran compilers.
l If the class
structure of the object oriented environment is not well chosen, there is little advantage over
structured programming.
l Programs
developed using the structured programming style are numerically more efficient than the object
oriented counterparts.
It is the intent of this contribution to describe the object oriented environment
PZ applied to finite element
programming.
The term ‘environment’
is emphasized
because in more than just one application,
the
environment is an extendible software package which allows the user to implement any finite element algorithm
which fits within the proposed class structure.
The motivation for the development of the PZ environment found its roots in Texas, where the author was
involved in research on adaptive finite element algorithms. Adaptivity applied to compressible fluid flows, to
elasticity, to three dimensional
Poisson operators, etc. each generated a separate code even though the
underlying data structures were basically the same [l-6]. Therefore, an optimization which was implemented in
a particular code was not available in the other codes. The PZ environment started as a postprocessor for finite
element data with a graphical user interface [7] and has gradually evolved in the environment
presented
[g-11,26]. The code has been rewritten several times, to incorporate the experience of the author in object
oriented programming [ 121. The PZ environment relies on the TMatrix class library which implements a variety
of matrix storage patterns and solvers [ 131.
004%7825/97/$17&l
0 1997 Elsevier Science S.A. All rights reserved
PII SOO45-7825(97’100097-2
134
P.R.B. Devloo
I Comput. Methods Appl. Mech. Engrg. 150 (1997) 133-153
Many object oriented implementations
of the finite element method have been documented in the literature.
Most of the works dedicated to this topic are oriented to solve specific problems [14-18,291, to extend the user
interface to the finite element software [7,19,20,28,29] or to add artificial intelligence interface [21,22]. In this
contribution, the emphasis will lie in the description of how current algorithms can be implemented within a
unique class structure (in similarity with the objectives of Donescu and Laursen [23]), even if some of the
mentioned algorithms haven’t been implemented yet. It must be emphasized that the PZ environment is part of
an evolving project, in which the research efforts of graduate and undergraduate students will be incorporated.
My gratitude extends to Professor J.T. Oden to whom this paper is dedicated. The time I studied will never be
forgotten. Under his impulse and leadership I learned everything I know about the finite element method.
2. PZ: An object oriented environment for finite element programming
2.1. Scope of the project
The PZ environment intends to offer the user a class structure which is general enough to include a large
number of finite element algorithms. The environment
includes one- and two-dimensional
elements and
three-dimensional
elements are being developed. The class structure, which determines most of the characteristics of the environment, is modeled after the mathematical description of the finite element method. There is a
strong separation between the description of the geometry of the domain, the definition of the interpolation space
and the description of the bilinear form being approximated.
This implies that interpolation
spaces and
descriptions of differential equations are more important than the description of the implemented finite elements.
Thus far, the following problems are implemented within the PZ environment
A. Geometric elements
l One-dimensional
linear and quadratic elements.
l Two-dimensional
triangular and quadrilateral element with linear/ bilinear and quadratic / biquadratic
mapping.
B. Differential equations
l Beams
and frames.
l One-dimensional
convection diffusion equations.
l Two-dimensional
second-order systems of equations.
0 Plane elasticity.
l Flow
through porous media.
l Reissner-Mindlin
plate descriptions.
C. Interpolation spaces
l One-dimensional
elements with arbitrary order of interpolation.
l Two
dimensional quadrilaterals and triangles with arbitrary orders of interpolation.
Elements and differential equations are implemented separately. Therefore any differential equation which is
added to the environment
can automatically
be approximated
using the available interpolation
spaces.
Considering that the interpolation order of the elements is arbitrary, it also does not make sense to enumerate
linear, quadratic, cubic, etc. elements. Once adaptivity will be available, it will automatically be applicable to all
implemented differential equations. This implies that the environment will be able to unify adaptivity applied to
thermal problems, elliptic problems, elasticity, plate problems, etc.
2.2. Global class structure
The PZ environment is split in three class trees: the description of the geometry of the domain, the definition
of the interpolation space and the description of the differential equation. This mimics the development of a
finite element approximation:
l The domain
is partitioned in finite elements, which can be triangular or quadrilateral (two dimensions).
l A
set of functions is defined over each element, with CO continuity.
(All shapefunctions
have CO
continuity.)
P.R.B. Devloo I Comput. Methods Appl. Mech. Engrg. 150 (1997) 133-153
l
l
135
The set of shape functions is used to form a finite dimensional approximation of the original differential
equation.
The variational statement of the differential equation is used to define the residual and tangent stiffness
matrix.
2.3. Description
of the geometry
2.3.1. Class TGeoGrid
The TGeoGrid class is a container class which holds a set of geometric nodes (TGeoNod),
elements (TGeoEl) and boundary nodes (TGeoNodBc). At this stage elements, nodes and boundary
stored in binary trees (AVL trees) to minimize the search time for an element identified by it Id. The
grid contains the colmplete tree of refined elements (i.e. elements, sons, grandchildren, etc.). Therefore,
geometric elements does not necessarily partition the domain of the differential equation.
geometric
nodes are
geometric
the set of
2.3.1.1. Element connectivity
An efficient procedure is included to compute the connectivity of the elements which is compatible with any
grid dimension. A ‘brute force’ approach would be to compare the side nodes of each element with the side
nodes of all other elements in the grid. The operation count of such an algorithm is O(nel*) where nel is the
number of elements in the grid. The following algorithm computes the connectivity of the element in O(ne1)
operations:
Build an auxiliary integer vector of the size of the number of nodes in the grid. This vector is initialized with
- 1 and is named fclotprintvector. Each element tries to put the footprint of its side within the vector or verifies
whether a loaded footprint matches its own. An element loads a footprint of a side by loading its Id number in
the auxiliary vector at the indices corresponding to the Id’s of its sidenodes (Fig. 1).
Using the notion of footprint, the algorithm is described as follows:
1. while element :sides have been detected, loop over the elements
2. loop over the sides of the element
3. if the side is undefined
4. check the footprint of the element ith the footprintvector
5. if the footprint matches another element Id
l a connectivity
has been detected
l flag the element
side as defined
6. if the footprint matches the same element Id
l A complete
loop has been performed
l load
- 1 in the footprint of the element side
l flag the element
side as defined
7. if the footprint matches - 1
l load the element
side in the footprint
8. end if
9. end loop
10. end while
14
El5
10
ci
9
Fig. 1. Footprint
of a side of Element 5.
P.R.B. De&o
136
I Comput. Methods
Appl. Mech. Engrg. 150 (1997) 133-153
k
1
2
3
Fig. 2. Multipleconnectivityof a planargrid.
When building the connectivity, the following care has been taken:
l For planar
grids in three dimensions (modeling of plates and shells), an element can have more than one
connectivity on a side (Fig. 2). In this case element 1 will point to element 2. Element 2 will point to
element 3 and element 3 will point to element 1.
l The
algorithm will not function properly to find the connectivity of grids with both one-, two- and/or
three-dimensional
elements. The reason is that the footprint of a lower-dimensional
element will fit within
the footprint of the higher-dimensional
element, thus creating an apparent match. Connectivity
between
elements of different dimension is not defined as yet.
l The algorithm
can be optimized for grids with single connectivity: after a match has been detected between
two elements, the footprint of the matching side can be reset to - 1 (i.e. it is guaranteed that no other
element will match this footprint).
2.3.2. Class TGeoEl
The objects of type TGeoEl define the map between the element in
element. TGeoEl is an abstract class which defines the behaviour of
derived from it to implement the mapping from a deformed surface
Three methods define the mapping between the deformed element
l
l
l
void
void
its deformed configuration and the master
the derived classes. Different classes are
to a specific master element (Fig. 3).
and the master element
Jacobian
(DoubleAVec
&coordinate,
TFMatrix &jac, TFMatrix &axes) ;
X ( DoubleAVec &coordinate,
DoubleAVec &result ) ;
void
NormalVector
TFMatrix &axes) ;
The method Jacobian
(int
computes
side,
DoubleAVec
&loc,
DoubleAVec
&normal,
the Jacobian of the mapping of the master-element
TFMatrix
&jac,
to the reference space
defined by axes. Contrary to traditional finite element, the mapping is not necessarily defined between the
parameter space and the Cartesian space. The generalization of the Jacobian presented allows to define oneand two-dimensional
The one-dimensional
x=
elements in space. (Fig. 4).
Jacobian is computed as follows.
Given:
C xiPj(5)
I
Y
=c
z=
Y;PA&
I
and
E ZiVj(5)
I
Fig. 3. TGeoEl and derived classes.
Fig. 4. One- and two-dimensional
map
P.R.B. Devloo
The Jacobian
is computed
The vector Vl is defned
I Comput. Methods Appl. Mech. Engrg. 150 (1997) 133-153
as
as
A vector 42 needs to be supplied
externally
to define
A third vector V3 is computed as the cross product of Vl and V2.
Using these definitions, the values returned by the method Jacobian
value .I in the argument jac (a one by one matrix)
The vectors Vl , V2, V3 in the argument axes.
The two-dimensional
Jacobian is computed as follows:
as
137
are
such
Using the definition
of Vl, V2, V3, the Jacobian
matrix is equal to
The main use of this three-dimensional
extension of one- and two-dimensional
Jacobians lies in the simulation
of spatial thrusses and plates and shells where the orientation of the axes is determining for the computation of
the stiffness matrix. ‘This formulation also allows to compute heat conduction problems over surfaces in space or
over beams.
The method X(DoubleAVec
&coordinate,
DoubleAVec
&result)
; returns in the parameter result the
(x,y,z) coordinates of a point in parameter space
The method NormalVector
( int side,
DoubleAVec
&loc,
DoubleAVec
&normal,
TFMatrix
&jac,
TFMatrix
&axes) ; returns the surface/line Jacobian in the parameter jac along the side side at the location
lot in parameter space. It also returns the normal to the boundary in parameter normal and the axes as defined
in the method Jacobian.
The generality of the geometric element is based on the methods which declare its topology to the grid. As
such the following methods are defined :
l
int NumberOfNodes
( ) ; returns the number of nodes of the element
l TGeoNod*
NodePtr
lint
i) ; returns a pointer to the ith node of the element
The connectivity
of the element is defined by:
l
int
NumSides
l
int
NumSideNodes(int
( ) ; returns the number of sides of the element
side);
returns the number of nodes on a particular
side
138
P.R.B. Devloo I Comput. Methods Appl. Mech. Engrg. 150 (1997) 133-153
TGeoNod *SideNode (int side,
int nodenum) ; returns a pointer to the node nodenum along side
long SideNodeId(
int side,
int nodenum) ; returns the id of the node nodenum along side
l
int WhichSide(TLongVec
&SideNodeIds);
returns the side number which is connected to the
SideNodes and returns - 1 if no side is found
l
int SideIsUndefined(int
side);
returns 1 if the side has no defined boundary condition and no
element connected to it
l void
SetSide
(int i, int siden)
; set the sidenumber of the neighbour which is connected to side i
equal to siden
l void
SetConnectivity(int i, TGeoEl* el ) ; set the element neighbour along side i equal to el
l TGeoEl
*Neighbour
(short
is) ; returns the pointer to the neighbouring element along side is
l
shortNeighbourSide(short
is); returns the side number of the neighbour which is connected to the
element along side is
l
short
Bc (short
side) ; returns the boundary condition number associated with side. If the method
returns 0, there is no boundary condition associated with the side
The methods associated with adaptivity are:
l voidDivide(VoidPtrVec
&pv); divides the element (if necessary) and puts the resulting elements in the
point vector pv
l TGeoEl
*SubEl (short
is) ; returns a pointer to the subelement is
Each geometric element contains a material id which can be accessed by the method
l
long MaterialNumber
( ) ; returns the id of the material associated with the geometric element
The creation of computational
elements, based on the geometric element is made easier by the method
CreateCompEl();
l TCompEl
*CreateCompEl();
creates a computational
element compatible with the geometric element
(e.g. a triangular geometric element creates a triangular computational element, a brick geometric element
creates a brick computational element, etc.)
Note that the methods associated with the geometric element can be called independently
of whether the
element is one- two- or three-dimensional.
When including a new element, is it sufficient to implement the
above methods in the derived class for the environment to function properly.
l
l
2.3.3. Class TGeoNod
The class TGeoNod defines a point in three-dimensional
Euclidean space and associates
As such its main methods are
returns the Id of the node
l
long IdO;
l double
Coord (int i) ; returns the ith coordinate of the node
a unique Id with it.
2.3.4. Extensions
The separation of the geometric mapping from the interpolation allows to derive classes from the geometric
elements with special purpose maps, without affecting the interpolation space. This is not possible within the
concept of isoparametric elements where the function space used for mapping the geometry is necessarily the
same as the function space used for interpolation. High precision geometric mapping is often a must when high
order interpolation is used. As immediate extensions of the implemented elements are elements which map a
boundary exactly to a circle and mappings which use Coons patches to fit a boundary. Further extensions
include spherical and cylindrical shell elements which can be mapped exactly to their surface.
2.4. Dejinition
of the approximation
space
The approximation space is the set of functions which form the basis for the finite dimensional subspace in
which the finite element approximation
will be defined. As in traditional finite element technology, the
shapefunctions
are defined elementwise. Each computational element contains a set of shapefunctions
which,
through its node connectivity is Co continuous with corresponding
shapefunctions
of neighbouring elements.
The shapefunctions
are defined with respect to a master element. In order to compute the derivatives of the
P.R.B. Devloo
I Comput. Methods Appl. Mech. Engrg. 150 (1997) 133-1.53
shapefunctions
on the deformed configuration of the element, the computational
element to which it refers to compute the Jacobian at given points.
139
element
uses the geometric
2.4.1. Class TCompGrid
Like its TGeoGrid counterpart, the TCompGrid class is mainly a container class which holds lists of
l computational
Elements (TCompEl)
l degree
of freedom nodes (TDofNod)
l material
objects (TMaterial)
l boundary
condition objects (TBndCond)
l nodal
boundary conditions (TDofNodBc)
The computational
grid implements, beside the methods necessary to access the above lists:
void
ComputeNodeSequence(long
ElementId);
Mesh renumbering
scheme based on the Cuthill
McKee algorithm. The interface of the program with the Metis package is being investigated.
virtual
void ComputeNodElCon(
); method to compute the number of elements connected to each
degree of freedom node
void Assemble (TBlock &block,
TMatrix
&rhs) ; loops over the elements, computes the element
stiffness matrix and right-hand side and assembles the result in the matrices referred to by block and rhs
void LoadSo Lution ( TBlock &sol ) ; Loads the current solution into the degree of freedom nodes.
Makes sol the current solution, which can be used for post-processing
or for the computing of residuals
when approximating nonlinear differential equations
A computational
grid is always derived from a geometric grid because each computational element needs a
geometric counterpart to compute the Jacobian matrix. On the other hand, several computational grids can be
derived from a single geometric grid. Such can be the case where, using h-adaptivity, computational grids of
different levels are derived from the same geometric grid.
There are cases where computational elements need to be referenced from geometric elements (e.g. during the
refinement process). This behaviour is implemented by including a pointer to a computational element in each
geometric element. This implies that each geometric element can point to only one computational element at a
time. Therefore, only one computational
grid can be referenced from a geometric grid at any one time. The
computational
grid is defined as ‘loaded’ within the geometric grid when the computational
elements are
referenced by the geometric elements (Fig. 5).
The method to load the computational
grid in the geometric grid is
l void
LoadReferences{);
Contrary to the geometric grid, the elements of the computational
grid
partition the domain of the differential equation.
Grid1
I
1
1
I
I
1 TCompElld
Compu~omd
1 TCompEl
I
1 1
Grid2
Fig. 5. Computational
Grid1 is loaded in the geometric
Fig. 6. Class tree structllre of TCompEl.
grid.
I
TCompel2d 1
140
P.R.B. Devloo
I Comput. Methods Appl. Mech. Engrg. 150 (1997) 133-153
2.4.2. Class TCompEl and TDojNod
The computational element defines the interpolation functions over the element and contains the method to
integrate the stiffness matrix over the area of the element. As computational
elements come in one, two and
three dimensions, the TCompEl class is an abstract class from which specific classes are derived (see Fig. 6).
The main methods of the computational element are:
l void
CalcStiff(TElementMatrix
&ek, TElementMatrix
&ef); computes the element stiffness matrix and
right-hand side. The class TElementMatrix
is a structure which contains the stiffness matrix itself and
sufficient information to assemble itself within the global stiffness matrix without consulting the element.
(this feature is useful when transmitting the computed stiffness matrix over a network of computers)
l void
ApplyBc(TElementMatrix
&ek, TElementMatrix
&ef, TBndCond
&bc, int lado); computes the
stiffness matrix and right hand side contribution of a boundary condition identified by bc along side side.
All boundary conditions within the PZ environment are implemented as an integral contribution over the
boundary. This implies that Dirichlet boundary conditions are applied using the penalty method. The reason
for this option is that, when using high order interpolation with hierarchic shapefunctions,
it is difficult to
specify the value of the contribution of the individual shapefunctions which will approximate the imposed
value of the state variable.
l void ProjectFlux(TElementMatrix
&ek, TElementMatrix
&ef) = 0; computes the L’ projection matrix in ek
and the L’ projection of the fluxes on the interpolation functions in ef. This method implements the first
step of the Zienkiewicz and Zhu [24] error estimator. Considering that most problems in computational
mechanics have vector or tensor valued flux functions, the ef matrix will contain several columns.
l void EvaluateError(void
(*fp)( DoubleAVec &loc, DoubleAVec &val, TMatrix &deriv), double &true-error,
double &L2_error, TBlock *flux, double &estimate); Computes the true energy error, L’-error, and error
estimate based on the projected flux function and a function which computes the exact solution.
l void
Solution(DoubleAVec
&qsi, int var, DoubleAVec &sol, int &numvar);
computes the solution
associated with index var at the point in the master element space qsi.
p-Adaptivity is implemented by the method
l void SetInterpolationOrder(ShortAVec
&ord); Sets the interpolation order of the shapefunctions to ord. The
computational elements within the PZ environment support arbitrary orders of interpolation. Compatibility
of order of interpolation over the sides of the element is guaranteed by the knowledge of the neighbours of
the elements.
TCompElld
implements one-dimensional
computational
elements. Their order of interpolation is arbitrary
and the set of shape functions is derived from the Chebyshev functions. The first two shape functions are linear.
The next shape functions are obtained by taking the product of the quadratic function by a Chebyshev function
(see Fig. 7). This procedure has the advantage that it does not depend on the value of the orthogonal family of
polynomials : the continuity of the high-order polynomials is ensured by the multiplication
with the quadratic
shapefunction.
TCompE12d implements generic two-dimensional
computational elements: it implements all methods which
two dimensional computational elements have in common. Its main methods are:
l void
CalcStiff(TElementMatrix
&ek, TElementMatrix
&ef); computes the element stiffness matrix
l void
ApplyBc(TElementMatrix
&ek, TElementMatrix
&ef, TBndCond
&bc, int side); integrates the
boundary condition specified by the object bc along the side of the element
l void
Solution(DoubleAVec
&qsi, int var, DoubleAVec &sol, int &numvar);
computes the solution
associated with index var at the point qsi in coordinates of the master element and stores the result in sol.
numvar indicates the number of variables stored in sol.
TCompE12d computes/integrates
the stiffness matrix of an arbitrary two-dimensional
element by relying on
virtual methods which compute the shapefunctions,
integration points and Jacobian.
TCompElQ2d
implements
all methods which are specific to quadrilateral bi-dimensional
computational
elements. Its main method is:
l void Shape2d(DoubleAVec
&pt, TFMatrix &dphi, TFMatrix &dphi); computes the value of the shapefunctions and its derivatives at the coordinate pt defined on the master element. As in one-dimensional
elements, the high-order shapefunctions of the two-dimensional
quadrilateral element are obtained as the
product of the ‘standard’ finite element shapefunctions and high-order Chebyshev functions. There is no
limit with regard to the interpolation order of the element. For instance, the shapefunction along side 0 of
the quadrilateral element is defined as (Fig. 8) fid( 6,~) = &(~)@,,(~). High-order shape functions along
P.R.B. Devloo
I Comput. Methods Appl. Mech. Engrg. 150 (1997) 133-153
Fig. 7. Basic one-dimensional
141
shapefunctions.
Fig. 8. Quadratic shapefunction along side 0.
side 0 are defined as fiz( 5,~) = qGb(5, r])T,_,( 5) where T,,( 6) is the Chebyshev function of order p. Similar
procedure is applied to compute the shape functions of the internal shape functions.
One particular feature of the computational
element is that quadratic and higher-order shapefunctions
along a side are grouped within a single node. As such, node 4 of Fig. 8 will not only be associated with the
quadratic shapefunction, but also with all higher order shape functions denoted by $z( 5,~) (see Fig. 10). It
follows that the computational
element has a fixed number of nodes (see Fig. 9), but each node can be
associated with: a variable number of shapefunctions and corresponding number of equations. For a cubic
interpolation of a thermal problem, the number of variables associated with each node are
l Quadrilateral
[l, 1, 1, 1, 2,2, 2,2,4]
l Triangular
[I., 1, 1,2,2,2,
I]
Increasing the order of interpolation of the element only affects the number of equations associated with the
degree of freedom nodes. If the interpolation of the element is linear, the node still exists, but is associated
with zero shape functions. The set of equations of a node are called a ‘block’ of equations. When
numbering the equations of the computational
grid, block numbers are assigned to the nodes. Then, the
equation number corresponding to a particular block is obtained by summing the block sizes of the lowered
numbered blocks. The number of equations of the node is equal to the number of shape functions associated
with it times the number of state variables of the system of differential equations being approximated.
TCompElT2d implements all methods which are specific to triangular bi-dimensional
computational elements.
Its main method is:
computes the shapefunction of
l void
Shape2d(DoubleAVec
&pt, TFMatrix&phi,
TFMatrix&dphi);
a triangular element in master coordinates. The side shapefunctions
are obtained in the same way as the
9 degree of freedom
nodes
Fig. 9. Computational
quadratic
cubic
elements
7 degree of freedom
nodes
have a fixed number of nodes.
4th order
5th order
Fig. 10. Shape functions of different order.
6th order
142
P.R.B. Devloo I Comput. Methods Appl. Mech. Engrg. 150 (1997) 133-153
quadrilateral element, which guarantees the compatibility between both families of elements. The interior
shape functions are obtained by multiplying the cubic bubble with appropriate Chebyshev functions. Unlike
quadrilateral elements, where integration rules up to order 23 are available, triangular elements are limited
by the highest order which is documented in finite element textbooks (7th order) [25].
TDofNod class implements a degree of freedom node. The degree of freedom node holds the block number
associated with the block of equations it represents and contains a vector with the current solution of the node.
TDofNod also contains a variable indicating the number of elements connected to it. Note that in the case of
linear interpolation, the number of equations associated with the degree of freedom node may be zero.
2.4.3. Class TMaterial and TBndCond
Computational
elements compute shapefunctions,
their derivatives and integrate functions over the actual
(deformed) shape of the element. The computational element is not aware which simulation is being performed.
The coefficients of the differential equation and/or bilinear form are the responsibility of the material class.
Consider the following abstract boundary value problem (ABVP): find u E H’(0) such that
a(u, v) =f(v)
Vu E H’(O)
Most differential equations can be reduced to the above bilinear form. a(u, v) andf(v) contain integrals over the
domain 0 and over the boundary 80. It is assumed that the space H’(n) includes the Dirichlet type boundary
conditions. In systems of differential equations both trial and test functions will belong to product spaces.
The finite element approximation of the above ABVP can be written as: find uh E Vh(a) such that
a@,
vh) =f(vh)
V uh E Vh(J2)Vh(LJ)
a(uh, uh) can be written
a(uh, uh) =
In
C H’(n)
as
p(uh, uh, VL/, Vuh, x) da
+
G(u~, uh, x) dw
The first part of the integral is implemented in the method CalcStiff and the second part in the method ApplyBc
of the computational element. The first integral is decomposed as
P(uh, uh, Vuh, Vuh, x) da
= x
P(uh, uh, Vuh, Vuh,x) dR
el I .n,
Considering uh = Xi uih.(x, y) and assuming
be written as
linearity,
the (i, j)th contribution
to the element stiffness matrix can
where, for the sake of clarity, the dependence of P(fi., I+$,V&, V$, x ) on the location of the integration point is
omitted. The derived class from the TMaterial class computes the contribution to the stiffness matrix in function
of the values of shape function, their derivatives, the weight of the integration point and the value of the
determinant
of the Jacobian matrix. This arrangement
separates the definition of the interpolation
space
(TCompEl class) from the definition of the differential equation which is being approximated (TMaterial class).
The TMaterial class is an abstract class which defines the behaviour of the derived classes. Considering the
fact that this class is the only part of the program where the physical problem is described, its features include
both computation of stiffness contributions and post-processing.
Its main methods are
Service
l
l
l
l
methods:
Long Id(); returns the Id of the material
short NumberOfFluxes();
returns the number of fluxes which need to be projected on the approximation
space to make the Zienkiewicz and Zhu error estimator effective
short NumVariables();
Number of variables associated with the physical problem (e.g. 1 for Poisson
problems, 2 for 2D plane elasticity, 3 for plate problems, etc.)
virtual void SetForcingFunction(void
(*fp)(DoubleAVec
&loc, DoubleAVec &result)); identifies the
function pointer fp as the one used to compute the volume forces of the physical problem.
P.R.B. Devloo I Comput. Methods Appl. Mech. Engrg. 150 (1997) 133-153
143
Methods to compute the contribution to the stiffness matrix and right-hand side are
l void
Contribute(DoubleAVec
&x, DoubleAVec &sol, double weight, TFMatrix &axes, TFMatrix &phi,
TFMatrix &dphi, TFMatrix &ek, TFMatrix &ef); Computes the contribution to the stiffness matrix at an
integration point. In order to include nonlinear problems, the current solution at the integration point is
passed through the parameter sol.
l virtual
void ContributeBc(DoubleAVec
&LX,DoubleAVec &sol, double weight, TFMatrix &axes, Doub1eAVec &normal, TFMatrix &phi, TFMatrix &ek, TFMatrix &ef, TBndCond
&bc); Computes the
contribution of the boundary integral to the stiffness matrix and right-hand side.
Post Processing methods. Their main function is to compute a scalar or vector value of a variable associated
with the material in function of the value of the state variables and their derivatives. All variables which can be
computed through post processing are accessible through an index.
l int VariableIndex(char
*name); returns the variable index associated with the variable name
l int NumSolutionVariables(int
var); returns the number of variables associated with the variable indexed by
var. var is obtained by calling VariableIndex
l void
Solution(TFMatrix
&Sol, TFMatrix &DSol, int var, DoubleAVec &Solout, int &numvar); returns the
solution associated with the var index based on the value of the state variables (Sol) and their derivative
(DSol)
The current class tr’ee structure of TMaterial is shown in Fig. 11.
The above approach has the advantage that, whatever new technology is introduced at the level of the
geometric or computational elements will automatically be available for all physical problems modeled by the
TMaterial class.
TBndCond class is a structure to hold the data necessary to define a boundary condition. It holds a boundary
condition number (int), a boundary condition type (int), a square matrix and a vector of the size of the number
of state variables of the corresponding material.
2.4.4. Class TSuperEl
A super element groups a number of elements into a single element [30]. The TSuperEl class is derived from
both the computational grid class and computational element class. This implies that TSuperEl implements the
behaviour of both a computational
grid (assembly, renumbering
of the equations, adaptivity, etc.) and of a
computational
element (computation of the stiffness matrix, projection of the flux function). Computing the
stiffness matrix, the super element assembles the stiffness matrices of the individual elements into its global
stiffness matrix. That global stiffness matrix is then the element stiffness matrix of the grid to which the super
element belongs. The TSuperEl class implements
a method ReduceIntemalNodes();
which sets up the
datastructure to perform a static reduction of the internal nodes of the group of elements: the nodes are
effectively deleted from the grid to which the superelement belongs and the superelement only declares its
interface nodes to the outside grid. Static reduction amounts to the following algebraic manipulation:
consider a
matrix which is partitioned by its internal and external nodes
I
TMaterial
TMatldLin
‘I
I~II~I
Fig. 11. Class tree structure
of the TMaterial
class.
P.R.B. Devloo I Comput. Methods Appl. Mech. Engrg. 150 (1997) 133-153
144
then the condensed
[&,,I = [&I
stiffness matrix and right-hand
-
[K2,1K,,l-‘W,Zl
side are computed
as
and L.,fJ= If,1 - [~,,l[~,ll~‘[hl
Note that whereas [K,,] is sparse, [K,,] is full. When the CalcStiff method is called for the object of type
TSuperEl, the static reduction is automatically
performed and the [&,I and l$J matrices are returned. A
superelement can contain other objects of type TSuperEl as its elements. In that case the static condensation is
performed recursively without any user intervention.
2.45. p-Adaptivity
p-type adaptivity is included within the structure of the PZ environment [lo]. The TCompEl class includes a
method SetInterpolationOrder(DoubleAVec
&ord); where ord specifies the interpolation order in the direction of
the axes of the element. Continuity between elements of different interpolation order is easily implemented
when using hierarchic shapefunctions: the number of equations associated with a degree of freedom node on the
side of the element determines the interpolation order of that side (Fig. 12).
The same philosophy will be applied when developing three dimensional elements.
2.4.6. Extensions
The definition of the approximation
space is already fairly complete. However, some interfaces which are
defined but only partially implemented
are obvious extensions which need to be implemented
to turn the
environment more attractive.
. Implement
h-adaptivity.
Everything
in the interface of the environment
foresees the availability
of
h-adaptivity.
. Develop the transfer operators between grids. Interface of the PZ environment with multi grid algorithms.
. Implementation
of auto-adaptive grid refinement strategies.
. Implementation
of three-dimensional
elements.
. Compatibilization
of one dimensional,
two dimensional
and three dimensional
elements. In structural
analysis, it is common to combine one dimensional beam elements with two-dimensional
shell elements.
. Parallelization of the environment. The objective is to extend the basic grid classes to allow to read the grid
distributed over different processors, implement efficient solution algorithms and parallelize adaptivity over
the subdomains.
. Reduce the number of dynamic memory allocations. Tests indicate that the most probable limitation of the
environment implement large scale simulation problems is the number of dynamic memory allocations. The
Ci- + language hides dynamic memory allocation and release behind the class interface and induces the
programmer to abuse of this facility [12].
2.5. Finite element analysis
The TAnalysis calls the apropriate sequence of methods to perform a finite element analysis. Its input data are
l a pointer
to a computational grid
0rdX=S
ordy = 2
Fig. 12. Orders of interpolation
of side nodes and internal nodes.
P.R.B. Devloo
I Comput. Methods Appl. Mech. Engrg. 150 (1997) 133-153
a pointer to the global stiffness matrix
a solver associated with the global stiffness
The TAnalysis class calls in sequence:
l The node
numbering scheme
l assembly
method
l solution
algorithm
l post-processing,
procedure
145
l
l
matrix
2.5.1. Solving a linear problem
All that is needed to solve a linear problem is create an object of type TAnalysis. The constructor of the
TAnalysis class receives a pointer to a computational grid (i.e. it is not possible to create an analysis without a
computational grid). Optionally the user provides an externally created matrix to be used as the global matrix
(the default of the TAnalysis class is a full matrix). Optionally, the user chooses the solution procedure to be
used (the default is the Jacobi iterative method). Choose the variable the post-processing classes will output and
the plot filename. Call the ‘Run’ method. The following is an extract of the code:
TAnalysis an( &comp);
TFMatrix *sti.ff = new
TFMatrix(an.NumEquations(),an.NumEquations(),O.);
an.SetMatrix(stiff);
an.SolverReference().SetDirect(ELU);
VoidPtrVec scalnames(O),vecnames(O);
an.Run(scalnames,vecnames,“noplot.plt”);
2.5.2. Solving a time-dependent problem
The TTimeAnalysis
Class is derived from the TAnalysis
class but solves time dependent problems.
TTimeAnalysis
redefines the Run method, assuming the time dependent problem is solved using the following
algorithm:
K”“u” = F’“’
Mu”+ ’ =Lu”+F
It therefore solves the time-dependent
problem in three steps:
(1) Computing the initial state Kiniuo = Fin’
(2) Compute the L matrix
(3) Compute the bf and F matrix and loop for n timesteps un+’ = M-‘{Lu” + F}.
The same object of the TMaterial class is used to compute the three global matrices. For each problem the
TMaterial object is put in a different state.
Other types of analyses can easily be implemented. It usually is sufficient to redefine the Run method to call
the appropriate sequence of methods.
2.6. Numerical
integration
Any finite element program needs to perform a one-, two- or three-dimensional
integration. Integration rules
of any dimension are encapsulated under the class TIntRule (Fig. 13). From this base class integration rules for
specific domains are derived.
I
TIntRule
Fig. 13. Integration
rules for different
I
domains.
146
P.R.B. Devloo I Comput. Methods Appl. Mech. Engrg. 150 (1997) 133-153
The constructor of TIntRule class receives as parameter a vector containing the order of the integration rule in
the different directions. The order of the integration rule is defined as the order of the polynomial which can be
integrated exactly by the integration rule. The main methods of the TIntRule class are
l int NumPoints();
Number of integration points associated with the integration rule
l void Point(int
ip, DoubleAVec &pos, double &w); returns the position pos and weight w corresponding to
the integration point ip
l void
SetOrder(ShortAVec &ord); modifies the order of the integration rule to the order specified in ord
The interface restricts the integration rules to be static. In a future development, the TIntRule class will be
extended to include an interface to integrate a function autonomously and return a result. This then will open the
possibility to include adaptive integration rules. Examples and further documentation of this class can be found
on the web at http://www.cenapad.unicamp.br/-phil/integv.html.
2.7. General purpose
classes generated
with g + +
Any software package uses vectors, lists, maps, etc. The PZ environment initially used homemade classes of
these type of objects. Later, it was converted to use the basic classes available through the g+ + library package.
The reason for using these classes is that they are generally available (who does not have a g++ compiler?),
they are documented, reliable and using a simple script are applicable to any class type. Finally, software using
these classes may be shared with other scientists (respecting the gnu license agreement) without infringing
copyright laws.
2.8.
Pre-processing
The PZ environment does not include a pre-processor with any level of sophistication. It was opted to write
interface classes which are capable of reading files produced by external pre-processors. Currently, three classes
are available to build Geometric Grids.
l
TModuZef: Reads grids generated by the package Modulef. Generates grids with triangular elements
l
TGenGrid: Builds a uniform grid of an arbitrary number of elements on a rectangle. The grid can use either
rectangular or triangular elements.
l TSyslD:
Builds a one-dimensional
grid from a file which is compatible with the Code1 input file
documented in Oden, Becker and Carey.
Boundary condition data usually needs to be supplied externally because most pre-processors generate boundary
condition information on a node basis. The PZ environment
implements boundary conditions at the element
level. A list of boundary elements between two points can be found by the method.
l TGeoGrid::
GetBoundaryElements(int
NodFrom, int NodTo, VoidPtrVec &ElementVec, TIntVec &Sides);
returns all elements beweeen NodFrom and NodTo counter clockwise. This method uses the connectivity of
the elements. Therefore, Buildconnectivity
should be called to initialize the connectivity information before
calling GetBoundaryElements.
This method will only work for grids with 2-D topology. The current
version will only work for a grid with only one level (i.e. no refinement)
New formats can be easily added. The current formats only reflect the software packages available to the author.
2.9.
Post-processing
Considering the fact that the PZ environment is suitable to approximate differential equations in one, two and
three dimensions with arbitrary orders of interpolation,
interpretation
of the results is almost restricted to
graphical analysis. PZ uses external post-processing software rather than have a graphics library associated with
it. There are several reasons for taking this option:
l Not including
a graphics interface makes the environment much more portable (graphical user interfaces
are not portable).
l Many
users are accustomed to a post-processor and are reluctant to switch.
l The
PZ environment
is focussed on scientific computing.
Adding a graphical user interface with
post-processing capability slows the development cycle down.
To visualize the results generated by PZ, output filters are implemented for three graphics package:
P.R.B. Devloo I Comput. Methods Appl. Mech. Engrg. 150 (1997) 133-153
147
(1) MView, deve:loped by the research
group TecGraf associated with PUG/Rio
de Janeiro (Catholic
University of Rio de Janeiro).
(2) View3d, developed by Prof. Fernando L. Ribeiro from COPPE/Civil UFRJ (Federal University of Rio de
Janeiro).
(3) DX@, Data Explorer developed by IBM@
Other file formats can easily be implemented.
The graphics output classes of the environment
superpose a
uniformly refined grid on top of the computational
grid. This uniform refinement allows to generate graphics
output compatible with the interpolation order of the elements. Graphical elements and graphical nodes are
grouped in a graphical grid.
2.9.1. TGrafGrid
TGrafGrid contains a map of graphical elements and nodes used to generate output to the graphics file which
will be read by the ‘external’ post processor. A graphical grid is built from a computational grid by creating a
graphical element from each computational
element and a graphical node from each degree of freedom node
(see Fig. 14). TGrafGrid is a virtual class from which output specific classes are derived (see Fig. 15).
The derived graphical grids draw the appropriate header files and drawing commands for the target graphics
package. Between the headers and drawing commands, the graphical nodes and elements spool their content to
the output file.
2.9.2. TGrafel
The graphical element TGrafel is the logical representation
of a uniformly
refined element. Like its
computational
element counterpart, TGrafEl is the base class of a family of topologically different graphical
elements (see Fig. 16). At this stage only quadrilateral and triangular elements are implemented.
The number of uniform refinements is stored as a parameter of TGrafGrid. The graphical element does not
store the subelements or their connectivity. The numbering of the nodes of the subelements and their location is
computed. Like the computational
elements, graphical elements have a fixed number of graphical nodes.
Graphical nodes use the graphical element to which they belong to compute their location and/or value.
When drawing the grid to a file, the following steps are taken:
Write the graphical nodes to the output file
l each node
corresponds to one or more physical nodes
l the graphical
node consults the graphical element to know the location of the physical nodes
Write the element connectivity to the output file
l the graphical
element computes the connectivity of the subelements based on the sequence number of the
graphical nodes
Write the solution to the output file node by node
l the
graphical node consults the graphical element to compute the solution. The graphical element
consults the computational element to compute the solution at a point within the parameter space of the
master element.
Using this scheme, each element can be represented by an arbitrary number of
subelements without any additional use of memory (e.g. all data of the refined data structure is computed, not
stored).
I
I
I
,
TGtzSrld
1TDXGrafGtid 1 1TWGrafGrld
Fig. 14. Sequence
of TGrafGrid,
TCompGrid
Fig. 15. Class tree structure
and TGeoGrid.
of TGrafGrid.
I
I
1 1TWdGmfGrid 1
148
P.R.B. Devloo
I Comput. Methods Appl. Mech. Engrg. 150 (1997) 133-153
6
3
I-: I-------C-. rL --m-i-2
I
r
I
r-
I’
k _’
____--_
I’: I
I
8
;/:
I
I5
I'
II
I
Fig. 16. Graphical
I
,
"
II
II
I
I_'1
1.. ,,.---d
1:: I
r-l
c,
I---
‘7
4,
1,
---‘:
__-_--
4
0
Fig. 17. Quadrilateral
2
1
element class tree structure.
graphical
element with 9 graphical
nodes.
2.10. Future developments
The PZ environment is showing itself to be an excellent testbed for algorithm development. Several areas are
being developed by the author, colaborators and students:
l h-Adaptivity:
inclusion of h-adaptive capabilities for the geometric and computational grids. Development
of error estimators and automatic grid refinement strategies.
l Parallelization
of the environment. In a parallel development, the author develops an environment for the
development of parallel scientific software (OOPAR). The objective of this line of research is to adapt the
classes of the PZ environment to the OOPAR environment.
l Study
of flexible mechanisms. Inclusion of dynamic analysis and use of a co-rotational approach of finite
element approximations
to simulate flexible mechanisms.
l Numerical
efficiency of the environment. Use the environment for large scale simulations. The intent is to
reduce the number of dynamic memory allocations to a minimum : no memory allocation during the
computation of the global stiffness matrix and allocation of elements and nodes in groups. Try to reach the
goal of less than one allocation per element and/or node.
l Use the sequence
of refined grids for multigrid acceleration. Development of the transfer operators of high
order interpolations.
Coupling of the environment with the MADPACK package of multigrid algorithms.
l Development
of iterative algorithms to efficiently solve substructured
grids. Substructuring
is already
included in the environment under the form of the TSuperEl class. The intent is to use substructuring as the
basis of parallel algorithms and to develop efficient/parallelizable
iterative solvers to reach efficient
parallelized codes.
4. Conclusion
The PZ environment is a set of classes which through mutual interaction constitute a powerful environment
for the development of finite element algorithms. The environment
is entirely based on the object oriented
programming philosophy which permits to combine generality and complexity with a relatively simple user
interface. The PZ environment intends to encompass a large number of finite element algorithms and to reach
this goal using a single class structure. The strong separation of geometrical approximation,
definition of the
interpolation
space and definition of the differential equation being modeled has allowed unprecedented
generality.
P.R.B. Devloo I Comput. Methods Appl. Mech. Engrg. 150 (1997) 133-153
149
Acknowledgement
The author gratefully acknowledges
the support of CNPq through grant No 522.135/94-3
and FAPESP
through grants No 94/4159-l,
96/3652-l,
96/7989-O. The author also extends his gratitude to CENAPAD and
its support staff and colleagues and students which contributed to the development of this work.
Appendix
A. TMatrix-A
set of matrix classes
A. I. Introduction
Any program of scientific computing and finite element computations
in particular uses matrix objects
extensively. Therefore, in a separate effort, a set of matrix classes was developed to facilitate the algebraic
manipulation of matrix object. The efficiency of the PZ environment is a direct consequence of the efficiency of
the TMatrix class I.ibrary. Therefore, special care was taken to give the user the option to perform algebraic
operations between matrix objects without allocating memory dynamically
[12]. Full documentation
with
examples of the TMatrix library class can be found at http: // www.cenapad.unicamp.br
/ -phi1 / tmatrix.html.
A.2. Global design philosophy
A TMatrix object is the representation of a linear transformation
between two vector spaces. The TMatrix
object should be able to perform this linear transformation efficiently, independently
of its internal representation. Even though any linear transformation has a matrix representation, the individual elements of the TMatrix
are not necessarily efficiently computed. Certain TMatrix allow their elements to be modified (e.g. full matrix,
banded matrix, sparse matrix), but others do not (e.g. element by element matrix storage pattern).
The TMatrix class is an abstract class which defines the methods which derived classes need to implement in
order to function properly. The following methods need to be defined:
l int PutVal(int
row, int col, REAL val ); // stores val at the location row, co1
l REAL
&GetVal(int row, int col); // retrieves the value at the location row, col
If these methods are not defined, the following method needs to be implemented for the derived class to make
sense:
l void
Multiply(TFMatrix
&A,TFMatrix &-es, int opt = 0); // computes the multiplication
of the current
object with the full matrix A and puts the result in res. The opt parameter indicates whether the
multiplication
should be performed with the transpose of the object or not.
The (abstract) TMatrix class implements the Multiply method using the methods GetVal and PutVal. It is
advisable to reimplement the Multiply method as well to obtain any efficiency. Multiplication
is only defined
between an abstract TMatrix object and objects of type TFMatrix (Full matrix storage pattern). This implies that
multiplication between matrix objects of different types is not defined. Besides the above methods, the TMatrix
class implements the following methods which declare its dimension, state, etc to the environment:
l virtual
int IsSimetric();
// indicates whether the storage pattern of the matrix is symmetric or not
l inline
int IsSquare(); // returns 1 if the matrix is square
l int IsDecomposed();
// returns the number of the decomposition type of the matrix and 0 if the matrix is
not decomposed
l void
SetIsDecomposed(int
val); // modifies the flag indicating the decomposition type
A.3. Implemented
The TMatrix
l TSpMatrix:
l TFBMatrix:
l TFMatrix:
l TMatRed:
l TSMatrix:
storage patterns
class is the base class from which a variety of matrix types are derived (see Fig. A.l).
implements a sparse matrix storage pattern
i:mplements a non-symmetric
banded storage pattern
implements a full matrix
implements a matrix which implements the algorithm of static reduction
base class for symmetric matrices
P.R.B. Devloo I Comput. Methods Appl. Mech. Engrg. 150 (1997) 133-153
150
Fig. A.l.
l
l
l
l
TMatrix class tree.
TSFMatrix: implements a full matrix with symmetric storage pattern
TSBMatrix: implements symmetric banded storage pattern
TSSpMatrix: implements a sparse symmetric storage pattern
TSkylMatrix: implements the skyline storage pattern
A.4. Direct
solvers
The TMatrix class implements three direct solution algorithms:
(1) LU decomposition:
this scheme only applies to non-symmetric
storage patterns
(2) Cholesky Decomposition:
this scheme only applies to positive definite matrices. When applied to
non-symmetric
storage patterns, only the upper triangular part is affected.
(3) LDL’ decomposition:
Crout decomposition
scheme, which is applicable to non-singular
symmetric
matrices.
All the direct solvers are implemented as methods of the TMatrix class, which uses the methods GetVal and
PutVal to access and modify the elements of the object. This is obviously not an economic approach when
applied to banded matrices. Therefore, the derived classes reimplement
the decomposition
schemes, using
pointer arithmetic to increase the efficiency.
A.5.
Iterative
solvers
The TMatrix class library is compatible with the Templates header files ([27]), and therefore implements the
iterative methods which are described therein. The templates header files were modified to avoid the use of
dynamic memory allocation: the header files distributed with the templates software on netlib use the ‘*’
overloaded operator to perform matrix multiplications.
The use of this operator implies the use of temporary
objects and dynamic memory allocations. The iterative solvers which are interfaced with the TMatrix class
library are
( 1) Jacobi iteration
(2) SOR iterative scheme
(3) SSOR iterative scheme
(4) Preconditioned
conjugate gradient
(5) Preconditioned
Richardson scheme
(6) GMRES iterative scheme
The Jacobi, SOR and SSOR schemes are reimplemented
in most subclasses to increase efficiency. The
conjugate gradient, Richardson and GMRES scheme only use the Multiply method to implement
their
behaviour.
P.K.B. Dcvloo I Camput. Methods Appl. Mech. Engrg. 150 (1997) 133-153
A. 6. Block structured
151
matrices
In finite element computations, most matrices are block oriented. A block of equations is coupled such that
the reordering of the equations within a block do not affect the structure of the matrix. As such the block size of
plane elasticity problems is at least two. When using high-order elements, all equations associated with internal
shapefunctions
(which are zero on the boundary of the element) form a block of equations. This feature is
explored within the PZ environment by the concept of degree of freedom node which groups all these equations
with a single node. To facilitate the use of blocks of equations, the TBlock class allows the user to logically
represent a matrix object as blocked. The TBlock constructor associates a TMatrix object with a vector of blocks
of different sizes. The main methods of the TBlock class are
l int SetNumBlocks(int
num_of_blocks);
// Sets the number of blocks of the TBlock object
l int Set(int
index, int dim, int pos = - 1 1; // Sets the dimensions of the block index to dim, indicating that
the block starts at position pos.
9 int SetAll(TlntAVec & dimensions); // modifies the dimension of all blocks by the dimension indicated by
the vector dimensions
l REAL&
GetVal(int bRow, int bCo1, int r, int c); // returns a reference to the (r,c) element of the block
(bRow,bCol)
l
int PutVal(int
bRow, int bCo1, int r, int c, REAL value);
//
puts value at the (r,c) element
of the block
(bRow,bCol)
A TBlock object can be associated with any matrix object. It does not affect the structure or storage pattern of
the matrix but allows to pick an element of the matrix in a blockwise fashion.
A.7.
The TSolver
class
The TSolver class allows the user to identify a solution process with an object. It associates one of the
implemented solution methods with a matrix object. This class is particularly useful to parametrize solution
processes. As such the preconditioner
of the conjugate gradient method is a TSolver object, which can be a
SSOR solution process, a number of Jacobi iterations, or any other kind of symmetric solution process. Its main
methods are
l void
SetMatrix.(TMatrix
*Refmat); // associates Refmat with the TSolver object
l void SetSOR(int
numiterations, REAL overrelax, REAL tol, int From&-rent);
// Sets the SOR method as
the solution process
l void SetSSOR(int
numiterations, REAL overrelax, REAL tol, int FromCurrent);
// Sets the SSOR method
as the solution process
l void SetJacobi(int
numiterations,
REAL tol, int FromCurrent);
// Sets the Jacobi method as the solution
process
l void SetCG(int
numiterations, TSolver &pre, REAL tol, int FromCurrent);
// Sets the conjugate gradient
method as the solution process using the TSolver object pre as the preconditioner
l void
SetDirect (DecomposeType
decomp); // Sets the solution process as the direct solver indicated by
decomp (ELU, ECholesky or ELDLt)
l void
Solve(TFMatrix
&F, Tl3latrix
&result, TFMatrix *residual = 0); // Solve the system of equations
using F as right-hand side, storing the result in result and returning the residual in residual if a non zero
pointer is provided.
A.8
Additional
features
for the full matrix
class
The purpose of the TFMatrix class is to implement small full matrices to be used as local objects in the code.
Whereas FORTRAN includes full matrix objects as a construct of the language, C and C+ + do not support this
feature. As such the TFMatrix class simulates its FORTRAN equivalent. The TFMatrix class documents and
uses the FORTRAN storage pattern for full matrices. Therefore, its elements can be accessed by pointer
arithmetic. As an additional feature the TFMatrix class can work with externally supplied storage space,
allowing to create TFMatrix objects on the stack, without any dynamic memory allocation. The numerical
efficiency of the PZ environment depends mainly on the efficiency of the TFMatrix class.
152
P.R.B. Devloo
I Comput. Methods Appl. Mech. Engrg. 150 (1997) 133-153
A.9. Future developments
(1) The TMatrix objects do not have any eigenvalue
methods associated with them. The next project on the
TMatrix class library is to include eigenvalue methods as well.
(2) The sparse matrix storage pattern which is implemented is inefficient. Include also other sparse matrix
storage patterns.
(3) Only the TFMatrix class is suitable for using small matrix objects. The extensions included in the
TFMatrix class need to be extended to the TSFMatrix class.
(4) Interface the TMatrix class library package with the BLAS routines. Implement the BLAS routines
efficiently and create compatible header files. This means improve on the BLAS distribution of Netlib,
which is a mere f2c run of the fortran version.
References
[ 11L. Demkowicz, P.R.B. Devloo and J.T. Oden, On an h-type mesh refinement strategy based on minimization
of interpolation errors,
Comput. Methods Appl. Mech. Engrg. 55(1-2) (1986) 63-87.
[2] J.T. Oden, L. Demkowicz, T. Strouboulis and P.R.B. Devloo, Adaptive methods for problems in solid and fluid mechanics, I. Babuska.
O.C. Zienkiewicz, J.P. d S.R. Gago and A. de Oliveira, eds., Accuracy Estimates and Adaptive Refinement in Finite Element
Computations (John Wiley and Sons, Ltd., London, 1986).
[3] T. Strouboulis, P.R.B. Devloo and J.T. Oden, Adaptive finite element methods for the analvsis of inviscid comoressible flow 1: Fast
refinement/unrefinement
and moving mesh methods for unstructured meshes, Comput. Methods Appl. Mech. ‘Engrg. 50(3) (I 986)
327-362.
[41 J.T. Oden. T. Strouboulis and P.R.B. Devloo, Adaptive finite element methods for high speed compressible Hews, in: R.H. Gallagher,
RI. Glowinski, P.M. Gresho. J.T. Oden, O.C. Zienkiewicz, eds., Finite element in fluids V7 (John Wiley and Sons Ltd., 1988)
223-240.
PI P.R.B. Devloo, J.T. Oden and P. Pattani, An h-p adaptive finite element method for the numerical
Comput. Methods Appl. Mech. Engrg. 70 (1988) 204-235.
simulation
of compressible
flow,
I61 P.R.B. Devloo, A three dimensional adaptive finite element strategy, Comput. Struct. 38(2) (1991) 121-130.
(Phase 1) a system independent
[71 P.R.B. Devloo and J.S.R. Alves, An object oriented approach to finite element programming
windowing environment for developing interactive scientific programs, Adv. Engrg. Software Workstations 14( 1) (1992) 161-169.
WI J.S.R. Alves and P.R.B. Devloo, An object oriented approach to finite element programming, Matemkica Aplicada e Computational
lO(3) (1992) 17-26.
191 P.R.B. Devloo and J.S.R. Alves. On the development of a finite element program based on the object oriented programming
philosophy, Proc. First European Conference on Numerical Methods in Engineering, Brussels, Belgium (Sept. 1992) 39-42.
of the p-adaptive finite element methods using the object
[lOI P.R.B. Devloo. C.A. Magalhaes and A.T. Noel, On the implementation
oriented programming philosophy, Proc. Int. Congress on Numerical Methods in Engineering and Applied Science, Conception, Chile
(Nov. 1992) 137-143.
I?R.B. Devloo, On the development of a finite element program based on the object oriented programming philosophy. OONSKI 93,
Proc. Object Oriented Numerics Conference, Sunriver, OR (Apr. 1993) 113-123.
1121 P.R.B. Devloo, Efficiency issues in an object oriented programming environment, Proc. Second Int. Conf. on Computational Structures
Technology, Athens, Greece (Aug. 1994) 147-152.
[I31 M.L.M. Santana and P.R.B. Devloo, Object oriented matrix classes, Proc. XVI CILAMCE, Italian Group of Computational Mechanics
and Ibero-Latin American Association of Computational Methods in Engineering, Padova, Italy, Sept. 1996.
MTT-S
[I41 M. Yuri, M. Lea1 and L.A. Bermudez. MEF system: an object oriented finite element package. Proc. 1995 SBMO/IEEE
International Microwave and Optoelectronics
Conference. Part 2 (of 2).
rw G. Archer, C. Thewalt and G.L. Fenves, New software architecture for finite element analysis, Comput. Civil Engrg. 683-689.
1161 A.P. Peskin and G.R. Hardin, Object oriented approach to general purpose fluid dynamics software. Comput. Chem. Engrg. 20(8)
(1996) 1043-10.58.
[I71 X.A. Kong and D.P. Chen. Object oriented design of FEM programs, Comput. Struct. 57( 1) (1995) 157-166.
field
[I81 E.J. Silva, R.C. Mesquita, R.R. Saldanha and P.F.M. Palmeira, Object oriented finite element program for electromagnetic
computations,
Proc. 9th Conf. on the Computation of Electromagnetic
Fields (COMPUMAG’93).
Miami, FL. USA, IEEE Trans.
Magnet. 30(5) pt 2 (Sept. 1994).
[19] H. Ohstubo, Y. Kawamura and A. Kubota, Development of the object oriented finite element modeling system-MODIFY,
Engrg.
Comput. 9(4) (1993) 187-197.
[20] Jianing Ju and M.U. Hosain, Finite element graphic objects in C+ +, J. Comput. Civil Engrg. lO(3) (1996) 258-260.
[21] Th. Breitfeld and B. Kroeplin, Expert System for the verification of finite element calculations, Proc. 1996 4th Int. Symposium on
Assessment of Software Tools. IEEE, Los Alamitos, CA, USA, pp. 18-23.
[22] Th. Zimmerman, Y. Dubois-Pelerin and P. Bomme, Object-oriented finite element programming (part I): governing principles, Comput.
Methods Appl. Mech. Engrg. 98(2) (1992) 291-303.
[Ill
153
P.R.B. Devloo I Comput. Methods Appl. Mech. Engrg. 150 (1997) 133-153
[23] F! Donescu and T.A. Laursen, Generalized object oriented approach to solving ordinary and partial differential equations using finite
elements, Finite Elcm. Anal. Des. 22(l) (1996) 93-107.
[24] O.C. Zienkiewicz and J.Z. Zhu, A simple error estimator and adaptive procedure for particle engineering analysis, Int. J. Numer.
Methods Engrg. 24 (1987) 333-357.
[25] T.J.R. Hughes, The finite element method, Linear Static and Dynamic Finite Element Analysis (Prentice Hall International Editions.
1987).
[26] J.S.R. Alves and P.1R.B. Devloo, Object oriented programming in scientific computing, the beginning of a new era, Engrg. Comput. 8
(1991) 81-87.
[27] R. Barrett, M. Berry, T.F. Ch, J. Demmel, J. Donato, J. Dongarra, V Eijkhout, R. Pozzo, C. Romine and H. Vandervorst, Templates for
the solution of linear systems: building blocks for iterative methods (SIAM, 1994.
[28] Yao Zheng, R.W. Lewis and D.T. Gethin, FEView: an interactive visualization tool for finite elements, Finite Elem. Anal. Des. 19(4)
(1995) 261-294.
[29] A. Cardona, I. Klapka
and M. Geradin,
Design of a new finite element programming
365-381.
[30] P.R.B. Devloo and M.L.M. Santana, Desenvolvimento
LATCYM 96, Floriarmpolis, SC, Brasil, Nov. 1996.
de algoritmo
de subestructuracao
environment,
para elementos
Engrg.
Comput.
ll(4)
(1994)
finitos, Proc. ENCIT 96 and