Merge branch 'fix/PerAlign_PrismWithPyramids' into 'master'

PerAlign extension for all types of meshes

See merge request nektar/nektar!1702

.gitlab-ci.yml

@@ -154,7 +154,7 @@ include:
     - apt-get update && apt-get install -y clang-format-$CL_F_VER git-core
   script:
     - cd "$CI_PROJECT_DIR"
-    - find $TARGETS -iname "*.cpp" -o -iname "*.h" -o -iname "*.hpp" | 
+    - find $TARGETS -iname "*.cpp" -o -iname "*.cpp.in" -o -iname "*.h" -o -iname "*.hpp" |
       xargs clang-format-$CL_F_VER --dry-run -Werror
   allow_failure: true
   timeout: 15m

CHANGELOG.md

@@ -47,7 +47,7 @@ v5.6.0
 - Add feature for r-adaption on user-defined CAD curves (!1349)
 - Add unit testing infrustructure and initial example (!1753)
 - Added a custom cmake cache file to load defaults for building only NekMesh without the solvers (!1641)
-
+- Extend peralign module to all types of meshes (!1702)
 
 **IncNavierStokesSolver**
 - Matrix-Free LinearADR operator for VCSImplicit and others (!1627)

library/LibUtilities/BasicConst/GitRevision.cpp.in

@@ -44,6 +44,6 @@ namespace Nektar::NekConstants
 extern const std::string kGitSha1   = "@GIT_SHA1@";
 extern const std::string kGitBranch = "@GIT_REFSPEC@";
 
-}
+} // namespace Nektar::NekConstants
 
 #endif

library/MatrixFreeOps/OperatorImp.cpp.in

@@ -6,6 +6,8 @@
 // shape so to reduce compile time timeout failures due to the large
 // number of template implementations.
 
+// clang-format off
+
 // These defines must be before the operator header file as they are
 // used in the header files.
 #define SHAPE_TYPE_@TYPE@

library/NekMesh/Module/Module.cpp

@@ -624,6 +624,7 @@ void Module::ReorderPrisms(PerMap &perFaces)
     std::unordered_set<int>::iterator fIt[2], fIt2;
 
     // Loop over prisms until we've found all lines of prisms.
+
     while (prismsDone.size() > 0)
     {
         vector<ElementSharedPtr> line;
@@ -758,12 +759,14 @@ void Module::ReorderPrisms(PerMap &perFaces)
                     if (e1->m_n1 == e2->m_n1 && e1->m_n2 == e2->m_n2)
                     {
                         e2->m_edgeNodes = e1->m_edgeNodes;
+                        e2->m_curveType = e1->m_curveType;
                     }
                     else if (e1->m_n1 == e2->m_n1 && e1->m_n2 == e2->m_n2)
                     {
                         e2->m_edgeNodes = e1->m_edgeNodes;
                         std::reverse(e2->m_edgeNodes.begin(),
                                      e2->m_edgeNodes.end());
+                        e2->m_curveType = e1->m_curveType;
                     }
                 }
             }
@@ -792,33 +795,45 @@ void Module::ReorderPrisms(PerMap &perFaces)
                        << "some prisms; this will be ignored in further module "
                        << "evaluations." << endl;
     }
-
+    int maxCouples = 3;
     // Loop over periodic faces, enumerate vertices.
