minor bug fixing

library/FieldUtils/ProcessModules/ProcessQualityMetric.cpp

@@ -45,6 +45,7 @@ using namespace std;
 #include <StdRegions/StdPrismExp.h>
 #include <StdRegions/StdQuadExp.h>
 #include <StdRegions/StdTetExp.h>
+#include <StdRegions/StdHexExp.h>
 #include <StdRegions/StdTriExp.h>
 
 namespace Nektar
@@ -295,6 +296,82 @@ inline vector<DNekMat> MappingIdealToRef(SpatialDomains::GeometrySharedPtr geom,
             }
         }
     }
+    else if(geom->GetShapeType() == LibUtilities::eHexahedron)
+    {
+        vector<Array<OneD, NekDouble> > xyz;
+        for (int i = 0; i < geom->GetNumVerts(); i++)
+        {
+            Array<OneD, NekDouble> loc(3);
+            SpatialDomains::PointGeomSharedPtr p = geom->GetVertex(i);
+            p->GetCoords(loc);
+            xyz.push_back(loc);
+        }
+
+        Array<OneD, const LibUtilities::BasisSharedPtr> b = chi->GetBase();
+        Array<OneD, NekDouble> eta1                       = b[0]->GetZ();
+        Array<OneD, NekDouble> eta2                       = b[1]->GetZ();
+        Array<OneD, NekDouble> eta3                       = b[2]->GetZ();
+
+        for (int k = 0; k < b[2]->GetNumPoints(); k++)
+        {
+            for (int j = 0; j < b[1]->GetNumPoints(); j++)
+            {
+                for (int i = 0; i < b[0]->GetNumPoints(); i++)
+                {
+                    NekDouble a1  = 0.5 * (1 - eta1[i]);
+                    NekDouble a2  = 0.5 * (1 + eta1[i]);
+                    NekDouble b1  = 0.5 * (1 - eta2[j]),
+                              b2  = 0.5 * (1 + eta2[j]);
+                    NekDouble c1  = 0.5 * (1 - eta3[k]),
+                              c2  = 0.5 * (1 + eta3[k]);
+
+                    DNekMat dxdz(3, 3, 1.0, eFULL);
+
+                    dxdz(0,0) = - 0.5 * b1 * c1 * xyz[0][0] + 0.5 * b1 * c1 * xyz[1][0]
+                             + 0.5 * b2 * c1 * xyz[2][0] - 0.5 * b2 * c1 * xyz[3][0]
+                             - 0.5 * b1 * c2 * xyz[5][0] + 0.5 * b1 * c2 * xyz[5][0]
+                             + 0.5 * b2 * c2 * xyz[6][0] - 0.5 * b2 * c2 * xyz[7][0];
+                    dxdz(1,0) = - 0.5 * b1 * c1 * xyz[0][1] + 0.5 * b1 * c1 * xyz[1][1]
+                             + 0.5 * b2 * c1 * xyz[2][1] - 0.5 * b2 * c1 * xyz[3][1]
+                             - 0.5 * b1 * c2 * xyz[5][1] + 0.5 * b1 * c2 * xyz[5][1]
+                             + 0.5 * b2 * c2 * xyz[6][1] - 0.5 * b2 * c2 * xyz[7][1];
+                    dxdz(2,0) = - 0.5 * b1 * c1 * xyz[0][2] + 0.5 * b1 * c1 * xyz[1][2]
+                             + 0.5 * b2 * c1 * xyz[2][2] - 0.5 * b2 * c1 * xyz[3][2]
+                             - 0.5 * b1 * c2 * xyz[5][2] + 0.5 * b1 * c2 * xyz[5][2]
+                             + 0.5 * b2 * c2 * xyz[6][2] - 0.5 * b2 * c2 * xyz[7][2];
+
+                    dxdz(0,1) = - 0.5 * a1 * c1 * xyz[0][0] - 0.5 * a2 * c1 * xyz[1][0]
+                             + 0.5 * a2 * c1 * xyz[2][0] + 0.5 * a1 * c1 * xyz[3][0]
+                             - 0.5 * a1 * c2 * xyz[5][0] - 0.5 * a2 * c2 * xyz[5][0]
+                             + 0.5 * a2 * c2 * xyz[6][0] + 0.5 * a1 * c2 * xyz[7][0];
+                    dxdz(1,1) = - 0.5 * a1 * c1 * xyz[0][1] - 0.5 * a2 * c1 * xyz[1][1]
+                             + 0.5 * a2 * c1 * xyz[2][1] + 0.5 * a1 * c1 * xyz[3][1]
+                             - 0.5 * a1 * c2 * xyz[5][1] - 0.5 * a2 * c2 * xyz[5][1]
+                             + 0.5 * a2 * c2 * xyz[6][1] + 0.5 * a1 * c2 * xyz[7][1];
+                    dxdz(2,1) = - 0.5 * a1 * c1 * xyz[0][2] - 0.5 * a2 * c1 * xyz[1][2]
+                             + 0.5 * a2 * c1 * xyz[2][2] + 0.5 * a1 * c1 * xyz[3][2]
+                             - 0.5 * a1 * c2 * xyz[5][2] - 0.5 * a2 * c2 * xyz[5][2]
+                             + 0.5 * a2 * c2 * xyz[6][2] + 0.5 * a1 * c2 * xyz[7][2];
+
+                    dxdz(0,0) = - 0.5 * b1 * a1 * xyz[0][0] - 0.5 * b1 * a2 * xyz[1][0]
+                             - 0.5 * b2 * a2 * xyz[2][0] - 0.5 * b2 * a1 * xyz[3][0]
+                             + 0.5 * b1 * a1 * xyz[5][0] + 0.5 * b1 * a2 * xyz[5][0]
+                             + 0.5 * b2 * a2 * xyz[6][0] + 0.5 * b2 * a1 * xyz[7][0];
+                    dxdz(1,0) = - 0.5 * b1 * a1 * xyz[0][1] - 0.5 * b1 * a2 * xyz[1][1]
+                             - 0.5 * b2 * a2 * xyz[2][1] - 0.5 * b2 * a1 * xyz[3][1]
+                             + 0.5 * b1 * a1 * xyz[5][1] + 0.5 * b1 * a2 * xyz[5][1]
+                             + 0.5 * b2 * a2 * xyz[6][1] + 0.5 * b2 * a1 * xyz[7][1];
+                    dxdz(2,0) = - 0.5 * b1 * a1 * xyz[0][2] - 0.5 * b1 * a2 * xyz[1][2]
+                             - 0.5 * b2 * a2 * xyz[2][2] - 0.5 * b2 * a1 * xyz[3][2]
+                             + 0.5 * b1 * a1 * xyz[5][2] + 0.5 * b1 * a2 * xyz[5][2]
+                             + 0.5 * b2 * a2 * xyz[6][2] + 0.5 * b2 * a1 * xyz[7][2];
+
+                    dxdz.Invert();
+                    ret.push_back(dxdz);
+                }
+            }
+        }
+    }
     else
     {
         ASSERTL0(false, "not coded");
@@ -359,6 +436,10 @@ Array<OneD, NekDouble> ProcessQualityMetric::GetQ(LocalRegions::ExpansionSharedP
             chiMod = MemoryManager<StdRegions::StdPrismExp>::AllocateSharedPtr(
                 basisKeys[0], basisKeys[1], basisKeys[2]);
             break;
+        case LibUtilities::eHexahedron:
+            chiMod = MemoryManager<StdRegions::StdHexExp>::AllocateSharedPtr(
+                basisKeys[0], basisKeys[1], basisKeys[2]);
+            break;
         default:
             ASSERTL0(false, "nope");
     }

library/LibUtilities/Foundations/NodalHexElec.cpp

@@ -64,9 +64,9 @@ void NodalHexElec::CalculatePoints()
         {
             for (int i = 0; i < numPoints; i++, ct++)
             {
-                m_points[0][ct] = m_e0[k];
+                m_points[0][ct] = m_e0[i];
                 m_points[1][ct] = m_e0[j];
-                m_points[2][ct] = m_e0[i];
+                m_points[2][ct] = m_e0[k];
             }
         }
     }

