I am currently working on a CFD code over a cubic grid. Now, the number of elements used in the simulation is decomposed among the number of processors. Each of those processors (a section of the cube) is then divided up in a number of elements that corresponds to the polynomial order. I would like to understand the mathematical reason of doing so.
What you have is most likely a mesh for a Spectral Element Method or something similar). The solution is represented using a high order polynomial each section of the cube (referred to as an element). For convenience, these polynomials are represented as Lagrange nodal functions, which are associated with a set of points on each element (which are denoted by the vertices of the sub-mesh on each cube).
High order methods tend to be more accurate than a low order method, at least when measured in terms of error as a function of the element size, and are fairly simple to make efficient on modern hardwares.