-    for (pIt = perFaces.begin(); pIt != perFaces.end(); ++pIt)
+    for (int flag = 0; flag < maxCouples; flag++)
     {
-        FaceSharedPtr f2     = pIt->second.first;
-        FaceSharedPtr f1     = perFaces[f2->m_id].first;
-        vector<int> perVerts = pIt->second.second;
-        int nVerts           = perVerts.size();
-
-        // Number periodic vertices first.
-        for (j = 0; j < nVerts; ++j)
+        for (pIt = perFaces.begin(); pIt != perFaces.end(); ++pIt)
         {
-            NodeSharedPtr n1 = f1->m_vertexList[j];
-            NodeSharedPtr n2 = f2->m_vertexList[perVerts[j]];
+            FaceSharedPtr f2     = pIt->second.first;
+            FaceSharedPtr f1     = perFaces[f2->m_id].first;
+            vector<int> perVerts = pIt->second.second;
+            int nVerts           = perVerts.size() - 1;
+            int coupleFlags      = perVerts[perVerts.size() - 1];
 
-            if (n1->m_id == -1 && n2->m_id == -1)
-            {
-                n1->m_id = nodeId++;
-                n2->m_id = nodeId++;
-            }
-            else if (n1->m_id != -1 && n2->m_id != -1)
+            if (coupleFlags != flag)
             {
                 continue;
             }
-            else
+            // Number periodic vertices first.
+            for (j = 0; j < nVerts; ++j)
             {
-                ASSERTL0(false, "Periodic face renumbering error");
+                NodeSharedPtr n1 = f1->m_vertexList[j];
+                NodeSharedPtr n2 = f2->m_vertexList[perVerts[j]];
+
+                if (n1->m_id == -1 && n2->m_id == -1)
+                {
+                    n1->m_id = nodeId++;
+                    n2->m_id = nodeId++;
+                }
+                else if (n1->m_id != -1 && n2->m_id != -1)
+                {
+                    continue;
+                }
+                else
+                {
+                    m_log(WARNING)
+                        << "n1 " << n1->m_id << " " << n1->m_x << " " << n1->m_y
+                        << " " << n1->m_z << " n2=" << n2->m_id << " "
+                        << n2->m_x << " " << n2->m_y << " " << n2->m_z << endl;
+                    ASSERTL0(false, "Periodic face renumbering error");
+                }
             }
         }
     }
@@ -851,10 +866,22 @@ void Module::ReorderPrisms(PerMap &perFaces)
     {
         if ((*it)->m_id == -1)
         {
+            m_log(VERBOSE) << "Vertex that is not connected to Prism or Tet in "
+                              "PerAlign id = nodeId++"
+                           << endl;
             (*it)->m_id = nodeId++;
         }
     }
 
+    for (it = m_mesh->m_vertexSet.begin(); it != m_mesh->m_vertexSet.end();
+         ++it)
+    {
+        if ((*it)->m_id == -1)
+        {
+            m_log(FATAL) << "Vetex no ID " << endl;
+        }
+    }
+
     ProcessEdges();
     ProcessFaces();
     ProcessElements();
