Merge remote-tracking branch 'upstream/master' into feature/meshing_tools

CHANGELOG.md

@@ -31,6 +31,9 @@ v4.4.0
   file (!678)
 - Extend ExtractDataToCoeffs to support interpolation between basis types for
   quads and hexahedra (!682)
+- Enabled MUMPS support in PETSc if a Fortran compiler was found and added 3D
+  support to the Helmholtz smoother used e.g. in FieldConverts C0Projection
+  module (!714)
 
 **ADRSolver:**
 - Add a projection equation system for C^0 projections (!675)
@@ -44,10 +47,14 @@ v4.4.0
 **IncNavierStokesSolver:**
 - Add ability to simulate additional scalar fields (!624)
 - Improve performance when using homogeneous dealiasing (!622)
+- Fix linearised advection for full 3D cases (!708)
+- Added a weak pressure formulation following Guermond & Shen (!713)
+- Added a convective like outflow boundary condition from Dong (!713)
 
 **FieldConvert:**
 - Allow equi-spaced output for 1D and 2DH1D fields (!613)
 - Update quality metric to include scaled Jacobian output (!695)
+- Allow multiple XML files to be specified in InterpField module (!705)
 
 **NekMesh:**
 - Modify curve module to allow for spline input (!628)
@@ -68,6 +75,9 @@ v4.4.0
 - Add flag to `insertsurface` process for non-conforming geometries (!700)
 - Bug fix to get two meshgen regression tests working (!700)
 - Remove libANN in deference to boost::geometry (!703)
+- 2D to 3D mesh extrusion module (!715)
+- Add a mesh extract option to the linearise module to visualise the result 
+  (!712)
 
 **FieldConvert:**
 - Move all modules to a new library, FieldUtils, to support post-processing

CMakeLists.txt

-CMAKE_MINIMUM_REQUIRED(VERSION 2.8.7)
+CMAKE_MINIMUM_REQUIRED(VERSION 2.8.8)
 SET(CMAKE_BUILD_TYPE Release CACHE STRING "Choose the type of build,
 options are: None(CMAKE_CXX_FLAGS or CMAKE_C_FLAGS used) Debug Release
 RelWithDebInfo MinSizeRel.")
 
-PROJECT(Nektar++)
+PROJECT(Nektar++ C CXX)
+
+
+INCLUDE(CheckLanguage)
+
+CHECK_LANGUAGE(Fortran)
+IF(CMAKE_Fortran_COMPILER)
+  ENABLE_LANGUAGE(Fortran)
+ELSE()
+  MESSAGE(STATUS "No Fortran support")
+ENDIF()
 
 # Helps organize projects in IDEs.
 SET_PROPERTY(GLOBAL PROPERTY USE_FOLDERS ON)

cmake/ThirdPartyPETSc.cmake

@@ -30,16 +30,30 @@ IF (NEKTAR_USE_PETSC)
 
         SET(PETSC_C_COMPILER   "${CMAKE_C_COMPILER}")
         SET(PETSC_CXX_COMPILER "${CMAKE_CXX_COMPILER}")
+        SET(PETSC_Fortran_COMPILER "${CMAKE_Fortran_COMPILER}")
 
         IF (NEKTAR_USE_MPI)
             IF (NOT MPI_BUILTIN)
                 SET(PETSC_C_COMPILER   "${MPI_C_COMPILER}")
                 SET(PETSC_CXX_COMPILER "${MPI_CXX_COMPILER}")
+                SET(PETSC_Fortran_COMPILER "${MPI_Fortran_COMPILER}")
             ENDIF (NOT MPI_BUILTIN)
         ELSE (NEKTAR_USE_MPI)
             SET(PETSC_NO_MPI "--with-mpi=0")
         ENDIF (NEKTAR_USE_MPI)
 
+        IF(CMAKE_Fortran_COMPILER)
+            IF(NEKTAR_USE_MPI AND NOT MPI_Fortran_COMPILER)
+                MESSAGE(ERROR "MPI_Fortran_COMPILER not set")
+            ENDIF()
+            # we use a MUMPS build in ordering here, in the future it might make
+            # sense to hook it up with metis/scotch since this MIGHT be faster
+            SET(PETSC_MUMPS --download-scalapack --download-mumps)
+        ELSE()
+            MESSAGE(WARNING "No Fortran support. Building PETSc without MUMPS support")
+            SET(PETSC_Fortran_COMPILER "0")
+        ENDIF()
+
         EXTERNALPROJECT_ADD(
             petsc-3.7.2
             PREFIX ${TPSRC}
@@ -52,17 +66,21 @@ IF (NEKTAR_USE_PETSC)
             URL http://www.nektar.info/thirdparty/petsc-lite-3.7.2.tar.gz
             URL_MD5 "26c2ff8eaaa9e49aea063f839f5daa7e"
             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}
                 --with-shared-libraries=1
+                --with-pic=1
                 --with-x=0
                 --with-ssl=0
                 --prefix=${TPDIST}
                 --with-petsc-arch=c-opt
-                --with-fc=0
+                ${PETSC_MUMPS}
                 ${PETSC_NO_MPI}
             BUILD_COMMAND MAKEFLAGS= make)
 

