## 3 The Control Volume Finite Element Method  
CVFEM is a combination of the FVM and the finite element method (FEM). The standard FVM utilizes a cell-centered scheme, with control volume and the element set as identical, and the variables stored at the control volume centroid. This is shown in Figure 1 with control volume shaded in gray and the variable set at the element centroid represented by a red point. The CVFEM is a vertex-centered scheme that stores the variable of interest at the vertices (nodes) of the mesh elements. The control volume is based on the dual mesh and is defined around the element nodes. The control volumes are thus not identical with the mesh elements; this is portrayed in Figure 1b, where the control volume, gray-shaded, is centered at the vertex. In the CVFEM, control volumes are also constructed at boundary nodes, because variables are stored on nodes; hence, these nodes are part of the solution vector. Boundary control volumes for CVFEM, on the other hand, are shown in Figure 2, where it is obvious that control volumes only partially encapsulate the nodes.  
Integration points (\(ip\)) are defined for the CVFEM on the surface of the dual volume, just as for the CCFVM, where they are also defined at the element face centroid. Integration points, shown at the centroids of the sub-control surfaces in Figure 3a, are used for the numerical integration of surface integrals with the use of arbitrary interpolation functions, unique to each element topology. The value of a generic quantity \(\phi_{ip}\) at an integration point is therefore computed based on the interpolation weights corresponding to the nodal values of \(\phi_{k}\).  
Figure 1: Control volumes representation for (**a**) cell-centered scheme and (**b**) vertex-centered scheme.  
Figure 2: Control volumes for vertex-centered scheme at (**a**) boundary and (**b**) a corner.  
\[\phi_{ip}=\sum_{k=1}^{n}N_{k}|_{ip}\phi_{k} \tag{8}\]  
\(n\) is the number of nodes of the element containing \(ip\) and \(\phi_{k}\) the value at the node of index \(k\) straddling \(ip\), whereas the calculation of the gradient of the quantity \(\phi\) at the integration point is given by  
\[\
abla\phi_{ip}=\sum_{k=1}^{n}\
abla N_{k}|_{ip}\phi_{k} \tag{9}\]  
While interior integration points involve all nodes of the element, only the nodes of a boundary face are required to interpolate the nodal values to the boundary integration point. This is shown in Figure 3a. That is due to the definition of the shape function to be zero at edge-opposite nodes if the integration point is on an edge. Therefore, the value of \(\phi\) at a boundary integration point is  
\[\phi_{bip}=\sum_{k=1}^{n^{s}}N_{k}|_{bip}\phi_{k} \tag{10}\]  
where \(n^{s}\) is the number of boundary (side) nodes straddling \(bip\). However, this is not the case when computing a gradient at a boundary integration point, because there do exist non-zero values for the derivative of shape functions for interior nodes contributing to boundary integration points (Figure 3c). Therefore, the gradient of \(\phi\) at \(bip\) is computed as  
\[\
abla\phi_{bip}=\sum_{k=1}^{n}\
abla N_{k}|_{bip}\phi_{k} \tag{11}\]  
Figure 3: Demonstration of the interpolation process from the nodes of the element to an integration point associated with a dual volume (shaded in gray) for (**a**) interior integration point, (**b**) boundary integration point, and (**c**) boundary integration point associated with the evaluation of a gradient.