 |
NeoPZ
|
|
NeoPZ
|
NeoPZ integrates zero, one, two and three dimensional simulations into a single finite element library. It also incorporates non linear geometric maps, hp adaptive meshes and runs a large variety of finite element simulations. It should therefore come as no surprise that its structure is somewhat different from textbook finite element structures.
In this section we describe which finite element concepts were modified or extended in the NeoPZ environment and how these concepts translated in an object oriented framework
Within the geometric mesh, all geometric elements keep track of their neighbours along all the sides (see Topological Concepts associated with an Element) of the element. This information is usually not encountered in other finite element programs.
Within NeoPZ a neighbour is represented by objects of type TPZGeoElSide. An object of type TPZGeoElSide is almost a geometrical element:
The concept of Jacobian matrix in regular finite elements is the (n n) matrix which maps a parameter in the master element to a point in (R^3). In NeoPZ all elements are in (R^3). As such a line element is mapped to a segment in three dimensional space, etc. Therefore the gradient of the map can not be represented in a traditional way. The Jacobian matrix is computed as the QR decomposition of the gradient of ( x()). The upper triangular square matrix /(R/) is called the Jacobian matrix and the orthogonal matrix (Q) is called the axes througout the software. This has evoluated this way because at an earlier stage of the NeoPZ development it seemed to P. Devloo he had developed a new approach to computing the Jacobian matrix.
The (QR) decomposition is computed using a Gram Schmidt orthogonalization. It is worth mentioning that this approach makes the NeoPZ library ready to describe geometric maps in any space dimension.
Within NeoPZ an element is considered as the union of open sets of points. These sets of points are named sides. As such:
All geometries are grouped in the namespace pzgeom. The topology themselves are defined in the namespace pztopology.
Each topology is associated with an area within the dimension associated with the topology. For example the one dimensional line element is associated with the line segment ![$]-1,1[\subset R$](form_8.png)
. A quadrilateral element is associated with the area ![$]-1,1[\times]-1,1[\subset R^2$](form_9.png)
. The area associated with a topology is named parameter space. In finite element textbooks the parameter space is associated with the space of the master element. Theoretically each finite element code can define its own parameter space. In the NeoPZ environment the parameter space is defined and/or can be modified by specifying other topologies.
Each sides of an element associated with a topology (point, line, quadrilateral, etc). The closure of a side (remember that a side is an open set of points) includes its neighbouring topologies. For instance the closure of the line includes two point topologies, the closure of a quadrilateral topology includes the four lines and four points.
The topology associated with a side of a topology is returned in the method Type(int side). This method exists in all classes of the pztopology namespace
The sides included in the closure of a given side are returned in the method LowerDimensionSides.
As each side has its own parameter space, an affine parameter transformation can be defined between the lower dimension sides and the side itself. This affine transformation is returned in the SideToSideTransform method
A more intuitive approach to object oriented finite element programming is to associate a distinct class with each element topology. This was the first approach to programming the NeoPZ library. After incorporating the three dimensional elements it became apparent that very similar code was replicated in the seven traditional topologies. This code replication generated a huge overhead in code maintenance : improvement and/or corrections for one topology had be duplicated to all topologies. These replicated code remained frequently untested.
Migrating the code to a class template concept was a maior improvement in terms of code maintenance. The downside is that template oriented code is much more difficult to read and more distant from traditional (i.e. fortran) coding concepts.
The template oriented implementation of geometric elements allows us to incorporate nonstandard maps as a geometric element specialization. As such we have implemented maps from a line to circel segment (in 3D space), a quadratic map, a wavyline map and a map to a NACA profile. All nonstandard maps are put in the Special Maps directory.
Finite element programming is essentially matrix analysis. Matrix objects are an essential ingredient of finite element programming. Within the NeoPZ environment a matrix object is a linear transformation capable of mapping a vector (matrix) to a vector (matrix). Some types of matrices implement decomposition procedures but the availability of such procedure is not essential. The storage pattern of each matrix class should not be of concern to the user.
Matrices can be of any type defined by a template parameter. All matrix classes are instantiated with float, double, long double and the corresponding complex types.
There is one matrix class which is storage declared : TPZFMatrix<T> implements a full matrix and stores its elements in a column-maior order (http://en.wikipedia.org/wiki/Row-major_order ). This allows the other matrix classes to implement matrix vector multiplications efficiently by using pointer arithmetic.
One of the characteristics of Finite element research is the quest of the more efficient procedure for inverting the global system of equations. Inumerous possibilities present themselves: direct inversion (LU, Cholesky or LDLt), iterative inversion using a Krylov method (conjugate gradient, GMRES, CGStab, Bi-CGStab), preconditioning based on approximate inverses, multigrid etc. It is impractical to provide the user of NeoPZ with a comprehensive list of inversion procedures.
Instead, a TPZSolver<T> class (where T is the data type e.g. float, double, etc) was idealized which represents a transformation process applied to a matrix object and a right hand side. Its derived class TPZMatrixSolver<T> stands for the same concept, but has a matrix object associated with it. TPZStepSolver<T> presents an interface allowing to choose between a direct solver a regular iterative solver (e.g. Jacobi, SOR, SSOR) or de preconditioned Krylov solver. The preconditioner is represented by a Solver<T> object.
Most known solution procedures can and have been implemented using this class structure.
Ever since his first developments in adaptivity, one of the leading authors of NeoPZ has used the concept of /em hangin /em nodes. The name hanging node refers to the fact that the element mesh is nonconforming and that it can be seen that a node /em hangs on the side between two other nodes (in 2D). The concept was generalized to hp-meshes by the same author in 1987 and extended to three dimensions in the same year. Withing NeoPZ, hanging nodes are referred to as /em coefficient /em restraints because this is how they are effectively implemented. Moreover, shape functions are generally not associated with nodes.
Coefficient restraints convey the idea that, in order to obtain conforming approximation space (i.e. continuous), the multiplying coeficients of certain shape functions can not be choosen independently: they are linearly dependent on the value of the multiplying coefficients of other shape functions. This then introduces the concept of dependent and independent connects.
The values which define the depdendency of the coefficient is stored in the TPZConnect objects. The TPZConnect class implements a method HasDependency() that verifies whether the TPZConnect object has a datastructure which defines the relationship between its coefficients and the coefficients of other TPZConnect objects. It is conceivable that a TPZConnect object depends on a TPZConnect object that also has dependency.
The TPZConnect objects which are dependent on other objects have no correspondence in the global system of equations. There value is updated in the TPZCompMesh::LoadSolution method.
The TPZConnect class represents multiplier coefficients associated with a set of shape functions. When two elements share a vertex, the continuity of the solution is obtained by establishing that the multiplying coeficients associated with all elements which share the node are identical.
Withing PZ, H1 shapefunctions are associated with the sides of the elements (see Neigbouring Information). As a consequence a TPZConnect object is created for each side of H1 elements. In a discontinuous Galerkin approximation a unique TPZConnect object is associated with each element. In this case no other element will be associated with this connect leading a discontinuous approximation.
The TPZConnect object contains information related to the numbering of the global system of equations, the maximum order of the shape functions associated with the connect, the number of state variables associated with each connect, the number of shape functions associated with the connect, whether the connect is restrained, whether the connect is associated with a Lagrange multiplier and whether the connect has been condensed at the element level.
The structure of the global system of equations generated by the finite element approximation is determined by the global equation number associated with each shape function. (The set of shape functions determine the approximation space). By modifying the global equation number associated with the shape function, the finite element approximation remains the same. The computational effort needed to invert the global system of equations, on the other hand, is strongly dependent on the ordering of the shape functions. Within PZ, equations are renumbered by the TPZRenumbering class. The TPZRenumbering class relies on the Boost library (http://www.boost.org ) to optimize the bandwidth and/or fillin of the global system of equations.
Rather than associating a global equation number with each shape function, the global equation number are associated with a group of shape functions (which is a TPZConnect object). This is why each connect contains a sequence number as the fSequenceNumber variable.
If the TPZConnect object isn't associated with any element, its sequence number will be set to -1. A TPZConnect object with sequence number equal to -1 is called a "free" connect.
Within PZ the maximum polynomial order of the shapefunctions associated with each side of the elements can be choosen independently. The data structure which represents the polynomial order of approximation is the fOrder variable of the TPZConnect object. This variable is of type unsigned char. This means that the order can vary from 0 to 255.
In finite element approximations of systems of partial differential equations it is customary to associate the same shape function with each state variable. For instance, in two dimensional elasticity each shape function is associated with the horizontal and vertical displacement. In most finite element approximations the number of state variables associated with each shapefunction is constant. In these cases the number of state variables can be a value associated with a mesh. In multiphysics problems, the number of state variables can vary acording to the physical quantity being represented. For instance, in numerical approximations of flow through porous media, two (or three) state variables are associated with the shapefunctions which approximate the displacements of the porous matrix and a single shape function is associated with the pressure variable. This is the reason why each connect keeps track of the number of state variables associated with its shape functions
1.8.13