library/LibUtilities/Foundations/NodalQuadElec.cpp

@@ -62,8 +62,8 @@ void NodalQuadElec::CalculatePoints()
     {
         for (int i = 0; i < numPoints; i++, ct++)
         {
-            m_points[0][ct] = m_e0[j];
-            m_points[1][ct] = m_e0[i];
+            m_points[0][ct] = m_e0[i];
+            m_points[1][ct] = m_e0[j];
         }
     }
 }

library/LibUtilities/Foundations/NodalUtil.cpp

@@ -967,11 +967,11 @@ NodalUtilHex::NodalUtilHex(int degree,
 
     // Construct a mapping (i,j,k) -> m from the tensor product space (i,j,k) to
     // a single ordering m.
-    for (int i = 0; i <= m_degree; ++i)
+    for (int k = 0; k <= m_degree; ++k)
     {
         for (int j = 0; j <= m_degree; ++j)
         {
-            for (int k = 0; k <= m_degree - i; ++k)
+            for (int i = 0; i <= m_degree; ++i)
             {
                 m_ordering.push_back(Mode(i, j, k));
             }
@@ -1053,11 +1053,18 @@ NekVector<NekDouble> NodalUtilHex::v_OrthoBasisDeriv(
             ret[i] = jacobi_di[i] * jacobi_j[i] * jacobi_k[i];
         }
     }
+    else if (dir == 1)
+    {
+        for (int i = 0; i < m_numPoints; ++i)
+        {
+            ret[i] = jacobi_dj[i] * jacobi_i[i] * jacobi_k[i];
+        }
+    }
     else
     {
         for (int i = 0; i < m_numPoints; ++i)
         {
-            ret[i] = jacobi_i[i] * jacobi_dj[i] * jacobi_dk[i];
+            ret[i] = jacobi_i[i] * jacobi_j[i] * jacobi_dk[i];
         }
     }
 

utilities/NekMesh/ProcessModules/ProcessVarOpti/ElUtil.cpp

@@ -638,7 +638,7 @@ void ElUtil::Evaluate()
     m_res->worstJac = min(m_res->worstJac,(mn2 / mx2));
     mtx2.unlock();
 
-    m_scaledJac = (mn2/mx2) ;
+    m_scaledJac = (mn2/mx2);
 }
 
 void ElUtil::InitialMinJac()

utilities/NekMesh/ProcessModules/ProcessVarOpti/Evaluator.hxx

@@ -287,10 +287,10 @@ NekDouble NodeOpti::GetFunctional(NekDouble &minJacNew, bool gradient)
 
     minJacNew = std::numeric_limits<double>::max();
     NekDouble integral = 0.0;
-    //NekDouble ep = minJac < 0.0 ? sqrt(1e-12 + 0.04*minJac*minJac) : 1e-6;
-    //NekDouble ep = minJac < 0.0 ? sqrt(gam*(gam-minJac)) : gam;
+    NekDouble ep = m_minJac < 0.0 ? sqrt(1e-8 + 0.04*m_minJac*m_minJac) : 1e-4;
+    //NekDouble ep = m_minJac < 0.0 ? sqrt(gam*(gam-m_minJac)) : gam;
     //NekDouble ep = minJac < 0.0 ? sqrt(gam*gam + minJac*minJac) : gam;
-    NekDouble ep = 1e-2;
+    //NekDouble ep = 1e-2;
     NekDouble jacIdeal[DIM][DIM], jacDet;
     m_grad = Array<OneD, NekDouble>(DIM == 2 ? 5 : 9, 0.0);
 

utilities/NekMesh/ProcessModules/ProcessVarOpti/NodeOpti.cpp

@@ -55,7 +55,7 @@ NodeOptiFactory &GetNodeOptiFactory()
     return asd;
 }
 
