Merge branch 'fix/nm-reainput' into 'master'

Fix issue with older rea input files in 2D

See merge request !765

CHANGELOG.md

@@ -31,6 +31,7 @@ v4.4.1
 - Rework meshing control so that if possible viewable meshes will be dumped
   when some part of the system fails (!756)
 - Add manifold meshing option (!756)
+- Fix issue with older rea input files (!765)
 
 **FieldConvert:**
 - Fix issue with field ordering in the interppointdatatofld module (!754)

utilities/NekMesh/CMakeLists.txt

@@ -136,6 +136,7 @@ ADD_NEKTAR_TEST        (Gmsh/SquareTriLinear)
 # Nektar tests
 ADD_NEKTAR_TEST        (Nektar/HexLinear)
 ADD_NEKTAR_TEST        (Nektar/Tube45)
+ADD_NEKTAR_TEST        (Nektar/UKMesh)
 # Nek5000 tests
 ADD_NEKTAR_TEST        (Nek5000/3da)
 ADD_NEKTAR_TEST        (Nek5000/r1854a)

utilities/NekMesh/InputModules/InputNek.cpp

@@ -59,6 +59,8 @@ ModuleKey InputNek::className = GetModuleFactory().RegisterCreatorFunction(
 
 InputNek::InputNek(MeshSharedPtr m) : InputModule(m)
 {
+    m_config["scalar"] = ConfigOption(
+        true, "0", "If defined then assume input rea is for scalar problem");
 }
 
 InputNek::~InputNek()
@@ -91,6 +93,8 @@ void InputNek::Process()
     boost::unordered_map<int, int> elMap;
     vector<LibUtilities::ShapeType> elmOrder;
 
+    bool scalar = m_config["scalar"].as<bool>();
+
     // Set up vector of processing orders.
     elmOrder.push_back(LibUtilities::eSegment);
     elmOrder.push_back(LibUtilities::eTriangle);
@@ -176,12 +180,16 @@ void InputNek::Process()
 
     // Set up field names.
     m_mesh->m_fields.push_back("u");
-    m_mesh->m_fields.push_back("v");
-    if (m_mesh->m_spaceDim > 2)
+
+    if (!scalar)
     {
-        m_mesh->m_fields.push_back("w");
+        m_mesh->m_fields.push_back("v");
+        if (m_mesh->m_spaceDim > 2)
+        {
+            m_mesh->m_fields.push_back("w");
+        }
+        m_mesh->m_fields.push_back("p");
     }
-    m_mesh->m_fields.push_back("p");
 
     // Loop over and create elements.
     for (i = 0; i < nElements; ++i)
@@ -663,6 +671,15 @@ void InputNek::Process()
         vector<ConditionType> type;
         ConditionSharedPtr c = MemoryManager<Condition>::AllocateSharedPtr();
 
+        ElementSharedPtr elm = m_mesh->m_element[m_mesh->m_spaceDim][elId];
+
+        // Ignore BCs for undefined edges/faces
+        if ((elm->GetDim() == 2 && faceId >= elm->GetEdgeCount()) ||
+            (elm->GetDim() == 3 && faceId >= elm->GetFaceCount()))
+        {
+            continue;
+        }
+
         // First character on each line describes type of BC. Currently
         // only support V, W, and O. In this switch statement we
         // construct the quantities needed to search for the condition.
@@ -671,14 +688,22 @@ void InputNek::Process()
             // Wall boundary.
             case 'W':
             {
-                for (i = 0; i < m_mesh->m_fields.size() - 1; ++i)
+                if (scalar)
                 {
                     vals.push_back("0");
                     type.push_back(eDirichlet);
                 }
-                // Set high-order boundary condition for wall.
-                vals.push_back("0");
-                type.push_back(eHOPCondition);
+                else
+                {
+                    for (i = 0; i < m_mesh->m_fields.size() - 1; ++i)
+                    {
+                        vals.push_back("0");
+                        type.push_back(eDirichlet);
+                    }
+                    // Set high-order boundary condition for wall.
+                    vals.push_back("0");
+                    type.push_back(eHOPCondition);
+                }
                 break;
             }
 
@@ -687,7 +712,7 @@ void InputNek::Process()
             case 'V':
             case 'v':
             {
-                for (i = 0; i < m_mesh->m_fields.size() - 1; ++i)
+                if (scalar)
                 {
                     getline(m_mshFile, line);
                     size_t p = line.find_first_of('=');
@@ -695,23 +720,42 @@ void InputNek::Process()
                         boost::algorithm::trim_copy(line.substr(p + 1)));
                     type.push_back(eDirichlet);
                 }
-                // Set high-order boundary condition for Dirichlet
-                // condition.
-                vals.push_back("0");
-                type.push_back(eHOPCondition);
+                else
+                {
+                    for (i = 0; i < m_mesh->m_fields.size() - 1; ++i)
+                    {
+                        getline(m_mshFile, line);
+                        size_t p = line.find_first_of('=');
+                        vals.push_back(
+                            boost::algorithm::trim_copy(line.substr(p + 1)));
+                        type.push_back(eDirichlet);
+                    }
+                    // Set high-order boundary condition for Dirichlet
+                    // condition.
+                    vals.push_back("0");
+                    type.push_back(eHOPCondition);
+                }
                 break;
             }
 
             // Natural outflow condition (default value = 0.0?)
             case 'O':
             {
-                for (i = 0; i < m_mesh->m_fields.size(); ++i)
+                if (scalar)
                 {
                     vals.push_back("0");
                     type.push_back(eNeumann);
                 }
-                // Set zero Dirichlet condition for outflow.
-                type[m_mesh->m_fields.size() - 1] = eDirichlet;
+                else
+                {
+                    for (i = 0; i < m_mesh->m_fields.size(); ++i)
+                    {
+                        vals.push_back("0");
+                        type.push_back(eNeumann);
+                    }
+                    // Set zero Dirichlet condition for outflow.
+                    type[m_mesh->m_fields.size() - 1] = eDirichlet;
+                }
                 break;
             }
 
