The Implementation of Shape Functions for CFD

The goal of this project was to implement the shape and derivative functions for the Phasta CFD code in a more complete and consistent manner. Previous implementations of shape and derivative functions were achieved by two different parts of the code and in two separate ways. In addition, not all topologies were implemented to the same polynomial order. This project sought to bring the full implementation of all four topologies into a single implementation, with as much commonality within the code as possible.

There are four different topologies used by Phasta. These topologies – known as the hexahedron, tetrahedron, wedge and pyramid – provide a method for describing the geometry to be used for generating solutions as well as visualization. Traditionally, flow parameters such as momentum and pressure are modeled as constant or linear across a topological body. This assumption simplifies the math used for generating a solution, but also requires a higher concentration of smaller topological bodies in areas where flow parameters are changing rapidly. As an alternative, Phasta can generate more accurate solutions without the need for additional topological bodies by allowing flow parameters to be modeled as higher order polynomial functions.

The shape and derivative functions are used by various sections of Phasta as a means of evaluating flow parameters. Low order functions are easily evaluated across each of the four topological bodies. Shape and derivative functions are needed when it is necessary to evaluate high order functions within a topological body. The shape and derivative functions are not only different for each type of topology, but are also different for the various polynomial orders of the same topology. Since the locations within each topology of the evaluations are known ahead of time and fixed, the shape and derivative functions must be evaluated only once for each topology and each order of accuracy.

The method for calculating these shape functions was laid out in “A Straightforward Structure to Construct Shape Functions for Variable P-Order Meshes”. Unfortunately, the implementation of the functions specified in the paper did not prove as straightforward. The methodology described in the paper models high order shape functions as the contribution from the entities of vertices, edges, faces and regions. In addition, it breaks each shape function down into the pure contribution from the topological entity, and a blending function – unique for each of the four topological bodies – that enforces continuity at the boundaries of the entity. The paper goes on to give mathematical formulas for each of the entity and blending functions on each of the topologies. In addition to the shape functions given in the paper, Phasta also requires that the three spatial partial derivatives of the shape functions be supplied.

To be fully functional, the shape and derivative functions must be able to evaluate and take the derivative of all of the entity and blending functions for all four topologies up to any given order. This was only accomplished for the hexahedron topology. For the tetrahedron, the blending functions were fully implemented. The entity functions on the tetrahedron also work, but not the evaluation of the entity derivatives.

For consistency as well as readability concerns, the new implementation of the shape and derivative functions closely follow the methods described in the paper. This includes not only grouping the shape functions into their contributing entities, but also breaking the shape functions down into entity and blending functions. The code for generating the shape and derivative functions is broken up for each topology, but all four follow a similar framework:

The work of calculating the entity functions is done separately and is discussed later.

An important part of the evaluation of shape and derivative functions is the parametric coordinates used to describe each of the entities. The entity parametric coordinates that are used in the implementation were chosen to remain consistent with the paper. Each parameterization was chosen to simplify the process of calculating the blending and entity functions for that particular topology. This means that an edge or face shared by different types of topologies could potentially be parameterized in two different ways. In order to ensure continuity across such shared edges and faces, the implementation uses a mapping function to convert one type of parameterization into the other before calculating the entity function. Specifically, the mapping function converts what the paper refers to as “parameterization I” into “parameterization II”. After the mapping function it then becomes possible to use the same function to calculate the entity function for both topologies.

Another consequence of the entity parameterization is the volume coordinates used in the tetrahedron parameterization. Although such a parameterization makes the calculation of the shape functions very easy, the calculation of their derivatives is not. The constraint applied to the four tetrahedral coordinates, , results in nontrivial differentiation of a shape function using all four tetrahedral coordinates. This difficulty is overcome by treating all four coordinates as independent and performing the calculation of false derivatives of the blending and entity functions as if no constraint existed. The constraint is finally applied only after the shape function derivatives have been calculated. This is done by setting the true derivatives to their unconstrained false equivalents and then subtracting the false derivative of the fourth coordinate:

In this way, the code that must differentiate the shape function can disregard the constraint caused by the volume coordinates used in the tetrahedron parameterization.

Due to the dissimilarities of the mathematics involved, the entity functions must be implemented separately for each entity on each topology. It is possible that some of this code for the hexahedron or the tetrahedron may be reused for implementation of the entity functions for the wedge or pyramid, but the nature of the parameterization for each of the topologies makes this doubtful. Two and three dimensional entities may contribute more than one entity function for each polynomial order. These multiple entity functions are known as modes and all possible modes must be generated. Once the polynomial order is high enough for a particular entity to contribute an entity function, all modes for all orders up to the one being used are included. The code for generating the entity functions, although separate, can be written in a consistent way:

Because it was based on the same paper, the old implementation of the tetrahedron shape and derivative functions shared a similar structure with the new implementation. This similarity means that some of the old code can be reused in places where it is need. The calculation of the entity function modes corresponds to the evaluation of specific mathematical equations. For the tetrahedron, the edge entity is given by equation 33, the triangular face by equation 36, and the body by equation 42. Since the evaluation of these equations is implementation independent, the new code uses existing code to perform these calculations. The edge, face and body entity functions are evaluated in “En.c”, “Fn.c” and “Bn.c” respectively. As stated earlier, the entity functions for the tetrahedron are fully functional and the integration of the old code works for all orders of accuracy.

Unfortunately, the entity derivative files, “EnDrv.c”, “FnDrv.c” and “BnDrv.c”, are not as easy to integrate into the new implementation. The old code appears to call the entity derivative functions in conjunction with something called “parDrv.v”. The description for the file says that it computes “the derivative of the parametric transformation from element to entity coordinate systems.” In order for the entity function derivatives to be evaluated correctly, it will probably be necessary to make use of this file.

Although the work of implementing the shape and derivative functions for all four topologies is far from complete, this project has laid the groundwork for their completion. Using the framework outlined above to implement the shape and derivative functions makes it easy to follow the methods described in the paper, and using the simplified approach for volume coordinates assists in this process. This project has made it possible to bring the full implementation of all four topologies into a single a file, with a large amount of commonality.