docs/refs.bib

@@ -345,6 +345,16 @@ year={2011}
   year={2014}
 }
 
+
+@article{Dong15,
+  title={A convective-like energy-stable open boundary condition for simulation of incompressible flows},
+  author={S. Dong},
+  journal={Journal of Computational Physics},
+  volume={302},
+  pages={300-328},
+  year={2015}
+}
+
 @article{Ko07,
   title        = {Vectorized Matlab codes for linear two-dimensional elasticity},
   author       = {Koko, Jonas},
@@ -441,3 +451,12 @@ year={2011}
   publisher    = {Springer London}
 }
 
+@Article{BaPlGrSh16,
+  author =       {Bao, Y.  and Palacios, R. and Graham, M. and Sherwin, S.J. },
+  title =        {Generalized “thick” strip modelling for vortex-induced vibration of long flexible cylinders},
+  journal =      {J. Comp. Phys},
+  year =         {2016},
+  volume =       {321},
+  pages =        {1079-1097},
+}
+

docs/user-guide/utilities/nekmesh.tex

@@ -632,6 +632,20 @@ match the interior domain:
 This module should be added to the module chain for any complex if the user
 suspects there may be a mesh issue. The module will print a warning if there is.
 
+\subsection{2D mesh extrusion}
+
+This module allows a 2D mesh, quads, triangles or both, to be extruded in the
+$z$ direction to make a simple 3D mesh made of prisms and hexahedra. It is also
+capable of extruding the high-order curvature within the 2D mesh. The module
+requires two parameters:
+
+\begin{lstlisting}[style=BashInputStyle]
+  NekMesh -m extrude:layers=n:length=l 2D.xml 3D.xml
+\end{lstlisting}
+
+length which determines how long the $z$ extrusion will be and layers, the
+number of elements in the $z$ direction.
+
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 

library/FieldUtils/ProcessModules/ProcessInterpField.cpp

@@ -90,7 +90,7 @@ void ProcessInterpField::Process(po::variables_map &vm)
     std::vector<std::string> files;
 
     // set up session file for from field
-    files.push_back(m_config["fromxml"].as<string>());
+    ParseUtils::GenerateOrderedStringVector(m_config["fromxml"].as<string>().c_str(), files);
     m_fromField->m_session =
         LibUtilities::SessionReader::CreateInstance(0, 0, files);
 

library/FieldUtils/ProcessModules/ProcessInterpPoints.cpp

@@ -445,8 +445,8 @@ void ProcessInterpPoints::Process(po::variables_map &vm)
     FieldSharedPtr fromField = boost::shared_ptr<Field>(new Field());
 
     std::vector<std::string> files;
+    ParseUtils::GenerateOrderedStringVector(m_config["fromxml"].as<string>().c_str(), files);
     // set up session file for from field
-    files.push_back(m_config["fromxml"].as<string>());
     fromField->m_session =
         LibUtilities::SessionReader::CreateInstance(0, 0, files);
 

library/LibUtilities/BasicUtils/ErrorUtil.hpp

@@ -75,7 +75,7 @@ namespace ErrorUtil
         NekError(const std::string& message) : std::runtime_error(message) {}
     };
         