@@ -877,7 +904,7 @@ void Module::PrismLines(int prism, PerMap &perFaces, set<int> &prismsDone,
     line.push_back(m_mesh->m_element[3][prism]);
 
     // Now find prisms connected to this one through a triangular face.
-    for (i = 1; i <= 3; i += 2)
+    for (i = 1; i <= 3; i += 2) // checks only face 1 and face 3
     {
         FaceSharedPtr f = m_mesh->m_element[3][prism]->GetFace(i);
         int nextId;

solvers/DummySolver/Tests/Dummy_3DCubeCwipi.tst

 <?xml version="1.0" encoding="utf-8"?>
 <test>
     <description>Run two CWIPI-coupled instances of the DummySolver, exchanging fields</description>
-    <segment>
+    <segment type="parallel">
         <executable> DummySolver </executable>
         <parameters> --verbose --cwipi 'Dummy0' Dummy_3DCubeCwipi_0.xml cube.left.xml </parameters>
         <processes> 1 </processes>
     </segment>
-    <segment>
+    <segment type="parallel">
         <executable> DummySolver </executable>
         <parameters> --verbose --cwipi 'Dummy1' Dummy_3DCubeCwipi_1.xml cube.right.xml </parameters>
         <processes> 1 </processes>

tests/TestData.cpp

@@ -153,13 +153,31 @@ Command TestData::ParseCommand(TiXmlElement *elmt) const
 
     cmd.m_pythonTest = false;
 
+    std::string commandType = "none";
+    if (elmt->Attribute("type"))
+    {
+        commandType = elmt->Attribute("type");
+
+        if (commandType == "none")
+        {
+            cmd.m_commandType = eNone;
+        }
+        else if (commandType == "parallel")
+        {
+            cmd.m_commandType = eParallel;
+        }
+        else if (commandType == "sequential")
+        {
+            cmd.m_commandType = eSequential;
+        }
+    }
+
     // Parse executable tag. Do not enforce a check because this might be
     // overridden on the command line.
     if (elmt->FirstChildElement("executable"))
     {
         tmp              = elmt->FirstChildElement("executable");
         cmd.m_executable = fs::path(tmp->GetText());
-
         // Test to see if this test requires Python
         std::string needsPython;
         tmp->QueryStringAttribute("python", &needsPython);
@@ -235,6 +253,12 @@ void TestData::Parse(TiXmlDocument *pDoc)
             m_commands.push_back(ParseCommand(tmp));
             tmp = tmp->NextSiblingElement("segment");
         }
+
+        for (int i = 1; i < m_commands.size(); ++i)
+        {
+            ASSERTL0(m_commands[i].m_commandType == m_commands[0].m_commandType,
+                     "All segment commands should be of the same type.");
+        }
     }
 
     ASSERTL0(m_commands.size() > 0,

tests/TestData.h

@@ -54,12 +54,20 @@ struct DependentFile
     std::string m_filename;
 };
 
+enum CommandType
+{
+    eNone,
+    eSequential,
+    eParallel
+};
+
 struct Command
 {
     fs::path m_executable;
     std::string m_parameters;
     unsigned int m_processes;
     bool m_pythonTest;
+    CommandType m_commandType;
 };
 
 /**

tests/Tester.cpp.in

@@ -44,11 +44,12 @@
  * and errors are appended to @p master.err. These files are sent to all of the
  * metrics for analysis. If the test fails, the output and error files are
  * dumped to the terminal for debugging purposes.
- * 
+ *
  * @see Metric
  * @see Metric#Test:
  */
 
+#include <algorithm>
 #include <chrono>
 #include <fstream>
 #include <iostream>
@@ -76,10 +77,12 @@ namespace po = boost::program_options;
 int setenv(const char *name, const char *value, int overwrite)
 {
     int errcode = 0;
-    if(!overwrite) {
+    if (!overwrite)
+    {
         size_t envsize = 0;
-        errcode = getenv_s(&envsize, NULL, 0, name);
-        if(errcode || envsize) return errcode;
+        errcode        = getenv_s(&envsize, NULL, 0, name);
+        if (errcode || envsize)
+            return errcode;
     }
     return _putenv_s(name, value);
 }
@@ -92,14 +95,12 @@ int main(int argc, char *argv[])
 
     // Set up command line options.
     po::options_description desc("Available options");
-    desc.add_options()
-        ("help,h",                  "Produce this help message.")
-        ("verbose,v",               "Turn on verbosity.")
-        ("generate-metric,g",       po::value<vector<int>>(),
-                                    "Generate a single metric.")
-        ("generate-all-metrics,a",  "Generate all metrics.")
-        ("executable,e",            po::value<string>(),
-                                    "Use specified executable.");
+    desc.add_options()("help,h", "Produce this help message.")(
+        "verbose,v", "Turn on verbosity.")("generate-metric,g",
+                                           po::value<vector<int>>(),
+                                           "Generate a single metric.")(
+        "generate-all-metrics,a", "Generate all metrics.")(
+        "executable,e", po::value<string>(), "Use specified executable.");
 
     po::options_description hidden("Hidden options");
     hidden.add_options()("input-file", po::value<string>(), "Input filename");
@@ -361,7 +362,8 @@ int main(int argc, char *argv[])
 
 #ifdef NEKTAR_TEST_FORCEMPIEXEC
 #else
-                if (cmd.m_processes > 1 || file.GetNumCommands() > 1)
+                if (cmd.m_processes > 1 || (file.GetNumCommands() > 1 &&
+                                            cmd.m_commandType == eParallel))
 #endif
                 {
                     if (mpiAdded)
@@ -397,13 +399,26 @@ int main(int argc, char *argv[])
 
                 // If running with multiple commands simultaneously, separate
                 // with colon.
-                if (j > 0)
+                if (j > 0 && cmd.m_commandType == eParallel)
                 {
                     command += " : ";
                 }
+                else if (j > 0 && cmd.m_commandType == eSequential)
+                {
+                    command += " && ";
+                    if (cmd.m_processes > 1)
+                    {
+                        command += "\"@MPIEXEC@\" ";
+                        if (std::string("@NEKTAR_TEST_USE_HOSTFILE@") == "ON")
+                        {
+                            command += "-hostfile hostfile ";
+                        }
+                    }
+                }
 
                 // Add -n where appropriate.
-                if (file.GetNumCommands() > 1 || cmd.m_processes > 1)
+                if (cmd.m_processes > 1 || (file.GetNumCommands() > 1 &&
+                                            cmd.m_commandType == eParallel))
                 {
                     command += "@MPIEXEC_NUMPROC_FLAG@ ";
                     command += std::to_string(cmd.m_processes) + " ";
@@ -424,7 +439,8 @@ int main(int argc, char *argv[])
                     command += "@PYTHON_EXECUTABLE@ ";
                 }
 
-                command += LibUtilities::PortablePath(execPath);
+                std::string pathString = LibUtilities::PortablePath(execPath);
+                command += pathString;
                 if (HaveOptFile)
                 {
                     command += " --use-opt-file test.opt ";

utilities/NekMesh/CMakeLists.txt

@@ -14,6 +14,9 @@ ADD_NEKTAR_TEST(Nektar++/extrude)
 ADD_NEKTAR_TEST(Nektar++/jac_list_tet_face)
 ADD_NEKTAR_TEST(Nektar++/linearise_invalid_quad)
 ADD_NEKTAR_TEST(Nektar++/linearise_invalid_tet)
+IF (NOT WIN32)
+    ADD_NEKTAR_TEST(Nektar++/peralign_double_periodic_hybrid)
+ENDIF ()
 
 # MeshGen test
 IF(NEKTAR_USE_MESHGEN)
@@ -35,7 +38,6 @@ IF(NEKTAR_USE_MESHGEN)
     IF(NEKTAR_USE_CCM)
         ADD_NEKTAR_TEST(MeshGen/StarCCM/projectcad_ahmed)
         ADD_NEKTAR_TEST(MeshGen/StarCCM/projectcad_pyramids)
-
     ENDIF()
 ENDIF()
 

utilities/NekMesh/Tests/Nektar++/peralign_double_periodic_hybrid.tst

+<?xml version="1.0" encoding="utf-8" ?>
+<test>
+    <description>Test PerAlign module from XML file</description>
+    <segment type="sequential">
+        <executable>NekMesh</executable>
+        <parameters>-m peralign:surf1=3,8:surf2=5,9:dir=y\;0.0,0.0,1.0:orient peralign_hybrid.xml peralign_double_periodic_hybrid.xml </parameters>
+        <processes> 1 </processes>
+    </segment>
+    <segment type="sequential">
+        <executable>../../solvers/IncNavierStokesSolver/IncNavierStokesSolver</executable>
+        <parameters>peralign_double_periodic_hybrid.xml peralign_hybrid_session.xml</parameters>
+        <processes> 1 </processes>
+    </segment>
+    <files>
+        <file description="Input File">peralign_hybrid.xml</file>
+        <file description="Input File">peralign_hybrid_session.xml</file>
+    </files>
+    <metrics>
+        <metric type="L2" id="1">
+            <value variable="u" tolerance="1e-1">0.55</value>
+        </metric>
+    </metrics>
+</test>

utilities/NekMesh/Tests/Nektar++/peralign_hybrid_session.xml

+<?xml version="1.0" encoding="utf-8" ?>
+<NEKTAR>
+    <EXPANSIONS>
+        <E COMPOSITE="C[0,2]" NUMMODES="3" TYPE="MODIFIED" FIELDS="u,v,w" />
+        <E COMPOSITE="C[0,2]" NUMMODES="2" TYPE="MODIFIEDQUADPLUS1" FIELDS="p" />
+    </EXPANSIONS>
+    <CONDITIONS>
+        <SOLVERINFO>
+            <I PROPERTY="SolverType"                 VALUE="VelocityCorrectionScheme"/>
+            <I PROPERTY="EQTYPE"                     VALUE="UnsteadyNavierStokes"    />
+            <I PROPERTY="EvolutionOperator"          VALUE="Nonlinear"               />
+            <I PROPERTY="Projection"                 VALUE="Galerkin"                />
+            <I PROPERTY="GlobalSysSoln"              VALUE="IterativeStaticCond"     />
+            <I PROPERTY="TimeIntegrationMethod"      VALUE="IMEXOrder2"              />
+            <I PROPERTY="SPECTRALHPDEALIASING"       VALUE="True"                    />
+            <I PROPERTY="Driver"                     VALUE="Standard"                />
+            <I PROPERTY="SpectralVanishingViscosity" VALUE="DGKernel"                />
+        </SOLVERINFO>
+
+        <PARAMETERS>
+            <P> TimeStep      = 5e-4 </P>
+            <P> T             = 5    </P>
+            <P> NumSteps      = 2    </P>
+            <P> IO_CheckSteps = 40   </P>
+            <P> IO_InfoSteps  = 10   </P>
+            <P> IO_CFLSteps   = 1    </P>
+            <P> Re            = 1e3  </P>
+            <P> Kinvis        = 1/Re </P>
+        </PARAMETERS>
+
+        <VARIABLES>
+            <V ID="0"> u </V>
+            <V ID="1"> v </V>
+            <V ID="2"> w </V>
+            <V ID="3"> p </V>
+        </VARIABLES>
+
+        <BOUNDARYREGIONS>
+            <B ID="0"> C[7] </B>   <!-- Viscous wall blade        -->
+            <B ID="1"> C[6] </B>   <!-- Inflow                    -->
+            <B ID="2"> C[4] </B>   <!-- Outflow                   -->
+            <B ID="3"> C[8] </B>   <!-- Periodic - A              -->
+            <B ID="4"> C[9] </B>   <!-- Periodic - B              -->
+            <B ID="5"> C[3] </B>   <!-- Periodic - lower          -->
+            <B ID="6"> C[5] </B>   <!-- Periodic - upper          -->
+        </BOUNDARYREGIONS>
+
+        <BOUNDARYCONDITIONS>
+            <REGION REF="0">
+                <D VAR="u" VALUE="0" />
+                <D VAR="v" VALUE="0" />
+                <D VAR="w" VALUE="0" />
+                <N VAR="p" USERDEFINEDTYPE="H" VALUE="0" />
+            </REGION>
+            <REGION REF="1">
+                <D VAR="u" VALUE="0.7" />
+                <D VAR="v" VALUE="0.7" />
+                <D VAR="w" VALUE="0.7" />
+                <N VAR="p" USERDEFINEDTYPE="H" VALUE="0" />
+            </REGION>
+            <REGION REF="2">
+                <N VAR="u" USERDEFINEDTYPE="HOutflow" VALUE="0" />
+                <N VAR="v" USERDEFINEDTYPE="HOutflow" VALUE="0" />
+                <N VAR="w" USERDEFINEDTYPE="HOutflow" VALUE="0" />
+                <D VAR="p" USERDEFINEDTYPE="HOutflow" VALUE="0" />
+            </REGION>
+            <REGION REF="3">  <!-- Periodic (A) -->
+                <P VAR="u" VALUE="[4]" />
+                <P VAR="v" VALUE="[4]" />
+                <P VAR="w" VALUE="[4]" />
+                <P VAR="p" VALUE="[4]" />
+            </REGION>
+            <REGION REF="4">  <!-- Periodic (B) -->
+                <P VAR="u" VALUE="[3]" />
+                <P VAR="v" VALUE="[3]" />
+                <P VAR="w" VALUE="[3]" />
+                <P VAR="p" VALUE="[3]" />
+            </REGION>
+            <REGION REF="5">  <!-- Periodic (up) -->
+                <P VAR="u" VALUE="[6]" />
+                <P VAR="v" VALUE="[6]" />
+                <P VAR="w" VALUE="[6]" />
+                <P VAR="p" VALUE="[6]" />
+            </REGION>
+            <REGION REF="6">  <!-- Periodic (low) -->
+                <P VAR="u" VALUE="[5]" />
+                <P VAR="v" VALUE="[5]" />
+                <P VAR="w" VALUE="[5]" />
+                <P VAR="p" VALUE="[5]" />
+            </REGION>
+        </BOUNDARYCONDITIONS>
+
+        <FUNCTION NAME="InitialConditions">
+            <E VAR="u" VALUE="0.1" />
+            <E VAR="v" VALUE="0.1" />
+            <E VAR="w" VALUE="0.0" />
+            <E VAR="p" VALUE="0.0" />
+        </FUNCTION>
+    </CONDITIONS>
+    <FILTERS>
+        <FILTER TYPE="ModalEnergy">
+            <PARAM NAME="OutputFile">EnergyFile</PARAM>
+            <PARAM NAME="OutputFrequency">1</PARAM>
+        </FILTER>
+    </FILTERS>
+</NEKTAR>