- Email: [email protected]
Object-oriented Fortran 90 P-adaptive finite element method
Advances in Engineering Software 33 (2002) 461–468 www.elsevier.com/locate/advengsoft
Object-oriented Fortran 90 P-adaptive finite element method J.E. Akina,*, M. Singhb a
Rice University, Houston, TX 77005, USA b COADE Inc., Houston, TX 77070, USA
Abstract Over the last decade, there has been an increased awareness of the benefits of employing Object-Oriented (OO) design and methodologies for development of software. Among the various languages available for OO development, Fortran 95 has some clear advantages for scientific and engineering programming. It offers features similar to other OO languages like Cþ þ and Smalltalk as well as extensive and efficient numerical abilities. This paper will describe the OO design and implementation of P-adaptive finite element analysis (FEA) using Fortran. We will demonstrate how various OO principles were successfully employed to achieve greater flexibility, easier maintainability and extensibility. This is helpful for a complex program like an adaptive finite element implementation. q 2002 Civil-Comp Ltd and Elsevier Science Ltd. All rights reserved. Keywords: Object-oriented design; P-adaptive finite element analysis; Fortran
1. Introduction The idea of object-oriented programming can be traced back to Simula 67 language. It was realized that OO techniques could solve the problems of code reusability and maintainability encountered in engineering software development. Consequently, there has been considerable interest over the last decade in this alternate philosophy of programming. We can view OO programming as an analog model of real world phenomenon. A program thus builds an environment containing objects (such as beams, floors) with attributes analogous to real world counterparts (such as weight, lengths and properties) and analogous behavior (deformation under loads). An object can also be an abstract entity (such as matrix, vector, etc.). Akin [2], Szymanski [14,15] and his student Norton [11] have elaborated on the implementation of the OO concepts in F90. Rehak et al. [12] and Forde et al. [8] first applied OO principles to finite element programming. Zimmermann et al. [6,7,19] in a series of papers presented two OO finite element toolkits, implemented in Smalltalk and Cþ þ . Dubois-P‘elerin et al. [5] proposed further modularity by delineating two kinds of behavior of a finite element software: analysis type (solvers used, static or dynamic analysis) and domain information (element, grid points, boundary conditions). Liu et al. [9] introduced an OO * Corresponding author. E-mail address: [email protected] (J.E. Akin).
implementation of adaptive finite element and finite volume methods. They used the Cþ þ language to implement the program, which performed h-refinement and offered a choice of direct or iterative solvers. While much of the previous OO work is in Cþ þ , Smalltalk and Pascal, no adaptive finite element analysis (FEA) system has been reported in Fortran 90 (F90). Fortran 90 and later 95 introduced dramatic new features over Fortran 77. These features helped in flexible program design, dynamic memory management and powerful built in numerical capabilities. The Type command, used for creation of user-defined types in the Module construct allows easier implementation of the OO style of programming. We will illustrate various concepts of OO programming in an FEA flavor with examples from the actual implementation. Since its birth in the 1940s, FEA has matured and is now employed in a wide range of problems, from solid mechanics to heat transfer, fluid mechanics, acoustics, and electromagnetism. The main advantage of FEA is its ability to model complex and arbitrary shaped domains. It also supports modeling of general boundary conditions and non-homogeneous anisotropic materials. To achieve a greater degree of confidence in an FEA, adaptive techniques have been developed. At the core of an adaptive FEA program is the error-estimator, which provides numerical estimates of the error in the solution. Using this error information the mesh is then accordingly modified. This cycle is repeated until an error less than
0965-9978/02/$ - see front matter q 2002 Civil-Comp Ltd and Elsevier Science Ltd. All rights reserved. PII: S 0 9 6 5 - 9 9 7 8 ( 0 2 ) 0 0 0 4 8 - 0
462
J.E. Akin, M. Singh / Advances in Engineering Software 33 (2002) 461–468
Fig. 1. A Point class demonstrating various object-oriented concepts.
the specified threshold is achieved. There are three basic approaches for mesh adaptation: H-refinement, which subdivides the selected elements into smaller elements; P-refinement, which reduces the error in the solution by increasing the polynomial order, p; of the interpolating functions; and R-refinement, which relocates the position of nodes in a finite element mesh. This paper will describe the implementation of the P-adaptive FEA. More details are available in Refs. [1,13]. An adaptive FEA program places high demands on software development for effectively managing the complexities resulting from error-estimation and mesh
adaptation which requires modifying the underlining data structures. Traditionally, FEA programs have been programmed procedurally in Fortran 77 [1]. In a procedural style, also termed as water-fall technique [10], a program is viewed as a set of algorithms. Therefore, the program depends closely on the algorithm being implemented, creating a rigid framework where a change becomes costly and sometimes unmanageable unless the algorithm components are highly modular. OO programming is capable of overcoming some of the limitations of procedural programming, but does so at additional cost. We will demonstrate this in subsequent sections.
J.E. Akin, M. Singh / Advances in Engineering Software 33 (2002) 461–468
463
Fortran 90 supports the following types of inheritance:
Fig. 2. Example of inheritance by Composition.
2. Main OO features in F90 2.1. Class Fortran 90 allows the creation of user-defined types to encapsulate the related concepts into a single form. The type thus created is associated with the related methods within a module block, which provides an interface for communication with external objects. This is similar to the notion of Class in Cþ þ and Java. A module can contain more than one logically related class. Fig. 1 shows a Point class within a module. This class encapsulates the coordinates and spatial dimensions of the point. The USE statement provides higher class objects in different modules with access to the methods of lower class objects (this is not merely a text substitution, as done by the INCLUDE command). This can be seen in the case of the Point class located in Point Module, as the Point Class gets access to the data and methods in the Dbl_Precision_Module. We can declare the contents in a module as private or public to provide concepts the outside user access only to relevant information and to hide the complexity within the module. 2.2. Polymorphism It is the property of OO languages to allow the use of identical function names for logical similar methods in different types. This property can be seen in some intrinsic functions like int(a ), where a can be of integer, real or complex type. This can be extended to user-defined objects using generic programming, operator overloading and the interface construct. For example, in Fig. 1 Init_Point is bound to Init using the interface construct. This binding can be performed in all the classes in the program so that the user only utilizes one call to Init instead of individual calls to Init_Element, Init_Point, etc. At compile time, the compiler uses the information about the argument types to uniquely determine which procedure to execute. Similarly, operators can be overloaded for user-defined types (already existing for implicit types), and additional operators can be defined. Fig. 1, shows how the ‘ ¼ ¼ ’, equality operator is defined for the Point class. The assignment operator ‘ ¼ ’ between two objects of the same type is implicitly implemented by the language.
† Inheritance by granting access. This is the simplest form of inheritance. The use statement in Fortran 90 allows objects in one module to access the public methods and data in the base module. † Inheritance by composition. Also known as ‘has-a’ inheritance, this type of inheritance uses other classes to build larger and more complex classes. As shown in Fig. 2, the Gauss_Type class contains an instance of Point class to define its position in the space. The Gauss_Type class is granted access to contents of the Point class’s public declarations and USE statements. † Inheritance by subtyping. This type of inheritance can be defined as ‘is-a’ inheritance. The ‘is-a’ inheritance is explicitly implemented in Cþ þ and Java languages, while in Fortran 90 it has to be simulated by the programmer using composition. Szymanski et al. [15] introduced the technique to implement this type of inheritance in Fortran 90. It can add dynamic capabilities to the program, but at the added expense of a small increase in the complexity of the program. The array syntax and intrinsic array operators in F90 provide advanced numerical capabilities very similar to those in the Matlab environment. Using the assumed shape array, automatic array, and the SIZE command, the procedures can be made independent of the length of the array arguments, lending more flexibility to the code. Fortran 90 also provides some safety features; for example, pointers can be checked to verify that they are associated (pointing to something valid), and the allocation of dynamic arrays can be verified.
3. Design of the finite element program The design of an OO program is crucial to its success. First, we determine the type of objects needed, depending on the specific application. Next, we define their interfaces, that is, what these objects need to do and how they will communicate with each other. After these specifications are ready, we can then start implementation. It is a good habit to test various pieces of the program as they are implemented, followed by final testing of the full system. For any application there is no one optimum design of the classes. A good design should exploit the benefits offered by the object-oriented approach as well as addressing issues such as efficiency and ease of implementation. The design of this prototype adaptive FEA program consists of several modules. Each module contains one or more logically related classes.
2.3. Inheritance Inheritance is the ability to reuse and share the data and methods of existing classes to construct hierarchical classes.
† There are eight major classes: Problem, Adaptor, Solver, Domain, Element, Edge, Constraint and Grid_Point. † There are eight auxiliary classes: Skyline, Spr
464
J.E. Akin, M. Singh / Advances in Engineering Software 33 (2002) 461–468
Fig. 3. Illustrating the program design.
(superconvergent patch recovery), Grid_Point_Sublist, Gauss, Material, Dof (degree of freedom), Point, and Solution. † There are six modules containing data and procedures, used by more than one class: the User_Specific, Interpolation, Element_Assembly, Control_Data, Utils and Precision classes. These common data and methods are grouped (encapsulated) on the basis of logical or functional similarities. Fig. 3 depicts the design of this FEA prototype. In this figure the
class from where the arrow starts is either directly contained in or its functions are used by the class that is pointed to by the arrowhead. For example, the highest level class, Problem, is at the top. Thus relationships can be of containment of lower classes, or only using their methods. Note that to remove confusion between a finite element node and a node in the linked (and circular) list data structures the finite element node is referred to as a grid point while a list node is simply referred to as a node. We will now discuss the role of these modules. In this design the information about the finite element
Fig. 4. The Problem, Adaptor, Solver and Domain classes.
J.E. Akin, M. Singh / Advances in Engineering Software 33 (2002) 461–468
465
Fig. 5. The Element, Edge, GridPt and DircBc classes.
model is separated from the details of finding the solution. In other words, data are separated from the analysis method, making the program flexible and easy to debug. The highest level class is the Problem class. Its structure is shown in Fig. 4. The model information is contained in the Domain class and the adaptation and solver functionality are in the Adaptor and Solver classes, respectively. The Adaptor and the Solver classes are shown in Fig. 4. These classes need access to the data and methods encapsulated in the Domain class in order to compute the solution. So, a pointer link to the Domain class is stored in these classes. The solver class contains the finite element stiffness matrix in variable Stiff_Mat, which is an instance of Skyline class. Dof_Arr, an array of type Dof, stores the mapping between the local and the system numbering schemes. Soln, an instance of Solution class, holds the solution values for each degree of freedom. The Domain class contains data structures for the Element, Edge and GridPt (grid point) classes and the iterators for these lists. The iterators are used for implementing efficient traverse and query algorithms discussed in Ref. [13]. The Element module contains the Element class, Element Link List class, Element List Iterator class and Element Pointer class. The Element List class and the Element Iterator class provide methods for link list manipulation and traversal. The Element class is shown in Fig. 5. In adaptive analysis information about the neighboring elements is required to
recover accurate gradients in an element. This information is stored in an array of pointers to the adjacent elements. Information about the edges is likewise stored in another array of pointers. The geometric grid point coordinates of an Element are calculated and stored in Gm_GrdPt_Crds, to avoid expensive re-computation for each iteration. Dof_Map contains the system degree of freedom numbers for the local element degree of freedom numbers. Smp_Data is a link list containing the Gauss sampling Point data for the Element, used for accurate gradient recovery. The Edge module contains the Edge class, Edge list class and the Edge list iterator class. The Edge class is depicted in Fig. 5. Each Edge contains a circular list, GrdPtSublist, which stores pointers to the grid points lying on the Edge. IsBndryFlux is a logical field, used to distinguish Edges with flux boundary condition from rest of the Edges. The Grid_Point module also has a similar structure as the Element and Edge modules; it contains the Grid_Point class, Grid_Point list class and Grid_Point list iterator class. The Grid_Point class is shown in Fig. 5. Pos, an instance of Point class, keeps track of the coordinates of the Grid Point. Dof_info is an integer flag that is used to classify analytic or geometric nature of the Grid Point. Due to the use of different interpolation functions for the geometry and the solution, a Grid Point could be either analytical, geometric or both. The logical function Get_Is_PureGeometric( ) and its companion function Get_Is_PureAnalytic( ) make this
466
J.E. Akin, M. Singh / Advances in Engineering Software 33 (2002) 461–468
nature transparent to the user. Dof_info also serves to uniquely identify each Grid Point. After each adaptation the Dof_info is regenerated. The Constraint module contains the Dirchlet and Boundary Flux classes. The Dirchlet class is shown in Fig. 5. The integer flag Param indicates to which parameter of each grid point the boundary condition is applied. The auxiliary classes also deserve some description. The Skyline class is used to implement the skyline storage of a square matrix. The generic functions like Get(row,col), hide the implementation complexity from the user. This allows us the flexibility to use a different matrix storage method, since by implementing the corresponding Get(row,col) method, the rest of the program will remain unchanged. To recover the superconvergent gradients we need a list to store the sampling point information of the element. Sampling Point Data class and the corresponding list and iterator classes in the Spr module provide this functionality. The class GridPtSublist forms a circular list of Grid Point pointers. This list is used to represent the Grid Points on an Edge. The Gauss class, shown in Fig. 2 stores the position and the weight associated with a Gauss point. The Material class is used to read and store the properties. The Dof (Degree of freedom) class is used to form mapping between system and local degree of freedoms and also provides equation renumbering algorithms. The Point class, which was shown in Fig. 1, implements a point in one-, two- or three-dimensional space. The Solution class is used to hold the system solution. This prototype adaptive FEA program also has some modules that encapsulate common attributes and methods. The Interpolation module provides methods for both geometric and solution interpolation. The User_Specific class holds procedures such as exact solution for benchmark test cases and source functions. The Element_Assembly module provides the structures needed to generate the element stiffness matrix and assemble it into the system matrix. The Control_Data module stores global level tags and flags, while Utils module stores globally used utility procedures for implementing validation and error checking. Precision module defines the precision type that is used in the program.
4. Discussion In this section we will discuss the strengths and weaknesses of using the OO approach and demonstrate that it is worth the effort. Some of the advantages of the following OO policies are † The main advantage of an OO scheme lies in data encapsulation. The information is decentralized and as such it may be stored at easily reachable places and processed at a suitable time. Secondly, since objects possess more than intrinsic identifiers, they can possess
†
†
†
†
†
†
other objects. For example, the Element class does not contain the number of material properties but contains the material properties themselves. In OO programming we group data (attributes) and functions (methods) in a class. We do not have to pass arguments, just by passing the object we can use its methods to find information about its attributes. For example, the element class has functions for constructing and initializing an element, generating the element stiffness matrix and element-level post processing. The program is more organized and intuitive, as objects present themselves to the external world in a more meaningful way (for example, “I am a beam on the third floor connecting A and B” as opposed to “I am array element 5 in the element array”). A finite element application maps very well into the computational domain using the OO design techniques. One part of the code can be debugged and modified without affecting the remainder. This is of great relevance to scientific research as the researcher can try new algorithms, new data structures and different problems with less programming effort. In this application the Solver is uncoupled from the domain information (Element, Edges, etc.). The Domain is also independent of the matrix storage technique and the data structures. These benefits where experienced in our research as we experimented with different element types, different solvers, different data structures and different precision with relative ease without affecting the rest of the program. Since we are separating “what we are doing” from “how we are doing it”, it is easier for other people to understand the code. An OO approach encourages localized, distributed grouping of data and tasks because it is relatively natural to see how data and tasks can be distributed over processors (useful for parallel machines) and among individuals (good for team environment). By making use of the inheritance, we can add new components and reuse previously written code.
Now we will study some of the actual and perceived weaknesses of using OO techniques. † One of the main weaknesses of OO programming is its overhead, which makes it inherently inefficient with respect to speed, primarily because of the addition of new procedure calls for implementing encapsulation and in some cases due to late (run-time) binding typically associated with OO environments. The increase in overhead due to more function calls can in many cases be offset by the reduction in human time and effort. The number of function calls can be reduced by the use of pointers. The encapsulation can be broken for efficiency for critical operations that are performed repeatedly. † Another criticism of OO methodologies is creation of
J.E. Akin, M. Singh / Advances in Engineering Software 33 (2002) 461–468
rigid hierarchies of classes. Too much emphasis on code reuse and the tendency to classify all objects into an existing class hierarchy can lead to deep hierarchies of unrelated classes. This type of framework can break down when an object that has multiple characteristics is encountered. This was also discussed in a recent posting on the Cþ þ users mailing list [3], as the Platypus effect. A platypus is a fur-bearing, egg-laying, aquatic mammal and to classify it in the class hierarchy of species using only ‘is-a’ inheritance leads to breakdown of the usual hierarchy. This is due to overuse of inheritance (‘is-a’) in the design and is not necessarily a weakness of OO method. With proper design, such as use of clear logical distinguishing of classes and judicial use of inheritance (‘is-a’) and composition (‘has-a’), a flexible infrastructure can be achieved. In this program we only used composition. † It should be noted that the OO analysis and OO design stages that must be completed before implementation can begin can be very time consuming and involved [2]. Fortran 90/95 provides good support for OO development, however, it lacks some of the features of other OO languages like Cþ þ and Java. One of the missing features is the lack of intrinsic support for ‘is-a’ inheritance. This type of inheritance can lead to code reuse and is more useful when there are multiple classes of similar kind. It can be implemented in F90/F95 but a simple procedure is not part of the language standard. However, it is expected in the Fortran 2000 standard. This type of inheritance is used successfully in places like graphical user interface (GUI) development, but is not as useful in scientific computing. In the present application the OO infrastructure was constructed only using composition. One place where ‘is-a’ inheritance could be used is for representing the element hierarchy to add multiple types of elements without having to modify part of the code. However, this is achieved in the present code without inheritance, by using parameters to distinguish the various element types. The parameters N_SPACE (dimension of physical space), GRDPT_PER_EL (number of grid points in a element) and NO_PARAMS_PER_GrdPt (number of degrees of freedom per grid point) are used to make the methods of element class general enough to be used by element objects with different geometric shape functions. We also use a function that generates variable degree element shape functions depending upon these parameters. Szymanski’s technique [15] implements ‘is-a’ inheritance in Fortran 90 with a little more coding and extra function calls. Another related feature absent in Fortran 90 is run-time polymorphism. This is a related concept to ‘is-a’ inheritance. It can be useful in some situation though adds expense and makes it difficult for the compiler to optimize. This can also be implemented in F90 using Szymanski’s method. Fortran 90 does not have templates like Cþ þ , though they can again be implemented using a preprocessor
467
[2]. Fortran 90 has many other useful features that are missing in other languages, especially for scientific computing. Different OO languages have their own focus: Cþ þ is a general purpose language, Fortran 90 is more equipped for scientific computation and Java is more adept for developing portable and multi-threaded applications. Today’s engineers need to know more than one programming language. In this research, we were able to reuse large amount of previously written Fortran 77 routines since it remains a subset of F95. To automate some of the repeated tasks of pre- and post-processing, there was a need for a GUI. Therefore, a GUI was developed in Java and was interlinked with the Fortran 90 program.
5. Results The P-adaptive finite element application developed during this research has shown good results. Benchmark problems were run to verify the method. They provided further insight into the adaptive refinement process and were used to test some variations to the original superconvergent patch method. Here we show the results of one of the benchmark cases: Poisson’s equation for potential flow around a cylinder, which is a two-dimensional problem of an ideal incompressible fluid. The domain was discretized into an initial Q8 16 £ 24 mesh (8-noded elements), with combinations of natural and essential boundary conditions. Fig. 6 shows results in terms of rate of convergence and the number of unknowns in the problem. The unknowns increase as the polynomial degree of independent edges are increased by adding additional nodes in regions selected by the error estimator. The line with p ticks is the standard energy norm error estimator and superconvergent patch recovery process [16 –18], while the line with K ticks shows Blacker’s variation of adding the residual error in the partial differential equation to the error estimator [4]. Both curves involve adaptive meshes containing both even and odd degree polynomial elements. Another test was performed to determine the effect of using only even degree polynomials versus odd or mixed degree polynomials in the element interpolation functions. Zienkiewicz and Zhu [18] suggested that even degree polynomials give a higher rate of convergence for their gradients. Our results agree with their observation. The line with þ ticks in Fig. 6 show the convergence rate for the standard SPR where the adapting edges were restricted to even degree polynomials, while the line with the p ticks used odd degree polynomials. Likewise, the A symbol denotes even polynomial solutions of SPR plus equilibrium residuals. The even degree solution is cheaper. Early stagnation in the even degree is due to the reduction in the number of degrees of freedom as a result of the forced shift from odd to lower order even degrees.
468
J.E. Akin, M. Singh / Advances in Engineering Software 33 (2002) 461–468
Fig. 6. Convergence comparisons for an adaptive potential flow problem.
6. Conclusion Much has been said about object orientation making tasks easy for the software developer. However, OO programming does not obviate the need for forward thinking nor does it make it impossible to introduce a bug. The success of an OO program depends heavily on its design. Since the OO program is supposed to be flexible and easily extensible, the developer needs to think beyond the present application. Initially the program design can consume a large amount of total application development time, as was the case with this adaptive FEA prototype. As the developer gains experience in OO principles, subsequent projects will require less time for analysis and design. The language used for the OO development depends upon the type of application. Fortran 90/95 is well equipped to provide OO software for scientific computing and it usually executes much faster than Cþ þ [11].
References [1] Akin JE. Finite elements for analysis and design. London: Academic Press; 1995. [2] Akin JE. Object oriented programming via Fortran 90/95. Cambridge, Cambridge University Press, 2002. [3] Anonymous. Encapsulation, inheritance and the platypus effect. The Cþ þ developer discussion diary, www.advogato.org/article/83.html. [4] Blacker T, Belytschko T. Superconvergent patch recovery with equilibrium and conjoint interpolant enhancements. Int J Numer Meth Engng 1994;37:517–36. [5] Dubois-P‘elerin Y, Pegon P. Improving modularity in object-oriented finite element programming. Commun Numer Meth Engng 1997;13: 193–8.
[6] Dubois-P‘elerin Y, Zimmermann T, Bomme P. Object-oriented finite element programming. II. A prototype program in Smalltalk. Comput Meth Appl Mech Engng 1992;98:361 –97. [7] Dubois-P‘elerin Y, Zimmermann T. Object-oriented finite element programming. III. An efficient implementation in Cþ þ. Comput Meth Appl Mech Engng 1993;108:165 –83. [8] Forde BWR, Foschi RB, Stiemer SF. Object-oriented finite element analysis. Comput Struct 1990;34:355 –74. [9] Liu J-L, Lin I-J, Shih M-Z, Chen R-C, Hsieh M-C. Object-oriented programming of adaptive finite element and finite volume methods. Appl Numer Math 1996;21:439–67. [10] Gang YG. Object oriented models for numerical and finite element analysis. PhD Thesis. The Ohio State University, USA; 1994. [11] Norton CD. Object-oriented programming paradigms in scientific computing. PhD Thesis. Rensselaer Polytechnic Institute, USA; 1996. [12] Rehak DR, Baugh Jr JW. Alternative programming techniques for finite element program development. Proceedings of the IABSE Colloquium on Expert Systems in Civil Engineering, Bergamo, Italy; 1989. [13] Singh M. Object-oriented implementation of P-adaptive finite element method. MS Thesis. Mechanical Engineering and Material Science, Rice University, Houston, TX, USA; 1999. [14] Szymanski BK, Decyk VK, Norton CD. Expressing object-oriented concepts in Fortran 90. ACM Fortran Forum 1997;16:13– 18. [15] Szymanski BK, Decyk VK, Norton CD. How to support inheritance and run-time polymorphism in Fortran 90. Comput Phys Commun 1998;115:9–17. [16] Zienkiewicz OC, Zhu JZ. A simple error estimator and adaptive procedure for practical engineering analysis. Int J Numer Meth Engng 1987;24:337–57. [17] Zienkiewicz OC, Zhu JZ. The superconvergent patch recovery (SPR) and adaptive finite element refinement. Comput Meth Appl Mech Engng 1992;101:207 –24. [18] Zienkiewicz OC, Zhu JZ. The superconvergent patch recovery and a posteriori error estimates. Part I, II. Int J Numer Meth Engng 1992;33: 1331–82. [19] Zimmermann T, Dubois-P‘elerin Y, Bomme P. Object-oriented finite element programming. I. Governing principles. Comput Meth Appl Mech Engng 1992;98:291 –303.
Object-oriented Fortran 90 P-adaptive finite element method J.E. Akina,*, M. Singhb a
Rice University, Houston, TX 77005, USA b COADE Inc., Houston, TX 77070, USA
Abstract Over the last decade, there has been an increased awareness of the benefits of employing Object-Oriented (OO) design and methodologies for development of software. Among the various languages available for OO development, Fortran 95 has some clear advantages for scientific and engineering programming. It offers features similar to other OO languages like Cþ þ and Smalltalk as well as extensive and efficient numerical abilities. This paper will describe the OO design and implementation of P-adaptive finite element analysis (FEA) using Fortran. We will demonstrate how various OO principles were successfully employed to achieve greater flexibility, easier maintainability and extensibility. This is helpful for a complex program like an adaptive finite element implementation. q 2002 Civil-Comp Ltd and Elsevier Science Ltd. All rights reserved. Keywords: Object-oriented design; P-adaptive finite element analysis; Fortran
1. Introduction The idea of object-oriented programming can be traced back to Simula 67 language. It was realized that OO techniques could solve the problems of code reusability and maintainability encountered in engineering software development. Consequently, there has been considerable interest over the last decade in this alternate philosophy of programming. We can view OO programming as an analog model of real world phenomenon. A program thus builds an environment containing objects (such as beams, floors) with attributes analogous to real world counterparts (such as weight, lengths and properties) and analogous behavior (deformation under loads). An object can also be an abstract entity (such as matrix, vector, etc.). Akin [2], Szymanski [14,15] and his student Norton [11] have elaborated on the implementation of the OO concepts in F90. Rehak et al. [12] and Forde et al. [8] first applied OO principles to finite element programming. Zimmermann et al. [6,7,19] in a series of papers presented two OO finite element toolkits, implemented in Smalltalk and Cþ þ . Dubois-P‘elerin et al. [5] proposed further modularity by delineating two kinds of behavior of a finite element software: analysis type (solvers used, static or dynamic analysis) and domain information (element, grid points, boundary conditions). Liu et al. [9] introduced an OO * Corresponding author. E-mail address: [email protected] (J.E. Akin).
implementation of adaptive finite element and finite volume methods. They used the Cþ þ language to implement the program, which performed h-refinement and offered a choice of direct or iterative solvers. While much of the previous OO work is in Cþ þ , Smalltalk and Pascal, no adaptive finite element analysis (FEA) system has been reported in Fortran 90 (F90). Fortran 90 and later 95 introduced dramatic new features over Fortran 77. These features helped in flexible program design, dynamic memory management and powerful built in numerical capabilities. The Type command, used for creation of user-defined types in the Module construct allows easier implementation of the OO style of programming. We will illustrate various concepts of OO programming in an FEA flavor with examples from the actual implementation. Since its birth in the 1940s, FEA has matured and is now employed in a wide range of problems, from solid mechanics to heat transfer, fluid mechanics, acoustics, and electromagnetism. The main advantage of FEA is its ability to model complex and arbitrary shaped domains. It also supports modeling of general boundary conditions and non-homogeneous anisotropic materials. To achieve a greater degree of confidence in an FEA, adaptive techniques have been developed. At the core of an adaptive FEA program is the error-estimator, which provides numerical estimates of the error in the solution. Using this error information the mesh is then accordingly modified. This cycle is repeated until an error less than
0965-9978/02/$ - see front matter q 2002 Civil-Comp Ltd and Elsevier Science Ltd. All rights reserved. PII: S 0 9 6 5 - 9 9 7 8 ( 0 2 ) 0 0 0 4 8 - 0
462
J.E. Akin, M. Singh / Advances in Engineering Software 33 (2002) 461–468
Fig. 1. A Point class demonstrating various object-oriented concepts.
the specified threshold is achieved. There are three basic approaches for mesh adaptation: H-refinement, which subdivides the selected elements into smaller elements; P-refinement, which reduces the error in the solution by increasing the polynomial order, p; of the interpolating functions; and R-refinement, which relocates the position of nodes in a finite element mesh. This paper will describe the implementation of the P-adaptive FEA. More details are available in Refs. [1,13]. An adaptive FEA program places high demands on software development for effectively managing the complexities resulting from error-estimation and mesh
adaptation which requires modifying the underlining data structures. Traditionally, FEA programs have been programmed procedurally in Fortran 77 [1]. In a procedural style, also termed as water-fall technique [10], a program is viewed as a set of algorithms. Therefore, the program depends closely on the algorithm being implemented, creating a rigid framework where a change becomes costly and sometimes unmanageable unless the algorithm components are highly modular. OO programming is capable of overcoming some of the limitations of procedural programming, but does so at additional cost. We will demonstrate this in subsequent sections.
J.E. Akin, M. Singh / Advances in Engineering Software 33 (2002) 461–468
463
Fortran 90 supports the following types of inheritance:
Fig. 2. Example of inheritance by Composition.
2. Main OO features in F90 2.1. Class Fortran 90 allows the creation of user-defined types to encapsulate the related concepts into a single form. The type thus created is associated with the related methods within a module block, which provides an interface for communication with external objects. This is similar to the notion of Class in Cþ þ and Java. A module can contain more than one logically related class. Fig. 1 shows a Point class within a module. This class encapsulates the coordinates and spatial dimensions of the point. The USE statement provides higher class objects in different modules with access to the methods of lower class objects (this is not merely a text substitution, as done by the INCLUDE command). This can be seen in the case of the Point class located in Point Module, as the Point Class gets access to the data and methods in the Dbl_Precision_Module. We can declare the contents in a module as private or public to provide concepts the outside user access only to relevant information and to hide the complexity within the module. 2.2. Polymorphism It is the property of OO languages to allow the use of identical function names for logical similar methods in different types. This property can be seen in some intrinsic functions like int(a ), where a can be of integer, real or complex type. This can be extended to user-defined objects using generic programming, operator overloading and the interface construct. For example, in Fig. 1 Init_Point is bound to Init using the interface construct. This binding can be performed in all the classes in the program so that the user only utilizes one call to Init instead of individual calls to Init_Element, Init_Point, etc. At compile time, the compiler uses the information about the argument types to uniquely determine which procedure to execute. Similarly, operators can be overloaded for user-defined types (already existing for implicit types), and additional operators can be defined. Fig. 1, shows how the ‘ ¼ ¼ ’, equality operator is defined for the Point class. The assignment operator ‘ ¼ ’ between two objects of the same type is implicitly implemented by the language.
† Inheritance by granting access. This is the simplest form of inheritance. The use statement in Fortran 90 allows objects in one module to access the public methods and data in the base module. † Inheritance by composition. Also known as ‘has-a’ inheritance, this type of inheritance uses other classes to build larger and more complex classes. As shown in Fig. 2, the Gauss_Type class contains an instance of Point class to define its position in the space. The Gauss_Type class is granted access to contents of the Point class’s public declarations and USE statements. † Inheritance by subtyping. This type of inheritance can be defined as ‘is-a’ inheritance. The ‘is-a’ inheritance is explicitly implemented in Cþ þ and Java languages, while in Fortran 90 it has to be simulated by the programmer using composition. Szymanski et al. [15] introduced the technique to implement this type of inheritance in Fortran 90. It can add dynamic capabilities to the program, but at the added expense of a small increase in the complexity of the program. The array syntax and intrinsic array operators in F90 provide advanced numerical capabilities very similar to those in the Matlab environment. Using the assumed shape array, automatic array, and the SIZE command, the procedures can be made independent of the length of the array arguments, lending more flexibility to the code. Fortran 90 also provides some safety features; for example, pointers can be checked to verify that they are associated (pointing to something valid), and the allocation of dynamic arrays can be verified.
3. Design of the finite element program The design of an OO program is crucial to its success. First, we determine the type of objects needed, depending on the specific application. Next, we define their interfaces, that is, what these objects need to do and how they will communicate with each other. After these specifications are ready, we can then start implementation. It is a good habit to test various pieces of the program as they are implemented, followed by final testing of the full system. For any application there is no one optimum design of the classes. A good design should exploit the benefits offered by the object-oriented approach as well as addressing issues such as efficiency and ease of implementation. The design of this prototype adaptive FEA program consists of several modules. Each module contains one or more logically related classes.
2.3. Inheritance Inheritance is the ability to reuse and share the data and methods of existing classes to construct hierarchical classes.
† There are eight major classes: Problem, Adaptor, Solver, Domain, Element, Edge, Constraint and Grid_Point. † There are eight auxiliary classes: Skyline, Spr
464
J.E. Akin, M. Singh / Advances in Engineering Software 33 (2002) 461–468
Fig. 3. Illustrating the program design.
(superconvergent patch recovery), Grid_Point_Sublist, Gauss, Material, Dof (degree of freedom), Point, and Solution. † There are six modules containing data and procedures, used by more than one class: the User_Specific, Interpolation, Element_Assembly, Control_Data, Utils and Precision classes. These common data and methods are grouped (encapsulated) on the basis of logical or functional similarities. Fig. 3 depicts the design of this FEA prototype. In this figure the
class from where the arrow starts is either directly contained in or its functions are used by the class that is pointed to by the arrowhead. For example, the highest level class, Problem, is at the top. Thus relationships can be of containment of lower classes, or only using their methods. Note that to remove confusion between a finite element node and a node in the linked (and circular) list data structures the finite element node is referred to as a grid point while a list node is simply referred to as a node. We will now discuss the role of these modules. In this design the information about the finite element
Fig. 4. The Problem, Adaptor, Solver and Domain classes.
J.E. Akin, M. Singh / Advances in Engineering Software 33 (2002) 461–468
465
Fig. 5. The Element, Edge, GridPt and DircBc classes.
model is separated from the details of finding the solution. In other words, data are separated from the analysis method, making the program flexible and easy to debug. The highest level class is the Problem class. Its structure is shown in Fig. 4. The model information is contained in the Domain class and the adaptation and solver functionality are in the Adaptor and Solver classes, respectively. The Adaptor and the Solver classes are shown in Fig. 4. These classes need access to the data and methods encapsulated in the Domain class in order to compute the solution. So, a pointer link to the Domain class is stored in these classes. The solver class contains the finite element stiffness matrix in variable Stiff_Mat, which is an instance of Skyline class. Dof_Arr, an array of type Dof, stores the mapping between the local and the system numbering schemes. Soln, an instance of Solution class, holds the solution values for each degree of freedom. The Domain class contains data structures for the Element, Edge and GridPt (grid point) classes and the iterators for these lists. The iterators are used for implementing efficient traverse and query algorithms discussed in Ref. [13]. The Element module contains the Element class, Element Link List class, Element List Iterator class and Element Pointer class. The Element List class and the Element Iterator class provide methods for link list manipulation and traversal. The Element class is shown in Fig. 5. In adaptive analysis information about the neighboring elements is required to
recover accurate gradients in an element. This information is stored in an array of pointers to the adjacent elements. Information about the edges is likewise stored in another array of pointers. The geometric grid point coordinates of an Element are calculated and stored in Gm_GrdPt_Crds, to avoid expensive re-computation for each iteration. Dof_Map contains the system degree of freedom numbers for the local element degree of freedom numbers. Smp_Data is a link list containing the Gauss sampling Point data for the Element, used for accurate gradient recovery. The Edge module contains the Edge class, Edge list class and the Edge list iterator class. The Edge class is depicted in Fig. 5. Each Edge contains a circular list, GrdPtSublist, which stores pointers to the grid points lying on the Edge. IsBndryFlux is a logical field, used to distinguish Edges with flux boundary condition from rest of the Edges. The Grid_Point module also has a similar structure as the Element and Edge modules; it contains the Grid_Point class, Grid_Point list class and Grid_Point list iterator class. The Grid_Point class is shown in Fig. 5. Pos, an instance of Point class, keeps track of the coordinates of the Grid Point. Dof_info is an integer flag that is used to classify analytic or geometric nature of the Grid Point. Due to the use of different interpolation functions for the geometry and the solution, a Grid Point could be either analytical, geometric or both. The logical function Get_Is_PureGeometric( ) and its companion function Get_Is_PureAnalytic( ) make this
466
J.E. Akin, M. Singh / Advances in Engineering Software 33 (2002) 461–468
nature transparent to the user. Dof_info also serves to uniquely identify each Grid Point. After each adaptation the Dof_info is regenerated. The Constraint module contains the Dirchlet and Boundary Flux classes. The Dirchlet class is shown in Fig. 5. The integer flag Param indicates to which parameter of each grid point the boundary condition is applied. The auxiliary classes also deserve some description. The Skyline class is used to implement the skyline storage of a square matrix. The generic functions like Get(row,col), hide the implementation complexity from the user. This allows us the flexibility to use a different matrix storage method, since by implementing the corresponding Get(row,col) method, the rest of the program will remain unchanged. To recover the superconvergent gradients we need a list to store the sampling point information of the element. Sampling Point Data class and the corresponding list and iterator classes in the Spr module provide this functionality. The class GridPtSublist forms a circular list of Grid Point pointers. This list is used to represent the Grid Points on an Edge. The Gauss class, shown in Fig. 2 stores the position and the weight associated with a Gauss point. The Material class is used to read and store the properties. The Dof (Degree of freedom) class is used to form mapping between system and local degree of freedoms and also provides equation renumbering algorithms. The Point class, which was shown in Fig. 1, implements a point in one-, two- or three-dimensional space. The Solution class is used to hold the system solution. This prototype adaptive FEA program also has some modules that encapsulate common attributes and methods. The Interpolation module provides methods for both geometric and solution interpolation. The User_Specific class holds procedures such as exact solution for benchmark test cases and source functions. The Element_Assembly module provides the structures needed to generate the element stiffness matrix and assemble it into the system matrix. The Control_Data module stores global level tags and flags, while Utils module stores globally used utility procedures for implementing validation and error checking. Precision module defines the precision type that is used in the program.
4. Discussion In this section we will discuss the strengths and weaknesses of using the OO approach and demonstrate that it is worth the effort. Some of the advantages of the following OO policies are † The main advantage of an OO scheme lies in data encapsulation. The information is decentralized and as such it may be stored at easily reachable places and processed at a suitable time. Secondly, since objects possess more than intrinsic identifiers, they can possess
†
†
†
†
†
†
other objects. For example, the Element class does not contain the number of material properties but contains the material properties themselves. In OO programming we group data (attributes) and functions (methods) in a class. We do not have to pass arguments, just by passing the object we can use its methods to find information about its attributes. For example, the element class has functions for constructing and initializing an element, generating the element stiffness matrix and element-level post processing. The program is more organized and intuitive, as objects present themselves to the external world in a more meaningful way (for example, “I am a beam on the third floor connecting A and B” as opposed to “I am array element 5 in the element array”). A finite element application maps very well into the computational domain using the OO design techniques. One part of the code can be debugged and modified without affecting the remainder. This is of great relevance to scientific research as the researcher can try new algorithms, new data structures and different problems with less programming effort. In this application the Solver is uncoupled from the domain information (Element, Edges, etc.). The Domain is also independent of the matrix storage technique and the data structures. These benefits where experienced in our research as we experimented with different element types, different solvers, different data structures and different precision with relative ease without affecting the rest of the program. Since we are separating “what we are doing” from “how we are doing it”, it is easier for other people to understand the code. An OO approach encourages localized, distributed grouping of data and tasks because it is relatively natural to see how data and tasks can be distributed over processors (useful for parallel machines) and among individuals (good for team environment). By making use of the inheritance, we can add new components and reuse previously written code.
Now we will study some of the actual and perceived weaknesses of using OO techniques. † One of the main weaknesses of OO programming is its overhead, which makes it inherently inefficient with respect to speed, primarily because of the addition of new procedure calls for implementing encapsulation and in some cases due to late (run-time) binding typically associated with OO environments. The increase in overhead due to more function calls can in many cases be offset by the reduction in human time and effort. The number of function calls can be reduced by the use of pointers. The encapsulation can be broken for efficiency for critical operations that are performed repeatedly. † Another criticism of OO methodologies is creation of
J.E. Akin, M. Singh / Advances in Engineering Software 33 (2002) 461–468
rigid hierarchies of classes. Too much emphasis on code reuse and the tendency to classify all objects into an existing class hierarchy can lead to deep hierarchies of unrelated classes. This type of framework can break down when an object that has multiple characteristics is encountered. This was also discussed in a recent posting on the Cþ þ users mailing list [3], as the Platypus effect. A platypus is a fur-bearing, egg-laying, aquatic mammal and to classify it in the class hierarchy of species using only ‘is-a’ inheritance leads to breakdown of the usual hierarchy. This is due to overuse of inheritance (‘is-a’) in the design and is not necessarily a weakness of OO method. With proper design, such as use of clear logical distinguishing of classes and judicial use of inheritance (‘is-a’) and composition (‘has-a’), a flexible infrastructure can be achieved. In this program we only used composition. † It should be noted that the OO analysis and OO design stages that must be completed before implementation can begin can be very time consuming and involved [2]. Fortran 90/95 provides good support for OO development, however, it lacks some of the features of other OO languages like Cþ þ and Java. One of the missing features is the lack of intrinsic support for ‘is-a’ inheritance. This type of inheritance can lead to code reuse and is more useful when there are multiple classes of similar kind. It can be implemented in F90/F95 but a simple procedure is not part of the language standard. However, it is expected in the Fortran 2000 standard. This type of inheritance is used successfully in places like graphical user interface (GUI) development, but is not as useful in scientific computing. In the present application the OO infrastructure was constructed only using composition. One place where ‘is-a’ inheritance could be used is for representing the element hierarchy to add multiple types of elements without having to modify part of the code. However, this is achieved in the present code without inheritance, by using parameters to distinguish the various element types. The parameters N_SPACE (dimension of physical space), GRDPT_PER_EL (number of grid points in a element) and NO_PARAMS_PER_GrdPt (number of degrees of freedom per grid point) are used to make the methods of element class general enough to be used by element objects with different geometric shape functions. We also use a function that generates variable degree element shape functions depending upon these parameters. Szymanski’s technique [15] implements ‘is-a’ inheritance in Fortran 90 with a little more coding and extra function calls. Another related feature absent in Fortran 90 is run-time polymorphism. This is a related concept to ‘is-a’ inheritance. It can be useful in some situation though adds expense and makes it difficult for the compiler to optimize. This can also be implemented in F90 using Szymanski’s method. Fortran 90 does not have templates like Cþ þ , though they can again be implemented using a preprocessor
467
[2]. Fortran 90 has many other useful features that are missing in other languages, especially for scientific computing. Different OO languages have their own focus: Cþ þ is a general purpose language, Fortran 90 is more equipped for scientific computation and Java is more adept for developing portable and multi-threaded applications. Today’s engineers need to know more than one programming language. In this research, we were able to reuse large amount of previously written Fortran 77 routines since it remains a subset of F95. To automate some of the repeated tasks of pre- and post-processing, there was a need for a GUI. Therefore, a GUI was developed in Java and was interlinked with the Fortran 90 program.
5. Results The P-adaptive finite element application developed during this research has shown good results. Benchmark problems were run to verify the method. They provided further insight into the adaptive refinement process and were used to test some variations to the original superconvergent patch method. Here we show the results of one of the benchmark cases: Poisson’s equation for potential flow around a cylinder, which is a two-dimensional problem of an ideal incompressible fluid. The domain was discretized into an initial Q8 16 £ 24 mesh (8-noded elements), with combinations of natural and essential boundary conditions. Fig. 6 shows results in terms of rate of convergence and the number of unknowns in the problem. The unknowns increase as the polynomial degree of independent edges are increased by adding additional nodes in regions selected by the error estimator. The line with p ticks is the standard energy norm error estimator and superconvergent patch recovery process [16 –18], while the line with K ticks shows Blacker’s variation of adding the residual error in the partial differential equation to the error estimator [4]. Both curves involve adaptive meshes containing both even and odd degree polynomial elements. Another test was performed to determine the effect of using only even degree polynomials versus odd or mixed degree polynomials in the element interpolation functions. Zienkiewicz and Zhu [18] suggested that even degree polynomials give a higher rate of convergence for their gradients. Our results agree with their observation. The line with þ ticks in Fig. 6 show the convergence rate for the standard SPR where the adapting edges were restricted to even degree polynomials, while the line with the p ticks used odd degree polynomials. Likewise, the A symbol denotes even polynomial solutions of SPR plus equilibrium residuals. The even degree solution is cheaper. Early stagnation in the even degree is due to the reduction in the number of degrees of freedom as a result of the forced shift from odd to lower order even degrees.
468
J.E. Akin, M. Singh / Advances in Engineering Software 33 (2002) 461–468
Fig. 6. Convergence comparisons for an adaptive potential flow problem.
6. Conclusion Much has been said about object orientation making tasks easy for the software developer. However, OO programming does not obviate the need for forward thinking nor does it make it impossible to introduce a bug. The success of an OO program depends heavily on its design. Since the OO program is supposed to be flexible and easily extensible, the developer needs to think beyond the present application. Initially the program design can consume a large amount of total application development time, as was the case with this adaptive FEA prototype. As the developer gains experience in OO principles, subsequent projects will require less time for analysis and design. The language used for the OO development depends upon the type of application. Fortran 90/95 is well equipped to provide OO software for scientific computing and it usually executes much faster than Cþ þ [11].
References [1] Akin JE. Finite elements for analysis and design. London: Academic Press; 1995. [2] Akin JE. Object oriented programming via Fortran 90/95. Cambridge, Cambridge University Press, 2002. [3] Anonymous. Encapsulation, inheritance and the platypus effect. The Cþ þ developer discussion diary, www.advogato.org/article/83.html. [4] Blacker T, Belytschko T. Superconvergent patch recovery with equilibrium and conjoint interpolant enhancements. Int J Numer Meth Engng 1994;37:517–36. [5] Dubois-P‘elerin Y, Pegon P. Improving modularity in object-oriented finite element programming. Commun Numer Meth Engng 1997;13: 193–8.
[6] Dubois-P‘elerin Y, Zimmermann T, Bomme P. Object-oriented finite element programming. II. A prototype program in Smalltalk. Comput Meth Appl Mech Engng 1992;98:361 –97. [7] Dubois-P‘elerin Y, Zimmermann T. Object-oriented finite element programming. III. An efficient implementation in Cþ þ. Comput Meth Appl Mech Engng 1993;108:165 –83. [8] Forde BWR, Foschi RB, Stiemer SF. Object-oriented finite element analysis. Comput Struct 1990;34:355 –74. [9] Liu J-L, Lin I-J, Shih M-Z, Chen R-C, Hsieh M-C. Object-oriented programming of adaptive finite element and finite volume methods. Appl Numer Math 1996;21:439–67. [10] Gang YG. Object oriented models for numerical and finite element analysis. PhD Thesis. The Ohio State University, USA; 1994. [11] Norton CD. Object-oriented programming paradigms in scientific computing. PhD Thesis. Rensselaer Polytechnic Institute, USA; 1996. [12] Rehak DR, Baugh Jr JW. Alternative programming techniques for finite element program development. Proceedings of the IABSE Colloquium on Expert Systems in Civil Engineering, Bergamo, Italy; 1989. [13] Singh M. Object-oriented implementation of P-adaptive finite element method. MS Thesis. Mechanical Engineering and Material Science, Rice University, Houston, TX, USA; 1999. [14] Szymanski BK, Decyk VK, Norton CD. Expressing object-oriented concepts in Fortran 90. ACM Fortran Forum 1997;16:13– 18. [15] Szymanski BK, Decyk VK, Norton CD. How to support inheritance and run-time polymorphism in Fortran 90. Comput Phys Commun 1998;115:9–17. [16] Zienkiewicz OC, Zhu JZ. A simple error estimator and adaptive procedure for practical engineering analysis. Int J Numer Meth Engng 1987;24:337–57. [17] Zienkiewicz OC, Zhu JZ. The superconvergent patch recovery (SPR) and adaptive finite element refinement. Comput Meth Appl Mech Engng 1992;101:207 –24. [18] Zienkiewicz OC, Zhu JZ. The superconvergent patch recovery and a posteriori error estimates. Part I, II. Int J Numer Meth Engng 1992;33: 1331–82. [19] Zimmermann T, Dubois-P‘elerin Y, Bomme P. Object-oriented finite element programming. I. Governing principles. Comput Meth Appl Mech Engng 1992;98:291 –303.