-    inline static void Error(ErrType type, const char *routine, int lineNumber, const char *msg, unsigned int level)
+    inline static void Error(ErrType type, const char *routine, int lineNumber, const char *msg, unsigned int level, bool DoComm = false)
     {
         // The user of outStream is primarily for the unit tests.
         // The unit tests often generate errors on purpose to make sure
@@ -94,11 +94,14 @@ namespace ErrorUtil
         // the correct rank. Messages are only printed on rank zero.
         int rank = 0;
 #if defined(NEKTAR_USE_MPI)
-        int flag;
-        MPI_Initialized(&flag);
-        if(flag)
+        int flag = 0;
+        if(DoComm)
         {
-            MPI_Comm_rank(MPI_COMM_WORLD,&rank);
+            MPI_Initialized(&flag);
+            if(flag)
+            {
+                MPI_Comm_rank(MPI_COMM_WORLD,&rank);
+            }
         }
 #endif
 
@@ -138,9 +141,12 @@ namespace ErrorUtil
                 }
             }
 #if defined(NEKTAR_USE_MPI)
-            if (flag)
+            if(DoComm)
             {
-                MPI_Barrier(MPI_COMM_WORLD);
+                if (flag)
+                {
+                    MPI_Barrier(MPI_COMM_WORLD);
+                }
             }
 #endif
             throw NekError(baseMsg);
@@ -185,6 +191,10 @@ namespace ErrorUtil
 #define NEKERROR(type, msg) \
     ErrorUtil::Error(type, __FILE__, __LINE__, msg, 0);
 
