Merge remote-tracking branch 'upstream/master' into feature/parseutils-simplify

CHANGELOG.md

@@ -14,7 +14,7 @@ v5.0.0
   functionality in FUNCTIONs defined for forcings (!759)
 - Add patch to tinyxml to fix size_t vs int bug (!820)
 - Add ARPACK thirdparty build capabilities (!828)
-- Added native support for csv files in addititon to pts (!760 !835)
+- Added native support for csv files in addititon to pts (!760, !835, !906)
 - Utilize LAPACK_DIR env variable to find the native blas/lapack install (!827)
 - Remove StdExpansion use from MultiRegion (use Expansions instead). (!831)
 - Move steady state check and CFL output from solvers to SolverUtils (!832)
@@ -34,9 +34,11 @@ v5.0.0
 - Fix mac OSX on buildbots (!876)
 - Fix error from (!826) (!876)
 - Fix minor bug in ARPACK thirdparty build cmake (!874)
-- Switch MeshGraph to use factory pattern and add HDF5 geometry support (!900,
-  !904)
-
+- Switch MeshGraph to use factory pattern and add HDF5 geometry support
+  (!900, !904)
+- Fix bug in MeshPartition.cpp which caused incorrect array access when
+  WeightPartitions was used in parallel (!923)
+	
 **NekMesh**:
 - Add feature to read basic 2D geo files as CAD (!731)
 - Add periodic boundary condition meshing in 2D (!733)
@@ -61,6 +63,10 @@ v5.0.0
 - Skip CFI groups of bodies and non-numbered nodes (!891)
 - Add ability to space out 2D BL nodes to better fit local target Delta (!890)
 - Fix automatic peralign call in 2D periodic meshing (!888)
+- Fix BL splitting call from MCF (!910)
+- Support CFI combined lines (!917)
+- Order nodes in Gmsh output (!912)
+- Fix manifold face curvature nodes (!913)
 
 **FieldConvert**:
 - Add input module for Semtex field files (!777)
@@ -86,6 +92,7 @@ v5.0.0
 - Allow performing axi-symmetric Euler simulations (!771)
 - Add ability to use an exponential filtering for stabilization with
   seg, quad and hex elements (!771, !862)
+- Introduce equations of state to account for real gas effects (!880)
 
 **APESolver:**
 - Added two new boundary conditions to the APE system: RiemannInvariantBC
@@ -106,6 +113,8 @@ v4.4.2
 - Fix calculation of BLPoints (!892)
 - Fix deadlock in DiffusionLDG (!885)
 - Fix uninitialised coefficients in DirectFull solver (!898)
+- Updated PETSc to 3.7.7 (!916)
+- Fix typcase to an integer which set Lz < 1 to zero when postprocess hdf5 output (!9922)
 
 **NekMesh**
 - Fix missing periodic boundary meshing and boundary layer mesh adjustment
@@ -115,6 +124,7 @@ v4.4.2
 **Documentation**:
 - Fix sign of the viscous term in the velocity correction scheme equations in
   the user guide (!856)
+- Fixed anonymous clone URL (!909)
 
 **FieldConvert**
 - Allow passing input name with trailing separator (!879)

cmake/ThirdPartyPETSc.cmake

