Merge branch 'feature/dump-mesh-on-failure' into 'master'

feature/dump-mesh-on-failure

See merge request !756

CHANGELOG.md

@@ -4,7 +4,6 @@ Changelog
 v4.5.0
 ------
 **NekMesh**:
-
 - Add periodic boundary condition meshing in 2D (!733)
 - Adjust boundary layer thickness in corners in 2D (!739)
 
@@ -13,8 +12,14 @@ v4.5.0
 
 v4.4.1
 ------
+**Library**
+- Remove the duplicate output of errorutil (!756)
+
 **NekMesh**:
 - Fix memory consumption issue with Gmsh output (!747)
+- 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)
 
 v4.4.0
 ------

library/NekMeshUtils/CADSystem/OCE/CADSurfOCE.cpp

@@ -90,11 +90,9 @@ Array<OneD, NekDouble> CADSurfOCE::locuv(Array<OneD, NekDouble> p)
     uvr[1]      = p2.Y();
 
     gp_Pnt p3 = m_sas->Value(p2);
-    if (p3.Distance(loc) > 1e-6)
-    {
-        cout << "large locuv distance " << p3.Distance(loc)/1000.0 << " " << m_id
-             << endl;
-    }
+    WARNINGL2(p3.Distance(loc) < 1e-3, "large locuv distance " +
+                boost::lexical_cast<string>(p3.Distance(loc)/1000.0) + " " +
+                boost::lexical_cast<string>(m_id));
 
     // if the uv returned is slightly off the surface
     //(which ShapeAnalysis_Surface can do sometimes)

library/NekMeshUtils/Octree/Octree.cpp

@@ -959,8 +959,8 @@ void Octree::CompileSourcePointList()
 
         // these are the acutal number of sample points in each parametric
         // direction
-        int nu = ceil(DeltaU / m_minDelta) * 2;
-        int nv = ceil(DeltaV / m_minDelta) * 2;
+        int nu = ceil(DeltaU / m_minDelta) * 40 * 2;
+        int nv = ceil(DeltaV / m_minDelta) * 40 * 2;
 
         for (int j = 0; j < nu; j++)
         {

utilities/NekMesh/InputModules/InputMCF.cpp

@@ -124,10 +124,10 @@ void InputMCF::ParseFile(string nm)
     }
 
     set<string> refinement;