+
+#define ROOTONLY_NEKERROR(type, msg)                                     \
+    ErrorUtil::Error(type, __FILE__, __LINE__, msg, 0,true);
+
 #define ASSERTL0(condition,msg) \
     if(!(condition)) \
 { \
@@ -211,7 +221,7 @@ namespace ErrorUtil
 #define WARNINGL1(condition,msg) \
     if(!(condition)) \
 { \
-    ErrorUtil::Error(ErrorUtil::ewarning, __FILE__, __LINE__, msg, 0); \
+    ErrorUtil::Error(ErrorUtil::ewarning, __FILE__, __LINE__, msg, 1); \
 }
 
 #else //defined(NEKTAR_DEBUG) || defined(NEKTAR_FULLDEBUG)
@@ -233,7 +243,7 @@ namespace ErrorUtil
 #define WARNINGL2(condition,msg) \
     if(!(condition)) \
 { \
-    ErrorUtil::Error(ErrorUtil::ewarning, __FILE__, __LINE__, msg, 0); \
+    ErrorUtil::Error(ErrorUtil::ewarning, __FILE__, __LINE__, msg, 2); \
 }
 
 #else //NEKTAR_FULLDEBUG

library/LibUtilities/BasicUtils/FieldIO.cpp

@@ -85,6 +85,7 @@ enum FieldIOType {
     eHDF5
 };
 
+
 /**
  * @brief Determine file type of given input file.
  *
@@ -234,7 +235,8 @@ FieldIOSharedPtr FieldIO::CreateForFile(
 void Write(const std::string &outFile,
            std::vector<FieldDefinitionsSharedPtr> &fielddefs,
            std::vector<std::vector<NekDouble> > &fielddata,
-           const FieldMetaDataMap &fieldinfomap)
+           const FieldMetaDataMap &fieldinfomap,
+           const bool backup)
 {
 #ifdef NEKTAR_USE_MPI
     int size;
@@ -256,7 +258,7 @@ void Write(const std::string &outFile,
 #endif
     CommSharedPtr c    = GetCommFactory().CreateInstance("Serial", 0, 0);
     FieldIOSharedPtr f = GetFieldIOFactory().CreateInstance("Xml", c, false);
-    f->Write(outFile, fielddefs, fielddata, fieldinfomap);
+    f->Write(outFile, fielddefs, fielddata, fieldinfomap, backup);
 }
 
 /**
@@ -395,7 +397,7 @@ void FieldIO::AddInfoTag(TagWriterSharedPtr root,
  *
  * @return Absolute path to resulting file.
  */
-std::string FieldIO::SetUpOutput(const std::string outname, bool perRank)
+std::string FieldIO::SetUpOutput(const std::string outname, bool perRank, bool backup)
 {
     ASSERTL0(!outname.empty(), "Empty path given to SetUpOutput()");
 
@@ -406,6 +408,37 @@ std::string FieldIO::SetUpOutput(const std::string outname, bool perRank)
     // serial.
     fs::path specPath(outname), fulloutname;
 
+    // in case we are rank 0 or not on a shared filesystem, check if the specPath already exists
+    if (backup && (rank == 0 || !m_sharedFilesystem) && fs::exists(specPath))
+    {
+        // rename. foo/bar_123.chk -> foo/bar_123.bak0.chk and in case
+        // foo/bar_123.bak0.chk already exists, foo/bar_123.chk -> foo/bar_123.bak1.chk
+        fs::path bakPath = specPath;
+        int cnt = 0;
+        while (fs::exists(bakPath))
+        {
+            bakPath = specPath.parent_path();
+            bakPath += specPath.stem();
+            bakPath += fs::path(".bak" + boost::lexical_cast<std::string>(cnt++));
+            bakPath += specPath.extension();
+        }
+        std::cout << "renaming " << specPath << " -> " << bakPath << std::endl;
+        try
+        {
+            fs::rename(specPath, bakPath);
+        }
+        catch (fs::filesystem_error &e)
+        {
+            ASSERTL0(e.code().value() == berrc::no_such_file_or_directory,
+                        "Filesystem error: " + string(e.what()));
+        }
+    }
+    // wait until rank 0 has backed up the old specPath
+    if (backup)
+    {
+        m_comm->Block();
+    }
+
     if (nprocs == 1)
     {
         fulloutname = specPath;

library/LibUtilities/BasicUtils/FieldIO.h

@@ -184,7 +184,8 @@ LIB_UTILITIES_EXPORT void Write(
     const std::string &outFile,
     std::vector<FieldDefinitionsSharedPtr> &fielddefs,
     std::vector<std::vector<NekDouble> > &fielddata,
-    const FieldMetaDataMap &fieldinfomap = NullFieldMetaDataMap);
+    const FieldMetaDataMap &fieldinfomap = NullFieldMetaDataMap,
+    const bool backup = false);
 LIB_UTILITIES_EXPORT void Import(
     const std::string &infilename,
     std::vector<FieldDefinitionsSharedPtr> &fielddefs,
@@ -234,7 +235,8 @@ public:
         const std::string &outFile,
         std::vector<FieldDefinitionsSharedPtr> &fielddefs,
         std::vector<std::vector<NekDouble> > &fielddata,
-        const FieldMetaDataMap &fieldinfomap = NullFieldMetaDataMap);
+        const FieldMetaDataMap &fieldinfomap = NullFieldMetaDataMap,
+        const bool backup = false);
 
     LIB_UTILITIES_EXPORT inline void Import(
         const std::string &infilename,
@@ -280,14 +282,15 @@ protected:
     }
 
     LIB_UTILITIES_EXPORT std::string SetUpOutput(
-        const std::string outname, bool perRank);
+        const std::string outname, bool perRank, bool backup = false);
 
     /// @copydoc FieldIO::Write
     LIB_UTILITIES_EXPORT virtual void v_Write(
         const std::string                      &outFile,
         std::vector<FieldDefinitionsSharedPtr> &fielddefs,
         std::vector<std::vector<NekDouble> >   &fielddata,
-        const FieldMetaDataMap                 &fieldinfomap) = 0;
+        const FieldMetaDataMap                 &fieldinfomap,
+        const bool                              backup = false) = 0;
 
     /// @copydoc FieldIO::Import
     LIB_UTILITIES_EXPORT virtual void v_Import(
@@ -317,9 +320,10 @@ typedef boost::shared_ptr<FieldIO> FieldIOSharedPtr;
 inline void FieldIO::Write(const std::string                      &outFile,
                            std::vector<FieldDefinitionsSharedPtr> &fielddefs,
                            std::vector<std::vector<NekDouble> > &fielddata,
-                           const FieldMetaDataMap                 &fieldinfomap)
+                           const FieldMetaDataMap                 &fieldinfomap,
+                           const bool                              backup)
 {
-    v_Write(outFile, fielddefs, fielddata, fieldinfomap);
+    v_Write(outFile, fielddefs, fielddata, fieldinfomap, backup);
 }
 
 /**

library/LibUtilities/BasicUtils/FieldIOHdf5.cpp

@@ -173,7 +173,8 @@ FieldIOHdf5::FieldIOHdf5(LibUtilities::CommSharedPtr pComm,
 void FieldIOHdf5::v_Write(const std::string &outFile,
                           std::vector<FieldDefinitionsSharedPtr> &fielddefs,
                           std::vector<std::vector<NekDouble> > &fielddata,
-                          const FieldMetaDataMap &fieldmetadatamap)
+                          const FieldMetaDataMap &fieldmetadatamap,
+                          const bool backup)
 {
     std::stringstream prfx;
     prfx << m_comm->GetRank() << ": FieldIOHdf5::v_Write(): ";
@@ -184,7 +185,7 @@ void FieldIOHdf5::v_Write(const std::string &outFile,
         tm0 = m_comm->Wtime();
     }
 
-    SetUpOutput(outFile, false);
+    SetUpOutput(outFile, false, backup);
 
     // We make a number of assumptions in this code:
     //   1. All element ids have the same type: unsigned int

library/LibUtilities/BasicUtils/FieldIOHdf5.h

@@ -241,7 +241,8 @@ private:
         const std::string &outFile,
         std::vector<FieldDefinitionsSharedPtr> &fielddefs,
         std::vector<std::vector<NekDouble> > &fielddata,
-        const FieldMetaDataMap &fieldinfomap = NullFieldMetaDataMap);
+        const FieldMetaDataMap &fieldinfomap = NullFieldMetaDataMap,
+        const bool backup = false);
 
     LIB_UTILITIES_EXPORT virtual void v_Import(
         const std::string &infilename,

library/LibUtilities/BasicUtils/FieldIOXml.cpp

@@ -87,7 +87,8 @@ FieldIOXml::FieldIOXml(LibUtilities::CommSharedPtr pComm, bool sharedFilesystem)
 void FieldIOXml::v_Write(const std::string &outFile,
                          std::vector<FieldDefinitionsSharedPtr> &fielddefs,
                          std::vector<std::vector<NekDouble> > &fielddata,
-                         const FieldMetaDataMap &fieldmetadatamap)
+                         const FieldMetaDataMap &fieldmetadatamap,
+                         const bool backup)
 {
     double tm0 = 0.0, tm1 = 0.0;
     if (m_comm->GetRank() == 0)
@@ -112,7 +113,7 @@ void FieldIOXml::v_Write(const std::string &outFile,
     // Prepare to write out data. In parallel, we must create directory and
     // determine the full pathname to the file to write out.  Any existing
     // file/directory which is in the way is removed.
-    std::string filename = SetUpOutput(outFile, true);
+    std::string filename = SetUpOutput(outFile, true, backup);
     SetUpFieldMetaData(outFile, fielddefs, fieldmetadatamap);
 
     // Create the file (partition)

library/LibUtilities/BasicUtils/FieldIOXml.h

@@ -242,7 +242,8 @@ private:
         const std::string &outFile,
         std::vector<FieldDefinitionsSharedPtr> &fielddefs,
         std::vector<std::vector<NekDouble> > &fielddata,
-        const FieldMetaDataMap &fieldinfomap = NullFieldMetaDataMap);
+        const FieldMetaDataMap &fieldinfomap = NullFieldMetaDataMap,
+        const bool backup = false);
 
     LIB_UTILITIES_EXPORT virtual DataSourceSharedPtr v_ImportFieldMetaData(
         const std::string &filename, FieldMetaDataMap &fieldmetadatamap);

library/LibUtilities/BasicUtils/PtsIO.cpp

@@ -54,55 +54,6 @@ namespace Nektar
 namespace LibUtilities
 {
 
-void Import(const string &inFile, PtsFieldSharedPtr &ptsField)
-{
-#ifdef NEKTAR_USE_MPI
-    int size;
-    int init;
-    MPI_Initialized(&init);
-
-    // If MPI has been initialised we can check the number of processes
-    // and, if > 1, tell the user he should not be running this
-    // function in parallel. If it is not initialised, we do not
-    // initialise it here, and assume the user knows what they are
-    // doing.
-    if (init)
-    {
-        MPI_Comm_size(MPI_COMM_WORLD, &size);
-        ASSERTL0(size == 1,
-                 "This static function is not available in parallel. Please "
-                 "instantiate a FieldIO object for parallel use.");
-    }
-#endif
-    CommSharedPtr c = GetCommFactory().CreateInstance("Serial", 0, 0);
-    PtsIO p(c);
-    p.Import(inFile, ptsField);
-}
-
-void Write(const string &outFile, const PtsFieldSharedPtr &ptsField)
-{
-#ifdef NEKTAR_USE_MPI
-    int size;
-    int init;
-    MPI_Initialized(&init);
-
-    // If MPI has been initialised we can check the number of processes
-    // and, if > 1, tell the user he should not be running this
-    // function in parallel. If it is not initialised, we do not
-    // initialise it here, and assume the user knows what they are
-    // doing.
-    if (init)
-    {
-        MPI_Comm_size(MPI_COMM_WORLD, &size);
-        ASSERTL0(size == 1,
-                 "This static function is not available in parallel. Please "
-                 "instantiate a FieldIO object for parallel use.");
-    }
-#endif
-    CommSharedPtr c = GetCommFactory().CreateInstance("Serial", 0, 0);
-    PtsIO p(c);
-    p.Write(outFile, ptsField);
-}
 
 PtsIO::PtsIO(CommSharedPtr pComm, bool sharedFilesystem)
     : FieldIOXml(pComm, sharedFilesystem)
@@ -202,12 +153,13 @@ void PtsIO::Import(const string &inFile,
  * @param ptsField  the pts field
  */
 void PtsIO::Write(const string &outFile,
-                  const Nektar::LibUtilities::PtsFieldSharedPtr &ptsField)
+                  const Nektar::LibUtilities::PtsFieldSharedPtr &ptsField,
+                  const bool backup)
 {
     int nTotvars = ptsField->GetNFields() + ptsField->GetDim();
     int np = ptsField->GetNpoints();
 
-    std::string filename = SetUpOutput(outFile, true);
+    std::string filename = SetUpOutput(outFile, true, backup);
     SetUpFieldMetaData(outFile);
 
     // until tinyxml gains support for line break, write the xml manually

library/LibUtilities/BasicUtils/PtsIO.h

@@ -61,12 +61,6 @@ using namespace std;
 typedef std::map<std::string, std::string> PtsMetaDataMap;
 static PtsMetaDataMap NullPtsMetaDataMap;
 
-LIB_UTILITIES_EXPORT void Import(const string &inFile,
-                                 PtsFieldSharedPtr &ptsField);
-
-LIB_UTILITIES_EXPORT void Write(const string &outFile,
-                                const PtsFieldSharedPtr &ptsField);
-
 class PtsIO : public FieldIOXml
 {
 public:
@@ -83,7 +77,8 @@ public:
         FieldMetaDataMap &fieldmetadatamap = NullFieldMetaDataMap);
 
     LIB_UTILITIES_EXPORT void Write(const string &outFile,
-                                    const PtsFieldSharedPtr &ptsField);
+                                    const PtsFieldSharedPtr &ptsField,
+                                    const bool backup = false);
 
     LIB_UTILITIES_EXPORT void ImportFieldData(TiXmlDocument docInput,
                                               PtsFieldSharedPtr &ptsField);