-const NekDouble NodeOpti::gam = numeric_limits<float>::epsilon();
+const NekDouble NodeOpti::gam = numeric_limits<float>::epsilon()*10;
 
 void NodeOpti::CalcMinJac()
 {

utilities/NekMesh/ProcessModules/ProcessVarOpti/PreProcessing.cpp

@@ -52,10 +52,14 @@ map<LibUtilities::ShapeType, DerivUtilSharedPtr> ProcessVarOpti::BuildDerivUtil(
     map<LibUtilities::ShapeType, DerivUtilSharedPtr> ret;
 
     map<LibUtilities::ShapeType, LibUtilities::PointsType> typeMap;
-    typeMap[LibUtilities::eTriangle] = LibUtilities::eNodalTriSPI;
+    if(m_mesh->m_nummode + o <= 11)
+    {
+        typeMap[LibUtilities::eTriangle] = LibUtilities::eNodalTriSPI;
+        typeMap[LibUtilities::eTetrahedron] = LibUtilities::eNodalTetSPI;
+        typeMap[LibUtilities::ePrism] = LibUtilities::eNodalPrismSPI;
+    }
+
     typeMap[LibUtilities::eQuadrilateral] = LibUtilities::eNodalQuadElec;
-    typeMap[LibUtilities::eTetrahedron] = LibUtilities::eNodalTetSPI;
-    typeMap[LibUtilities::ePrism] = LibUtilities::eNodalPrismSPI;
     //typeMap[LibUtilities::eHexahedron] = LibUtilities::eNodalHexElec;
 
     map<LibUtilities::ShapeType, LibUtilities::PointsType> typeMap2;

utilities/NekMesh/ProcessModules/ProcessVarOpti/ProcessVarOpti.cpp

@@ -262,8 +262,6 @@ void ProcessVarOpti::Process()
          << "Max set:\t\t" << mx << endl
          << "Residual tolerance:\t" << restol << endl;
 
-    //return;
-
     int nThreads = m_config["numthreads"].as<int>();
 
     int ctr = 0;