-    if(pSession->DefinesElement("NEKTAR/MESHING/REFINEMENT"))
+    if (pSession->DefinesElement("NEKTAR/MESHING/REFINEMENT"))
     {
         TiXmlElement *refine = mcf->FirstChildElement("REFINEMENT");
-        TiXmlElement *L     = refine->FirstChildElement("LINE");
+        TiXmlElement *L      = refine->FirstChildElement("LINE");
 
         while (L)
         {
@@ -179,18 +179,18 @@ void InputMCF::ParseFile(string nm)
     it = information.find("MeshType");
     ASSERTL0(it != information.end(), "no meshtype defined");
     m_makeBL = it->second == "3DBndLayer";
-    m_2D = it->second == "2D";
+    m_2D     = it->second == "2D";
+    m_manifold = it->second == "Manifold";
     if (it->second == "2DBndLayer")
     {
         m_makeBL = true;
-        m_2D = true;
+        m_2D     = true;
     }
-    if(!m_makeBL && !m_2D)
+    if (!m_makeBL && !m_2D && !m_manifold)
     {
         ASSERTL0(it->second == "3D", "unsure on MeshType")
     }
 
-
     it = parameters.find("MinDelta");
     ASSERTL0(it != parameters.end(), "no mindelta defined");
     m_minDelta = it->second;
@@ -217,13 +217,13 @@ void InputMCF::ParseFile(string nm)
         ASSERTL0(it != parameters.end(), "no BndLayerthick defined");
         m_blthick = it->second;
 
-        it = parameters.find("BndLayerLayers");
+        it        = parameters.find("BndLayerLayers");
         m_splitBL = it != parameters.end();
-        if(m_splitBL)
+        if (m_splitBL)
         {
             m_bllayers = it->second;
-            it = parameters.find("BndLayerProgression");
-            m_blprog = it != parameters.end() ? it->second : "2.0";
+            it         = parameters.find("BndLayerProgression");
+            m_blprog   = it != parameters.end() ? it->second : "2.0";
         }
 
         it = parameters.find("BndLayerAdjustment");
@@ -239,7 +239,7 @@ void InputMCF::ParseFile(string nm)
     }
 
     m_naca = false;
-    if(m_2D && m_cadfile.find('.') == std::string::npos)
+    if (m_2D && m_cadfile.find('.') == std::string::npos)
     {
         m_naca = true;
 
@@ -274,16 +274,16 @@ void InputMCF::ParseFile(string nm)
     m_adjustall = sit != boolparameters.end();
 
     m_refine = refinement.size() > 0;
-    if(m_refine)
+    if (m_refine)
     {
         stringstream ss;
-        for(sit = refinement.begin(); sit != refinement.end(); sit++)
+        for (sit = refinement.begin(); sit != refinement.end(); sit++)
         {
             ss << *sit;
             ss << ":";
         }
         m_refinement = ss.str();
-        m_refinement.erase(m_refinement.end()-1);
+        m_refinement.erase(m_refinement.end() - 1);
     }
 
     if (periodic.size() > 0)
@@ -307,114 +307,214 @@ void InputMCF::Process()
     m_mesh->m_spaceDim = 3;
     m_mesh->m_nummode  = boost::lexical_cast<int>(m_order) + 1;
 
-    vector<ModuleSharedPtr> mods;
+    ModuleSharedPtr module;
 
     ////**** CAD ****////
-    mods.push_back(GetModuleFactory().CreateInstance(
-        ModuleKey(eProcessModule, "loadcad"), m_mesh));
-    mods.back()->RegisterConfig("filename", m_cadfile);
+    module = GetModuleFactory().CreateInstance(
+        ModuleKey(eProcessModule, "loadcad"), m_mesh);
+    module->RegisterConfig("filename", m_cadfile);
 
-    if(m_2D)
+    if (m_2D)
     {
-        mods.back()->RegisterConfig("2D","");
+        module->RegisterConfig("2D", "");
     }
-    if(m_naca)
+    if (m_naca)
     {
-        mods.back()->RegisterConfig("NACA",m_nacadomain);
+        module->RegisterConfig("NACA", m_nacadomain);
     }
 
-    ////**** Octree ****////
-    mods.push_back(GetModuleFactory().CreateInstance(
-        ModuleKey(eProcessModule, "loadoctree"), m_mesh));
-    mods.back()->RegisterConfig("mindel", m_minDelta);
-    mods.back()->RegisterConfig("maxdel", m_maxDelta);
-    mods.back()->RegisterConfig("eps", m_eps);
+    module->SetDefaults();
+    module->Process();
+
+    ////**** OCTREE ****////
+    module = GetModuleFactory().CreateInstance(
+        ModuleKey(eProcessModule, "loadoctree"), m_mesh);
+    module->RegisterConfig("mindel", m_minDelta);
+    module->RegisterConfig("maxdel", m_maxDelta);
+    module->RegisterConfig("eps", m_eps);
     if (m_refine)
     {
-        mods.back()->RegisterConfig("refinement", m_refinement);
+        module->RegisterConfig("refinement", m_refinement);
     }
     if (m_woct)
     {
-        mods.back()->RegisterConfig("writeoctree", "");
+        module->RegisterConfig("writeoctree", "");
     }
 
-    if(m_2D)
+    module->SetDefaults();
+    module->Process();
+
+    ////**** LINEAR MESHING ****////
+    if (m_2D)
     {
         ////**** 2DGenerator ****////
         m_mesh->m_expDim = 2;
         m_mesh->m_spaceDim = 2;
-        mods.push_back(GetModuleFactory().CreateInstance(
-            ModuleKey(eProcessModule, "2dgenerator"), m_mesh));
+        module             = GetModuleFactory().CreateInstance(
+            ModuleKey(eProcessModule, "2dgenerator"), m_mesh);
         if (m_makeBL)
         {
-            mods.back()->RegisterConfig("blcurves", m_blsurfs);
-            mods.back()->RegisterConfig("blthick", m_blthick);
-
-            if (m_adjust)
-            {
-                mods.back()->RegisterConfig("bltadjust", m_adjustment);
+            module->RegisterConfig("blcurves", m_blsurfs);
+            module->RegisterConfig("blthick", m_blthick);
+        }
 
-                if (m_adjustall)
-                {
-                    mods.back()->RegisterConfig("adjustblteverywhere", "");
-                }
-            }
+        try
+        {
+            module->SetDefaults();
+            module->Process();
         }
-        if (m_periodic.size())
+        catch (runtime_error &e)
         {
-            mods.back()->RegisterConfig("periodic", m_periodic);
+            cout << "2D linear mesh generator failed with message:" << endl;
+            cout << e.what() << endl;
+            cout << "No mesh file has been created" << endl;
+            abort();
         }
     }
     else
     {
         ////**** SurfaceMesh ****////
-        mods.push_back(GetModuleFactory().CreateInstance(
-            ModuleKey(eProcessModule, "surfacemesh"), m_mesh));
+        module = GetModuleFactory().CreateInstance(
+            ModuleKey(eProcessModule, "surfacemesh"), m_mesh);
+
+        try
+        {
+            module->SetDefaults();
+            module->Process();
+        }
+        catch (runtime_error &e)
+        {
+            cout << "Surface meshing has failed with message:" << endl;
+            cout << e.what() << endl;
+            cout << "Any surfaces which were succsessfully meshed will be "
+                    "dumped as a manifold mesh"
+                 << endl;
+            m_mesh->m_expDim = 2;
+            ProcessVertices();
+            ProcessEdges();
+            ProcessFaces();
+            ProcessElements();
+            ProcessComposites();
+            return;
+        }
 
-        ////**** VolumeMesh ****////
-        mods.push_back(GetModuleFactory().CreateInstance(
-            ModuleKey(eProcessModule, "volumemesh"), m_mesh));
-        if(m_makeBL)
+        if(m_manifold)
+        {
+            //dont want to volume mesh
+            m_mesh->m_expDim = 2;
+        }
+        else
         {
-            mods.back()->RegisterConfig("blsurfs",m_blsurfs);
-            mods.back()->RegisterConfig("blthick",m_blthick);
-            mods.back()->RegisterConfig("bllayers",m_bllayers);
-            mods.back()->RegisterConfig("blprog",m_blprog);
+            ////**** VolumeMesh ****////
+            module = GetModuleFactory().CreateInstance(
+                ModuleKey(eProcessModule, "volumemesh"), m_mesh);
+            if (m_makeBL)
+            {
+                module->RegisterConfig("blsurfs", m_blsurfs);
+                module->RegisterConfig("blthick", m_blthick);
+                module->RegisterConfig("bllayers", m_bllayers);
+                module->RegisterConfig("blprog", m_blprog);
+            }
+
+            try
+            {
+                module->SetDefaults();
+                module->Process();
+            }
+            catch (runtime_error &e)
+            {
+                cout << "Volume meshing has failed with message:" << endl;
+                cout << e.what() << endl;
+                cout << "The linear surface mesh be dumped as a manifold mesh"
+                     << endl;
+                m_mesh->m_expDim = 2;
+                m_mesh->m_element[3].clear();
+                ProcessVertices();
+                ProcessEdges();
+                ProcessFaces();
+                ProcessElements();
+                ProcessComposites();
+                return;
+            }
         }
     }
 
     ////**** HOSurface ****////
-    mods.push_back(GetModuleFactory().CreateInstance(
-        ModuleKey(eProcessModule, "hosurface"), m_mesh));
+    module = GetModuleFactory().CreateInstance(
+        ModuleKey(eProcessModule, "hosurface"), m_mesh);
     if (m_surfopti)
     {
-        mods.back()->RegisterConfig("opti", "");
+        module->RegisterConfig("opti", "");
+    }
+
+    try
+    {
+        module->SetDefaults();
+        module->Process();
+    }
+    catch (runtime_error &e)
+    {
+        cout << "High-order surface meshing has failed with message:" << endl;
+        cout << e.what() << endl;
+        cout << "The mesh will be written as normal but the incomplete surface "
+                "will remain faceted"
+             << endl;
+        return;
     }
 
     ////*** VARIATIONAL OPTIMISATION ****////
-    if(m_varopti)
+    if (m_varopti)
     {
         unsigned int np = boost::thread::physical_concurrency();
-        if(m_mesh->m_verbose)
+        if (m_mesh->m_verbose)
+        {
+            cout << "Detecting 4 cores, will attempt to run in parrallel"
+                 << endl;
+        }
+        module = GetModuleFactory().CreateInstance(
+            ModuleKey(eProcessModule, "varopti"), m_mesh);
+        module->RegisterConfig("hyperelastic", "");
+        module->RegisterConfig("maxiter", "10");
+        module->RegisterConfig("numthreads",
+                                    boost::lexical_cast<string>(np));
+
+        try
+        {
+            module->SetDefaults();
+            module->Process();
+        }
+        catch (runtime_error &e)
         {
-            cout << "Detecting 4 cores, will attempt to run in parrallel" << endl;
+            cout << "Variational optimisation has failed with message:" << endl;
+            cout << e.what() << endl;
+            cout << "The mesh will be written as is, it may be invalid" << endl;
+            return;
         }
-        mods.push_back(GetModuleFactory().CreateInstance(
-            ModuleKey(eProcessModule, "varopti"), m_mesh));
-        mods.back()->RegisterConfig("hyperelastic","");
-        mods.back()->RegisterConfig("maxiter","10");
-        mods.back()->RegisterConfig("numthreads",boost::lexical_cast<string>(np));
     }
 
     ////**** SPLIT BL ****////
-    if(m_splitBL)
+    if (m_splitBL)
     {
-        mods.push_back(GetModuleFactory().CreateInstance(
-            ModuleKey(eProcessModule, "bl"), m_mesh));
-        mods.back()->RegisterConfig("layers",m_bllayers);
-        mods.back()->RegisterConfig("surf",m_blsurfs);
-        mods.back()->RegisterConfig("nq",boost::lexical_cast<string>(m_mesh->m_nummode));
-        mods.back()->RegisterConfig("r",m_blprog);
+        module = GetModuleFactory().CreateInstance(
+            ModuleKey(eProcessModule, "bl"), m_mesh);
+        module->RegisterConfig("layers", m_bllayers);
+        module->RegisterConfig("surf", m_blsurfs);
+        module->RegisterConfig(
+            "nq", boost::lexical_cast<string>(m_mesh->m_nummode));
+        module->RegisterConfig("r", m_blprog);
+
+        try
+        {
+            module->SetDefaults();
+            module->Process();
+        }
+        catch (runtime_error &e)
+        {
+            cout << "Boundary layer splitting has failed with message:" << endl;
+            cout << e.what() << endl;
+            cout << "The mesh will be written as is, it may be invalid" << endl;
+            return;
+        }
     }
 
     ////**** Peralign ****////
@@ -426,20 +526,16 @@ void InputMCF::Process()
         for (vector<string>::iterator il = lines.begin(); il != lines.end();
              ++il)
         {
-            mods.push_back(GetModuleFactory().CreateInstance(
-                ModuleKey(eProcessModule, "peralign"), m_mesh));
-            
+            module = GetModuleFactory().CreateInstance(
+                ModuleKey(eProcessModule, "peralign"), m_mesh);
+
             vector<string> tmp(2);
             boost::split(tmp, *il, boost::is_any_of(","));
-            mods.back()->RegisterConfig("surf1", tmp[0]);
-            mods.back()->RegisterConfig("surf2", tmp[1]);
+            module->RegisterConfig("surf1", tmp[0]);
         }
-    }
 
-    for(int i = 0; i < mods.size(); i++)
-    {
-        mods[i]->SetDefaults();
-        mods[i]->Process();
+        module->SetDefaults();
+        module->Process();
     }
 }
 }