@@ -750,7 +794,6 @@ void InputNek::Process()
         }
 
         int compTag, conditionId;
-        ElementSharedPtr elm = m_mesh->m_element[m_mesh->m_spaceDim][elId];
         ElementSharedPtr surfEl;
 
         // Create element for face (3D) or segment (2D). At the moment
@@ -805,7 +848,7 @@ void InputNek::Process()
                 surfEl->GetEdge(i)->m_curveType = f->m_edgeList[i]->m_curveType;
             }
         }
-        else
+        else if (faceId < elm->GetEdgeCount())
         {
             EdgeSharedPtr f = elm->GetEdge(faceId);
 
@@ -825,6 +868,11 @@ void InputNek::Process()
                 LibUtilities::eSegment, conf, nodeList, tags);
         }
 
+        if (!surfEl)
+        {
+            continue;
+        }
+
         LibUtilities::ShapeType surfElType = surfEl->GetConf().m_e;
 
         if (!found)

utilities/NekMesh/Tests/Nektar/UKMesh.tst

+<?xml version="1.0" encoding="utf-8" ?>
+<test>
+    <description>Nektar triangular mesh using older-style rea</description>
+    <executable>NekMesh</executable>
+    <parameters>-m jac:list UKMesh.rea.gz:rea:scalar UKMesh.xml:xml:test</parameters>
+    <files>
+        <file description="Input File">UKMesh.rea.gz</file>
+    </files>
+    <metrics>
+        <metric type="regex" id="1">
+            <regex>^Total negative Jacobians: (\d+)</regex>
+            <matches>
+                <match>
+                    <field id="0">0</field>
+                </match>
+            </matches>
+        </metric>
+    </metrics>
+</test>