First, I should mention that this is a basic package and is not limited to CFD, it can also be used for other problems.
Finite Element Method (FEM) is a powerful numerical method which has been used successfully for the solution of the existing problems in various scientific and engineering fields such as its application in CFD. Many algorithms have been expressed based on FEM, but none has been used in popular CFD software. In this section, full monopoly is according to Finite Volume Method (FVM) due to better efficiency and adaptability with the physics of problems in comparison with FEM. It doesn’t seem that FEM could compete with FVM unless it was fundamentally changed. In this paper, I am going to show those changes and its result will be a powerful method which has much better performance in all subjects in comparison with FVM and other computational method, I called it Reduced Finite Element Method (RFEM).
The general form of a linear differential equation with boundary conditions can be shown as follows [Zienkiewicz and Morgan (1983); Zienkiewicz and Taylor (2001)]:
(1)
The operators 
and 
in equations (1) can be zero-order, odd-order, even-order or a combination of two or all three
(2)
The weighted residual relationship for these relationships on any element is
(3)
And
(4)
By inserting (2) in (3) we have
(5)
The approximate relation of the field also is equal to
(6)
By inserting (6) in (5) we have
(7)
Together all the elements a comprehensive system of equations is obtained which can be written quite generally as
(8)
And
(9)
Finally by using of equation (7) can write
(10)
And this is the general form of approximation to a differential equation by the finite element method.
New Formulation for Finite Element Method (NFFEM)
Here I will introduce a new formulation for finite element method which its performance is much better than all conventional algorithms of finite element method in CFD. In NFFEM relationship (9) is written as follows:
(11)
And each of the matrix coefficients are
(12)
Where 
become
(13)
For 
we have
(14)
For 
we have
(15)
And for 
we have
(16)
These relations are used for boundary integrals the same way, these are relations completely of NFFEM.
In the equations (13) to (16), by choosing b = 1 full upwind difference scheme (FUDS) and by choosing b = 0 central difference scheme (CDS) is obtained (Flux Vector Splitting Methods (FVSM) can also use with NFFEM, in this form b = 1 is considered and equations are written separately for positive and negative-terms after that the sum of the equations is done). These relationships are also valid for non-linear operators (
), which leads to 
. In the equations presented above, 
is given by
(17)
For the central difference scheme the following equation can also be used instead of 
in equations (13) to (16)
(18)
Where 
is given by
(19)
In this form, the number of participating elements in the equation of node i are limited to two elements in 1D, 2D and 3D. The 
, 
and 
are equal to
(20)
Where
(21)
And αi is the sign of the unknown variable or derivations coefficients (any term that changes the sign of the matrix coefficients)
(22)
Equation (17) can be written with βi and – βi instead of the 
that its result will be 
, depending on the sign of βi that we have chosen, the 
gives the forward or backward differencing scheme.
Superscript k, l and m represent the terms are located on lines k, l and m (The node i in three dimensions is located on three lines, line k at x direction, line l at y direction, and line m at z direction, see figures (1) and (2)).
For hierarchical shape functions, the shape functions dependent to the sides and element as shown in figures (3) and (4) attributable to the points (these points as 
, 
and 
to calculate 
also to determine lines k, l and m are used).
When wind direction is constant before and after discontinuity (shock), as an alternative (or as another option for all position), 
can be used instead of 
, where 
is
. (23)
Where
(24)
Because βi is zero for internal nodes, there is no need to calculate 
, however it can be calculated by the following equation
(25)
Note that NFFEM only for quadrilateral and hexahedron elements can be used. Solving the system of equations of the NFFEM can be performed simply by element by element solving method and line by line sweeping method.
Reduced Finite Element Method (RFEM)
A new method for the finite element formulation presented in the previous section, In this section, I’m going to limit the participating nodes in equation for each node to the nodes that are located on lines k, l and m the node. For this purpose, I introduce reduced elements, characteristics of these elements is as follows:
1 – The different elements are used for different directions of operator, see figures (5) and (6).
2 – For each direction of operator, element only in the same direction has DOF (for one direction of shape function is used p – degree function and for other directions is used zero degree function, see figures (5), (6), (7) and (8)).
3 – For equation of node i, element only on the lines k, l and m the same node has the node, see figures (7) and (8).
In RFEM relationship (11) is written as follows:
(26)
And each of the matrix coefficients are
(27)
Where
(28)
For 
, 
and 
we have
(29)
For 
, 
and 
we have
(30)
For 
, 
and 
we have
(31)
Note: for the 
operator, unknown variable is written in all directions and each direction is no belongs to particular unknown variable. And for 
we have
(32)
Fig. 1. Lines k and l on two fourth-degree elements in two dimensions
Fig. 2. Lines k, l and m on a quadratic element in three dimensions
Fig. 3. Lines k and l also the places that β is calculated on for a hierarchical element in two dimensions
Fig. 4. Lines k, l and m also the places that β is calculated on for a fourth-degree hierarchical element in three dimensions
Fig. 5. A reduced standard quadratic element in two-dimensional, lines k and l also the participating nodes in the equation of node i for each of them
Fig. 6. A reduced hierarchical quadratic element in two-dimensional, lines k and l also the participating nodes in the equation of node i for each of them
Fig. 7. Nodes involved in the equation of node i for FEM and RFEM on several fourth-degree elements in two dimensions
Fig. 8. Nodes involved in the equation of node i for FEM and RFEM on two quadratic elements in three dimensions
Fig. 9. A reduced quadratic element in three-dimensional with two DOF for approximation the mixed derivatives ∂2ϕ/∂x∂y, ∂2ϕ/∂x∂z and ∂2ϕ/∂y∂z
Fig. 12. 1D advective equation ut + ux = 0
Fig. 13. Euler equations 1D based on sod’s shock tube problem (t = 0.2)
Fig. 14. Grid for flow around a cylinder
Fig. 15. Flow around a cylinder (steady-state Euler equations)
Fig. 16. Mach contours for flow around a cylinder (steady-state Euler equations)
Fig. 17. Grid for flow over a NACA 0012 airfoil
Fig. 18. Surface skin friction coefficient distributions on NACA 0012 airfoil (steady-state Navier-Stokes equations)
Figure 19. Surface pressure distribution for subcritical NACA 0012 airfoil (steady-state Euler equations)
Fig. 20. Solution domain and lines grid for flat plate flow
Fig. 21. Skin friction distributions for flat plate flow (steady-state Navier-Stokes equations)
Fig. 22. 1D convection equation with source term, auх = сos(x) on [0, π] and a = 1
For approximating the mixed derivatives we will need to more DOF for example; if in three dimensions, the mixed derivative be in two-direction, like following derivative
(33)
We use the shape function with two DOF in x and y directions, see figure (9). This elements are used only for mixed operator.
When we use the hierarchical shape functions, for some equations have two unknown variable (ϕi and ai, see figure (6)), the additional equation for additional unknown can be written as follows:
(34)
RFEM can be used to isogeometric analysis as well (see the following sections and examples).
A new hybrid difference scheme (HDS) for liner shape functions
In this section, I will introduce a technique to eliminate oscillations for using the FEM and RFEM on liner shape functions, in this case b is written as follows:
(35)
Where 
is “Upwind Parameter” and is chosen in the range 
. It is clear that there are many choices that can be used for. I used the following equation for it
(36)
Where
(37)
And
(38)
The performance of this scheme is much better than other second-order schemes that in the FEM are used for CFD, see figure (11).
A new shape functions for general problems
Equation (35) can be used only for the linear shape functions, in this section another technique is provided for eliminating oscillations that can be used for higher-order shape functions in any degree. This technique works directly on the shape functions (this technique is equivalent to creating new shape functions), for this purpose, the shape functions are written as following
(39)
Where δi is Element Degree Limiter, Field Variable Limiter (EDL, FVL), and N(p) can be any function of p ≥ 2 (for example, the functions of Lagrangian, Serendipity and Bezier for standard shape functions and the functions of Legendre, Chebyshev, Fourier, etc. for hierarchical shape functions). 
, 
and 
in equation (39) are reference functions, original function for them become
(40)
Where 
and 
are liner shape functions.
The δi in equation (39) has two range:
1 – Using the δi as EDL by choosing it in the range 
.
2 – Using the δi as FVL by choosing it in the range 
.
Note: if the reference function is greater than the shape function the sign of EDL and FVL will change.
The EDL can be used only in the directions of the shape function (directions of operator) that βi ≠ 0 is for, and the FVL is used only in the directions of the shape function that βi = 0 is for or inverse, or the FVL and EDL can be added on all nodes with identical or non-identical values, and their relationship is written as follows:
(41)
The value of the 
and 
depending on the application is, see examples. Note that δi is applied to shape functions as one-dimensional (separately for each direction of shape functions).
The following functions can also be used as reference function on the Lagrangian, Serendipity and Bezier functions to make the shape functions of non-oscillatory
(42)
And for boundary elements become
(43)
Where e is a constant coefficient and N0, Nj and NL can be the Lagrangian, Serendipity or Bezier functions and degree them can be chosen 1 or p (for 1 – degree and p – degree the equations and results are different). 
and 
are a constant coefficients and are given by
(44)
Where ξj is the location of the nodes (control points) and ω is the weight coefficient (optional). For p ≥ 3 the hierarchical shape functions with 
can be used. e is chosen as, e = 1 – p for ω = ∞, e = – 1 for ω = 0 and e = 1/1 – p for ω = – ∞.
Another reference function is
(45)
When αi = 1 equation (45) give backward difference approximation and when αi = – 1 give forward difference approximation. A reference function with high performance around discontinuities in derivative form is as follows:
(46)
The difference between equation (46) and equations (40), (42) and (45) on a fourth-degree element around the discontinuity in figure (10).
Fig. 10. The difference between equation (46) and equations (40), (42) and (45) on three fourth-degree elements around the discontinuity with FUDS
Fig. 11. Liner functions for quadratic and cubic IGA functions
The r for higher-order shape functions is written as re (the same for all nodes lines k, l and m)
(47)
Where ri is given by equation (38).
New isogeometric analysis and new functions
In this section, I’m going by reference functions presented in the previous section, I introduce a new isogeometric analysis. The isogeometric analysis [Hughes et al. (2005)] is done by the B-Spline and NURBS functions, although the reference functions presented above are applicable with B-Spline and NURBS functions, here I introduce a new Isogeometric Analysis Functions that are closer to the standard finite element method and more comfortable for the use in CFD and engineering, I made the following algorithm to make these functions on the Bezier functions
(48)
Note that 
are local and can add as hierarchical as well, for p = 2 these functions are equivalent to the B-Spline functions and for p = 3 are similar the PHT-spline functions work of [Deng et al. (2008)] and for more than a third degree are completely new. The equation (48) for first and last boundary elements become
(49)
(50)
For IGA functions, the liner functions to use in equations (40), (42), (45) and (46) are a liner functions between two control points, see figure (10), while the p-order functions for using in equations (42) and (45) are the same IGA functions.
Solutions and Examples
In this section, some examples that have been solved by standard shape functions that are including functions of the Lagrangian and Bezier and hierarchical shape functions that are including functions of the Legendre, Chebyshev (first and second kind) and Fourier Sine Series are presented. It explained that the volume examples were solved to test RFEM is very high and here are just a few of them select and presented. Examples include 1D advective equation, 1D and 2D Euler equations, 2D Navier – Stokes equations and1D convection equation. Here details of the solution of these equations due to the length of paper and be less important are not presented and only the results are shown. Note: All examples are RFEM (the results for the NFFEM is almost identical). I used the infinite elements for non-solid boundaries in all examples were solved.
Conclusions
The best criterion for evaluate a numerical method is the results of the method and as can be seen in the examples were solved, the result of RFEM comparable to the best results were obtained by other methods is (for these examples, because of the small size of the elements, the difference between the results from different methods can not be seen) and given that the system of equations resulting from this approach similar (in terms of density) Finite Difference Method (FDM) is its efficiency can be compared with finite difference methods (although the use of the Legendre shape functions gives an efficiency much higher than FDM) so do not think other methods that are used in CFD be able to compete with it. Also, due to the similarity of relations many of the techniques used in the FDM can be used to RFEM [Hoffmann, and Chiang (1998)].
The work is submitted to the International Scientific Conference “Engineering science and modern manufacture”, France, October, 18–25, 2015, came to the editorial office оn 13.08.2015.