library/LocalRegions/Expansion3D.cpp

@@ -1103,6 +1103,8 @@ namespace Nektar
                      "Not set up for non boundary-interior expansions");
             ASSERTL1(inoutmat->GetRows() == inoutmat->GetColumns(),
                      "Assuming that input matrix was square");
+            ASSERTL1(GetBasisType(0) == LibUtilities::eModified_A,
+                     "Method set up only for modified modal bases curretly");
             
             int i,j;
             int id1,id2;
@@ -1132,6 +1134,9 @@ namespace Nektar
             
             // - if inoutmat.m_rows() == v_NCoeffs() it is a full
             //   matrix system
+
+            // - if inoutmat.m_rows() == v_GetNverts() it is a vertex space
+            //  preconditioner. 
             
             // - if inoutmat.m_rows() == v_NumBndCoeffs() it is a
             //  boundary CG system
@@ -1144,6 +1149,32 @@ namespace Nektar
             {
                 GetFaceToElementMap(face,GetForient(face),map,sign);
             }
+            else if (rows == GetNverts())
+            {
+                int nfvert = faceExp->GetNverts();
+
+                // Need to find where linear vertices are in facemat
+                Array<OneD, unsigned int> linmap;
+                Array<OneD,          int> linsign;
+                
+                // Use a linear expansion to get correct mapping
+                GetLinStdExp()->GetFaceToElementMap(face,GetForient(face),linmap, linsign);
+
+                // zero out sign map to remove all other modes
+                sign = Array<OneD, int> (order_f,0);
+                map  = Array<OneD, unsigned int>(order_f,(unsigned int)0);
+                
+                int fmap; 
+                // Reset sign map to only have contribution from vertices
+                for(i = 0; i < nfvert; ++i)
+                {
+                    fmap = faceExp->GetVertexMap(i,true);
+                    sign[fmap] = 1;
+
+                    // need to reset map 
+                    map[fmap] = linmap[i]; 
+                }
+            }
             else if(rows == NumBndryCoeffs())
             {
                 int nbndry = NumBndryCoeffs();

library/LocalRegions/HexExp.cpp

@@ -567,6 +567,19 @@ namespace Nektar
         }
 
 