utilities/NekMesh/InputModules/InputMCF.h

@@ -66,7 +66,7 @@ private:
                 m_nacadomain, m_periodic, m_adjustment;
 
     bool m_makeBL, m_surfopti, m_varopti, m_refine, m_woct, m_2D, m_splitBL,
-         m_naca, m_adjust, m_adjustall;
+         m_naca, m_adjust, m_adjustall, m_manifold;
 
 };
 }

utilities/NekMesh/Tests/MeshGen/cylinder.mcf

@@ -7,9 +7,9 @@
         </INFORMATION>
 
         <PARAMETERS>
-            <P PARAM="MinDelta" VALUE="0.003"    />
+            <P PARAM="MinDelta" VALUE="0.03"    />
             <P PARAM="MaxDelta" VALUE="0.2"     />
-            <P PARAM="EPS"      VALUE="0.05"     />
+            <P PARAM="EPS"      VALUE="0.5"     />
 
             <P PARAM="Order"    VALUE="4"       />
 

utilities/NekMesh/Tests/MeshGen/rev.mcf

@@ -7,9 +7,9 @@
         </INFORMATION>
 
         <PARAMETERS>
-            <P PARAM="MinDelta" VALUE="0.012"    />
+            <P PARAM="MinDelta" VALUE="0.02"    />
             <P PARAM="MaxDelta" VALUE="0.5"     />
-            <P PARAM="EPS"      VALUE="0.005"     />
+            <P PARAM="EPS"      VALUE="0.05"     />
 
             <P PARAM="Order"    VALUE="4"       />
 

utilities/NekMesh/Tests/MeshGen/sphere.mcf

@@ -7,9 +7,9 @@
         </INFORMATION>
 
         <PARAMETERS>
-            <P PARAM="MinDelta" VALUE="0.01"    />
+            <P PARAM="MinDelta" VALUE="0.05"    />
             <P PARAM="MaxDelta" VALUE="0.2"     />
-            <P PARAM="EPS"      VALUE="0.01"     />
+            <P PARAM="EPS"      VALUE="0.03"     />
 
             <P PARAM="Order"    VALUE="4"       />
         </PARAMETERS>