@@ -76,23 +76,22 @@ IF (NEKTAR_USE_PETSC)
         ENDIF()
 
         EXTERNALPROJECT_ADD(
-            petsc-3.7.2
+            petsc-3.7.7
             DEPENDS ${PETSC_DEPS}
             PREFIX ${TPSRC}
             STAMP_DIR ${TPBUILD}/stamp
             DOWNLOAD_DIR ${TPSRC}
-            SOURCE_DIR ${TPBUILD}/petsc-3.7.2
-            TMP_DIR ${TPBUILD}/petsc-3.7.2-tmp
+            SOURCE_DIR ${TPBUILD}/petsc-3.7.7
+            TMP_DIR ${TPBUILD}/petsc-3.7.7-tmp
             INSTALL_DIR ${TPDIST}
-            BINARY_DIR ${TPBUILD}/petsc-3.7.2
-            URL http://www.nektar.info/thirdparty/petsc-lite-3.7.2.tar.gz
-            URL_MD5 "26c2ff8eaaa9e49aea063f839f5daa7e"
+            BINARY_DIR ${TPBUILD}/petsc-3.7.7
+            URL ${TPURL}/petsc-lite-3.7.7.tar.gz
+            URL_MD5 "7b476e38cfab4cddb26f73e5ef23538e"
             CONFIGURE_COMMAND
                 OMPI_FC=${CMAKE_Fortran_COMPILER}
                 OMPI_CC=${CMAKE_C_COMPILER}
                 OMPI_CXX=${CMAKE_CXX_COMPILER}
                 ${PYTHON_EXECUTABLE} ./configure
-                ./configure
                 --with-fc=${PETSC_Fortran_COMPILER}
                 --with-cc=${PETSC_C_COMPILER}
                 --with-cxx=${PETSC_CXX_COMPILER}
@@ -121,7 +120,7 @@ IF (NEKTAR_USE_PETSC)
         ENDIF (NOT PETSC_FOUND)
         SET(PETSC_CONFIG_INCLUDE_DIR ${PETSC_INCLUDES})
         INCLUDE_DIRECTORIES(${PETSC_INCLUDES})
-        ADD_CUSTOM_TARGET(petsc-3.7.2 ALL)
+        ADD_CUSTOM_TARGET(petsc-3.7.7 ALL)
     ENDIF (THIRDPARTY_BUILD_PETSC)
 
     ADD_DEFINITIONS(-DNEKTAR_USING_PETSC)

developer-guide @ b1461b45

-Subproject commit e128cfaffbbd37c734a667cdc2a07b6f06291615
+Subproject commit b1461b45abb08c48397fe6d046249703c4f8f160

docs/user-guide/installation/source.tex

@@ -24,7 +24,7 @@ There are two ways to obtain the source code for \nekpp:
 	public repository. Use this initially if you would like to try using
 	Nektar++ or make local changes to the code.
     \begin{lstlisting}[style=BashInputStyle]
-        git clone http://gitlab.nektar.info/clone/nektar/nektar.git nektar++
+        git clone https://gitlab.nektar.info/nektar/nektar.git nektar++
     \end{lstlisting}
 	\item Using authenticated access. This will allow you to directly contribute
 	back into the code.

docs/user-guide/solvers/compressible-flow.tex

@@ -243,6 +243,7 @@ Under this section it is possible to set the solver information.
   <I PROPERTY="UpwindType"            VALUE="ExactToro"           />
   <I PROPERTY="ProblemType"           VALUE="General"             />
   <I PROPERTY="ViscosityType"         VALUE="Constant"            />
+  <I PROPERTY="EquationOfState"       VALUE="IdealGas"            />
 </SOLVERINFO>
 \end{lstlisting}
 \begin{itemize}
@@ -303,6 +304,15 @@ we want to use for the advection operator:
 \item \inltt{Constant} (Constant viscosity);
 \item \inltt{Variable} (Variable viscosity through the Sutherland's law.);
 \end{itemize}
+\item \inltt{EquationOfState} allows selecting an equation of state for
+accounting for non-ideal gas behaviour:
+\begin{itemize}
+\item \inltt{IdealGas} (default option);
+\item \inltt{VanDerWaals} (requires additional parameters \inltt{Tcrit} and \inltt{Pcrit});
+\item \inltt{RedlichKwong} (requires additional parameters \inltt{Tcrit} and \inltt{Pcrit});
+\item \inltt{PengRobinson} (requires additional parameters \inltt{Tcrit}, \inltt{Pcrit} and
+\inltt{AcentricFactor});
+\end{itemize}
 \end{itemize}
 
 \subsection*{Boundary conditions}

library/Demos/MultiRegions/Tests/Helmholtz3D_CG_Prism_iter_ml_par3.tst

@@ -2,7 +2,7 @@
 <test>
     <description>Helmholtz 3D CG, prisms, Neumann BCs, iterative ML, Par(3)</description>
     <executable>Helmholtz3D</executable>
-    <parameters>--use-metis Helmholtz3D_Prism.xml</parameters>
+    <parameters>--use-metis -I WeightPartitions=Boundary Helmholtz3D_Prism.xml</parameters>
     <processes>3</processes>
     <files>
         <file description="Session File">Helmholtz3D_Prism.xml</file>

library/FieldUtils/InputModules/InputPts.cpp

@@ -86,13 +86,13 @@ void InputPts::Process(po::variables_map &vm)
     string inFile = m_config["infile"].as<string>();
 
     // Determine appropriate field input
-    if (m_f->m_inputfiles.count("pts") != 0)
+    if (m_f->m_inputfiles.count("csv") != 0)
     {
         LibUtilities::CsvIOSharedPtr csvIO =
             MemoryManager<LibUtilities::CsvIO>::AllocateSharedPtr(m_f->m_comm);
         csvIO->Import(inFile, m_f->m_fieldPts);
     }
-    else if (m_f->m_inputfiles.count("csv") != 0)
+    else if (m_f->m_inputfiles.count("pts") != 0)
     {
         LibUtilities::PtsIOSharedPtr ptsIO =
             MemoryManager<LibUtilities::PtsIO>::AllocateSharedPtr(m_f->m_comm);

library/FieldUtils/ProcessModules/ProcessHomogeneousPlane.cpp

@@ -151,6 +151,7 @@ void ProcessHomogeneousPlane::Process(po::variables_map &vm)
         }
         m_f->m_session = LibUtilities::SessionReader::CreateInstance(
             argc, (char **)argv, files, m_f->m_comm->GetRowComm());
+        m_f->m_session->InitSession();
     }
     else
     {

library/LibUtilities/BasicUtils/FieldIOHdf5.cpp

@@ -328,7 +328,7 @@ void FieldIOHdf5::v_Write(const std::string &outFile,
             homoLengths[f].resize(nSubFields);
             for (int sf = 0; sf < nSubFields; ++sf)
             {
-                uint64_t len = fielddefs[f]->m_homogeneousLengths[sf];
+                NekDouble len = fielddefs[f]->m_homogeneousLengths[sf];
                 hashStream << len;
                 homoLengths[f][sf] = len;
             }

library/LibUtilities/CMakeLists.txt

@@ -428,7 +428,7 @@ ENDIF()
 IF( NEKTAR_USE_PETSC )
     TARGET_LINK_LIBRARIES(LibUtilities LINK_PRIVATE ${PETSC_LIBRARIES})
     TARGET_LINK_LIBRARIES(LibUtilities LINK_PUBLIC ${CMAKE_DL_LIBS})
-    ADD_DEPENDENCIES(LibUtilities petsc-3.7.2)
+    ADD_DEPENDENCIES(LibUtilities petsc-3.7.7)
 ENDIF( NEKTAR_USE_PETSC )
 
 # HDF5

library/MultiRegions/CMakeLists.txt

@@ -147,7 +147,7 @@ ENDIF( NEKTAR_USE_MPI )
 
 IF( NEKTAR_USE_PETSC )
     TARGET_LINK_LIBRARIES(MultiRegions LINK_PRIVATE ${PETSC_LIBRARIES})
-    ADD_DEPENDENCIES(MultiRegions petsc-3.7.2)
+    ADD_DEPENDENCIES(MultiRegions petsc-3.7.7)
 ENDIF( NEKTAR_USE_PETSC )
 
 INSTALL(DIRECTORY ./ DESTINATION ${NEKTAR_INCLUDE_DIR}/MultiRegions COMPONENT dev FILES_MATCHING PATTERN "*.h" PATTERN "*.hpp")

library/NekMeshUtils/CADSystem/CFI/CADSystemCFI.cpp

@@ -187,7 +187,31 @@ bool CADSystemCFI::LoadCAD()
             {
                 vector<cfi::Oriented<cfi::TopoEntity *>> *edgeList =
                     face->getChildList();
+
+                // unrolling combined edges
+                vector<cfi::Oriented<cfi::TopoEntity *>> fullEdgeList;
                 for (it2 = edgeList->begin(); it2 != edgeList->end(); it2++)
+                {
+                    cfi::Line *edge = static_cast<cfi::Line *>(it2->entity);
+
+                    if (edge->getTopoSubtype() == cfi::SUBTYPE_COMBINED)
+                    {
+                        vector<cfi::Oriented<cfi::TopoEntity *>> *subEdgeList =
+                            edge->getChildList();
+
+                        for (it3 = subEdgeList->begin();
+                             it3 != subEdgeList->end(); it3++)
+                        {
+                            fullEdgeList.push_back(*it3);
+                        }
+                    }
+                    else
+                    {
+                        fullEdgeList.push_back(*it2);
+                    }
+                }
+
+                for (it2 = fullEdgeList.begin(); it2 != fullEdgeList.end(); it2++)
                 {
                     cfi::Oriented<cfi::TopoEntity *> orientatedEdge = *it2;
                     cfi::Line *edge =
@@ -311,15 +335,45 @@ void CADSystemCFI::AddSurf(int i, cfi::Face *in)
 
     vector<cfi::Oriented<cfi::TopoEntity *>> *edgeList = in->getChildList();
 
+    // unrolling combined edges
+    vector<cfi::Oriented<cfi::TopoEntity *>> fullEdgeList;
+    vector<cfi::Oriented<cfi::TopoEntity *>>::iterator it2, it3;
+    for (it2 = edgeList->begin(); it2 != edgeList->end(); it2++)
+    {
+        cfi::Line *edge = static_cast<cfi::Line *>(it2->entity);
+
+        if (edge->getTopoSubtype() == cfi::SUBTYPE_COMBINED)
+        {
+            vector<cfi::Oriented<cfi::TopoEntity *>> *subEdgeList =
+                edge->getChildList();
+
+            for (it3 = subEdgeList->begin(); it3 != subEdgeList->end(); it3++)
+            {
+                fullEdgeList.push_back(*it3);
+            }
+
+            if (it2->orientation == cfi::ORIENT_NEGATIVE)
+            {
+                reverse(fullEdgeList.begin() + fullEdgeList.size() -
+                            subEdgeList->size(),
+                        fullEdgeList.end());
+            }
+        }
+        else
+        {
+            fullEdgeList.push_back(*it2);
+        }
+    }
+
     vector<EdgeLoopSharedPtr> edgeloops;
     int done = 0;
-    while (done != edgeList->size())
+    while (done != fullEdgeList.size())
     {
         EdgeLoopSharedPtr edgeloop = EdgeLoopSharedPtr(new EdgeLoop);
         string firstVert;
         vector<cfi::Oriented<cfi::TopoEntity *>> *vertList =
-            edgeList->at(done).entity->getChildList();
-        if (edgeList->at(done).orientation == cfi::ORIENT_POSITIVE)
+            fullEdgeList.at(done).entity->getChildList();
+        if (fullEdgeList.at(done).orientation == cfi::ORIENT_POSITIVE)
         {
             firstVert = vertList->at(0).entity->getName();
             edgeloop->edgeo.push_back(CADOrientation::eForwards);
@@ -331,13 +385,13 @@ void CADSystemCFI::AddSurf(int i, cfi::Face *in)
         }
 
         edgeloop->edges.push_back(
-            m_curves[nameToCurveId[edgeList->at(done).entity->getName()]]);
+            m_curves[nameToCurveId[fullEdgeList.at(done).entity->getName()]]);
 
-        for (done++; done < edgeList->size(); done++)
+        for (done++; done < fullEdgeList.size(); done++)
         {
             bool end = false;
-            vertList = edgeList->at(done).entity->getChildList();
-            if (edgeList->at(done).orientation == cfi::ORIENT_POSITIVE)
+            vertList = fullEdgeList.at(done).entity->getChildList();
+            if (fullEdgeList.at(done).orientation == cfi::ORIENT_POSITIVE)
             {
                 if (vertList->at(1).entity->getName() == firstVert)
                 {
@@ -355,7 +409,7 @@ void CADSystemCFI::AddSurf(int i, cfi::Face *in)
             }
 
             edgeloop->edges.push_back(
-                m_curves[nameToCurveId[edgeList->at(done).entity->getName()]]);
+                m_curves[nameToCurveId[fullEdgeList.at(done).entity->getName()]]);
 
             if (end)
             {

library/NekMeshUtils/Module/Module.cpp

@@ -859,7 +859,7 @@ void Module::RegisterConfig(string key, string val)
         {
             it->second.value = it->second.defValue;
         }
-        else 
+        else
         {
             it->second.value = val;
         }

library/NekMeshUtils/SurfaceMeshing/HOSurfaceMesh.cpp

@@ -129,12 +129,21 @@ void HOSurfaceMesh::Process()
 
         FaceSharedPtr f = m_mesh->m_element[2][i]->GetFaceLink();
 
+        bool dumFace = false;
+
         if (!f)
         {
+            //This uses a fake face to build the high-order info
+            //in the case of 2D and manifold geometries without having to
+            //rewrite the 3D code
+            //important to note that face nodes need to be inserted into the
+            //volume nodes of the surface element or they will be forgotton
             f = std::shared_ptr<Face>(new Face(
                 m_mesh->m_element[2][i]->GetVertexList(),
                 vector<NodeSharedPtr>(), m_mesh->m_element[2][i]->GetEdgeList(),
                 LibUtilities::ePolyEvenlySpaced));
+
+            dumFace = true;
         }
 
         f->m_parentCAD = s;
@@ -506,6 +515,12 @@ void HOSurfaceMesh::Process()
             f->m_faceNodes = honodes;
             f->m_curveType = LibUtilities::eGaussLobattoLegendre;
         }
+
+        if(dumFace)
+        {
+            m_mesh->m_element[2][i]->SetVolumeNodes(f->m_faceNodes);
+            m_mesh->m_element[2][i]->SetCurveType(f->m_curveType);
+        }
     }
 
     if (m_mesh->m_verbose)

library/SolverUtils/CMakeLists.txt

@@ -10,7 +10,6 @@ SET(SOLVER_UTILS_SOURCES
   Diffusion/Diffusion.cpp
   Diffusion/Diffusion3DHomogeneous1D.cpp
   Diffusion/DiffusionLDG.cpp
-  Diffusion/DiffusionLDGNS.cpp
   Diffusion/DiffusionLFR.cpp
   Diffusion/DiffusionLFRNS.cpp
   Driver.cpp
@@ -55,7 +54,6 @@ SET(SOLVER_UTILS_HEADERS
   Diffusion/Diffusion.h
   Diffusion/Diffusion3DHomogeneous1D.h
   Diffusion/DiffusionLDG.h
-  Diffusion/DiffusionLDGNS.h
   Diffusion/DiffusionLFR.h
   Diffusion/DiffusionLFRNS.h
   Driver.h

library/SolverUtils/Core/SessionFunction.cpp

@@ -188,6 +188,10 @@ void SessionFunction::Evaluate(std::string pFieldName,
             EvaluateFld(pFieldName, pArray, pTime, domain);
         }
     }
+    else
+    {
+        ASSERTL0(false, "unknown eFunctionType");
+    }
 
     if (m_toCache)
     {
@@ -217,12 +221,16 @@ std::string SessionFunction::Describe(std::string pFieldName, const int domain)
         retVal = ffunc->GetExpression();
     }
     else if (vType == LibUtilities::eFunctionTypeFile ||
-             LibUtilities::eFunctionTypeTransientFile)
+             vType == LibUtilities::eFunctionTypeTransientFile)
     {
         std::string filename =
             m_session->GetFunctionFilename(m_name, pFieldName, domain);
         retVal = "from file " + filename;
     }
+    else
+    {
+        ASSERTL0(false, "unknown eFunctionType");
+    }
 
     return retVal;
 }

library/SolverUtils/Diffusion/DiffusionLDGNS.cpp

-///////////////////////////////////////////////////////////////////////////////
-//
-// File: DiffusionLDGNS.cpp
-//
-// For more information, please see: http://www.nektar.info
-//
-// The MIT License
-//
-// Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
-// Department of Aeronautics, Imperial College London (UK), and Scientific
-// Computing and Imaging Institute, University of Utah (USA).
-//
-// License for the specific language governing rights and limitations under
-// Permission is hereby granted, free of charge, to any person obtaining a
-// copy of this software and associated documentation files (the "Software"),
-// to deal in the Software without restriction, including without limitation
-// the rights to use, copy, modify, merge, publish, distribute, sublicense,
-// and/or sell copies of the Software, and to permit persons to whom the
-// Software is furnished to do so, subject to the following conditions:
-//
-// The above copyright notice and this permission notice shall be included
-// in all copies or substantial portions of the Software.
-//
-// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
-// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
-// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
-// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
-// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
-// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
-// DEALINGS IN THE SOFTWARE.
-//
-// Description: LDGNS diffusion class.
-//
-///////////////////////////////////////////////////////////////////////////////
-
-#include <SolverUtils/Diffusion/DiffusionLDGNS.h>
-#include <iostream>
-#include <iomanip>
-
-#include <LocalRegions/Expansion2D.h>
-
-namespace Nektar
-{
-    namespace SolverUtils
-    {
-        std::string DiffusionLDGNS::type = GetDiffusionFactory().
-        RegisterCreatorFunction("LDGNS", DiffusionLDGNS::create);
-        
-        DiffusionLDGNS::DiffusionLDGNS()
-        {
-        }
-        
-        void DiffusionLDGNS::v_InitObject(
-            LibUtilities::SessionReaderSharedPtr        pSession,
-            Array<OneD, MultiRegions::ExpListSharedPtr> pFields)
-        {
-            m_session = pSession;
-            m_session->LoadParameter ("Gamma",         m_gamma, 1.4);
-            m_session->LoadParameter ("GasConstant",   m_gasConstant, 287.058);
-            m_session->LoadParameter ("Twall",         m_Twall, 300.15);
-            m_session->LoadSolverInfo("ViscosityType", m_ViscosityType, 
-                                      "Constant");
-            m_session->LoadParameter ("mu",            m_mu, 1.78e-05);
-            m_session->LoadParameter ("thermalConductivity",
-                                      m_thermalConductivity, 0.0257);
-            m_session->LoadParameter ("rhoInf",        m_rhoInf, 1.225);
-            m_session->LoadParameter ("pInf",          m_pInf, 101325);
-            
-            // Setting up the normals
-            int i;
-            int nDim = pFields[0]->GetCoordim(0);
-            int nTracePts = pFields[0]->GetTrace()->GetTotPoints();
-            
-            m_spaceDim = nDim;
-            if (pSession->DefinesSolverInfo("HOMOGENEOUS"))
-            {
-                m_spaceDim = 3;
-            }
-            
-            m_diffDim = m_spaceDim - nDim;
-
-            m_traceVel = Array<OneD, Array<OneD, NekDouble> >(m_spaceDim);
-            m_traceNormals = Array<OneD, Array<OneD, NekDouble> >(m_spaceDim);
-            for(i = 0; i < m_spaceDim; ++i)
-            {
-                m_traceVel[i] = Array<OneD, NekDouble> (nTracePts, 0.0);
-                m_traceNormals[i] = Array<OneD, NekDouble> (nTracePts);
-            }
-            pFields[0]->GetTrace()->GetNormals(m_traceNormals);
-        }
-        
-        /**
-         * @brief Calculate weak DG Diffusion in the LDG form for the 
-         * Navier-Stokes (NS) equations:
-         *
-         * \f$ \langle\psi, \hat{u}\cdot n\rangle 
-         *   - \langle\nabla\psi \cdot u\rangle
-         *     \langle\phi, \hat{q}\cdot n\rangle - 
-         *     (\nabla \phi \cdot q) \rangle \f$
-         *
-         * The equations that need a diffusion operator are those related 
-         * with the velocities and with the energy.
-         *
-         */
-        void DiffusionLDGNS::v_Diffuse(
-            const int                                         nConvectiveFields,
-            const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
-            const Array<OneD, Array<OneD, NekDouble> >        &inarray,
-                  Array<OneD, Array<OneD, NekDouble> >        &outarray,
-            const Array<OneD, Array<OneD, NekDouble> >        &pFwd,
-            const Array<OneD, Array<OneD, NekDouble> >        &pBwd)
-        {
-            int i, j;
-            int nDim      = fields[0]->GetCoordim(0);
-            int nScalars  = inarray.num_elements();
-            int nPts      = fields[0]->GetTotPoints();
-            int nCoeffs   = fields[0]->GetNcoeffs();
-            int nTracePts = fields[0]->GetTrace()->GetTotPoints();
-
-            Array<OneD, NekDouble>               tmp1(nCoeffs);
-            Array<OneD, Array<OneD, NekDouble> > tmp2(nConvectiveFields);
-            
-            Array<OneD, Array<OneD, Array<OneD, NekDouble> > > 
-                                                    numericalFluxO1(m_spaceDim);
-            Array<OneD, Array<OneD, Array<OneD, NekDouble> > > 
-                                                    derivativesO1(m_spaceDim);
-
-            for (j = 0; j < m_spaceDim; ++j)
-            {
-                numericalFluxO1[j] = Array<OneD, Array<OneD, NekDouble> >(
-                                                                    nScalars);
-                derivativesO1[j]   = Array<OneD, Array<OneD, NekDouble> >(
-                                                                    nScalars);
-                
-                for (i = 0; i < nScalars; ++i)
-                {
-                    numericalFluxO1[j][i] = Array<OneD, NekDouble>(
-                                                                nTracePts, 0.0);
-                    derivativesO1[j][i]   = Array<OneD, NekDouble>(nPts, 0.0);
-                }
-            }
-            
-            // Compute the numerical fluxes for the first order derivatives
-            v_NumericalFluxO1(fields, inarray, numericalFluxO1, pFwd, pBwd);
-
-            for (j = 0; j < nDim; ++j)
-            {
-                for (i = 0; i < nScalars; ++i)
-                {
-                    fields[i]->IProductWRTDerivBase (j, inarray[i], tmp1);
-                    Vmath::Neg                      (nCoeffs, tmp1, 1);
-                    fields[i]->AddTraceIntegral     (numericalFluxO1[j][i],