+        StdRegions::StdExpansionSharedPtr HexExp::v_GetLinStdExp(void) const
+        {
+            LibUtilities::BasisKey bkey0(m_base[0]->GetBasisType(),
+                           2, m_base[0]->GetPointsKey());
+            LibUtilities::BasisKey bkey1(m_base[1]->GetBasisType(),
+                           2, m_base[1]->GetPointsKey());
+            LibUtilities::BasisKey bkey2(m_base[2]->GetBasisType(),
+                           2, m_base[2]->GetPointsKey());
+            
+            return MemoryManager<StdRegions::StdHexExp>
+                ::AllocateSharedPtr( bkey0, bkey1, bkey2);
+        }
+
         /**
 	 * \brief Retrieves the physical coordinates of a given set of 
          * reference coordinates.

library/LocalRegions/HexExp.h

@@ -160,6 +160,9 @@ namespace Nektar
     
             LOCAL_REGIONS_EXPORT virtual
                 StdRegions::StdExpansionSharedPtr v_GetStdExp(void) const;
+
+            LOCAL_REGIONS_EXPORT virtual
+                StdRegions::StdExpansionSharedPtr v_GetLinStdExp(void) const;
             
             LOCAL_REGIONS_EXPORT virtual void v_ExtractDataToCoeffs(
                 const NekDouble *data,

library/LocalRegions/NodalTriExp.cpp

@@ -669,6 +669,17 @@ namespace Nektar
                                         m_nodalPointsKey.GetPointsType());
         }
 
+        StdRegions::StdExpansionSharedPtr NodalTriExp::v_GetLinStdExp(void) const
+        {
+            LibUtilities::BasisKey bkey0(m_base[0]->GetBasisType(),
+                           2, m_base[0]->GetPointsKey());
+            LibUtilities::BasisKey bkey1(m_base[1]->GetBasisType(),
+                           2, m_base[1]->GetPointsKey());
+
+            return MemoryManager<StdRegions::StdNodalTriExp>
+                ::AllocateSharedPtr( bkey0, bkey1, m_nodalPointsKey.GetPointsType());
+        }
+
         DNekMatSharedPtr NodalTriExp::v_GenMatrix(const StdRegions::StdMatrixKey &mkey)
         {
             DNekMatSharedPtr returnval;

library/LocalRegions/NodalTriExp.h

@@ -178,6 +178,8 @@ namespace Nektar
                                        const StdRegions::StdMatrixKey &mkey); 
             
             virtual StdRegions::StdExpansionSharedPtr v_GetStdExp(void) const;
+
+            virtual StdRegions::StdExpansionSharedPtr v_GetLinStdExp(void) const;
             
             virtual DNekMatSharedPtr v_GenMatrix(const StdRegions::StdMatrixKey &mkey);
 

library/LocalRegions/PrismExp.cpp

@@ -456,6 +456,19 @@ namespace Nektar
                                     m_base[2]->GetBasisKey());
         }
 
+        StdRegions::StdExpansionSharedPtr PrismExp::v_GetLinStdExp(void) const
+        {
+            LibUtilities::BasisKey bkey0(m_base[0]->GetBasisType(),
+                           2, m_base[0]->GetPointsKey());
+            LibUtilities::BasisKey bkey1(m_base[1]->GetBasisType(),
+                           2, m_base[1]->GetPointsKey());
+            LibUtilities::BasisKey bkey2(m_base[2]->GetBasisType(),
+                           2, m_base[2]->GetPointsKey());
+            
+            return MemoryManager<StdRegions::StdPrismExp>
+                ::AllocateSharedPtr( bkey0, bkey1, bkey2);
+        }
+
         /**