Add fix for triangles with Lagrange basis

CHANGELOG.md

@@ -8,6 +8,7 @@ v5.8.0
 - Add FieldConvert module to perform local stability analysis for compressible flows (!1319)
 - Remove get() accessor from Array data structure (!1937)
 - Fix issue with `StdTetExp::v_LocCollapsedToLocCoord` (!1946)
+- Fix issue with `NodalTriExp::v_GetTracePhysVals` (!1951)
 
 **CI**
 -- fix CubeAllElements performance test tolerance (!1943)

library/LocalRegions/NodalTriExp.cpp

@@ -683,4 +683,64 @@ void NodalTriExp::v_ComputeTraceNormal(const int edge)
     }
 }
 
+void NodalTriExp::v_GetTracePhysVals(
+    const int edge, const StdRegions::StdExpansionSharedPtr &EdgeExp,
+    const Array<OneD, const NekDouble> &inarray,
+    Array<OneD, NekDouble> &outarray, StdRegions::Orientation orient)
+{
+    int nquad0 = m_base[0]->GetNumPoints();
+    int nquad1 = m_base[1]->GetNumPoints();
+    int nt     = 0;
+    // Extract in Cartesian direction because we have to deal with
+    // e.g. Gauss-Radau points.
+    switch (edge)
+    {
+        case 0:
+            Vmath::Vcopy(nquad0, &(inarray[0]), 1, &(outarray[0]), 1);
+            nt = nquad0;
+            break;
+        case 1:
+            Vmath::Vcopy(nquad1, &(inarray[0]) + (nquad0 - 1), nquad0,
+                         &(outarray[0]), 1);
+            nt = nquad1;
+            break;
+        case 2:
+            Vmath::Vcopy(nquad1, &(inarray[0]), nquad0, &(outarray[0]), 1);
+            nt = nquad1;
+            break;
+        default:
+            ASSERTL0(false, "edge value (< 3) is out of range");
+            break;
+    }
+
+    ASSERTL1(EdgeExp->GetBasis(0)->GetPointsType() ==
+                 LibUtilities::eGaussLobattoLegendre,
+             "Edge expansion should be GLL");
+
+    // Interpolate if required
+    if (m_base[edge ? 1 : 0]->GetPointsKey() !=
+        EdgeExp->GetBasis(0)->GetPointsKey())
+    {
+        Array<OneD, NekDouble> outtmp(max(nquad0, nquad1));
+
+        Vmath::Vcopy(nt, outarray, 1, outtmp, 1);
+
+        LibUtilities::Interp1D(m_base[edge ? 1 : 0]->GetPointsKey(), outtmp,
+                               EdgeExp->GetBasis(0)->GetPointsKey(), outarray);
+    }
+
+    if (orient == StdRegions::eNoOrientation)
+    {
+        orient = GetTraceOrient(edge);
+    }
+
+    // Reverse data if necessary
+    if (orient == StdRegions::eBackwards)
+    {
+        Vmath::Reverse(EdgeExp->GetNumPoints(0), &outarray[0], 1, &outarray[0],
+                       1);
+    }
+}
+
+
 } // namespace Nektar::LocalRegions

library/LocalRegions/NodalTriExp.h

@@ -369,6 +369,12 @@ private:
     }
 
     void v_ComputeTraceNormal(const int edge) override;
+
+    void v_GetTracePhysVals(
+        const int edge, const StdRegions::StdExpansionSharedPtr &EdgeExp,
+        const Array<OneD, const NekDouble> &inarray,
+        Array<OneD, NekDouble> &outarray, StdRegions::Orientation orient)
+        override;
 };
 
 typedef std::shared_ptr<NodalTriExp> NodalTriExpSharedPtr;

solvers/IncNavierStokesSolver/CMakeLists.txt

@@ -119,6 +119,7 @@ IF( NEKTAR_SOLVER_INCNAVIERSTOKES )
     ADD_NEKTAR_TEST(ChanFlow_m8_BodyForce)
     ADD_NEKTAR_TEST(ChanFlow_m8_singular)
     ADD_NEKTAR_TEST(ChanFlow_V8P7_Avg)
+    ADD_NEKTAR_TEST(ChanFlow_Tri_Lagrange)
     ADD_NEKTAR_TEST(Channel_Flow_3modes_rad)
     ADD_NEKTAR_TEST(channelTemp)
     ADD_NEKTAR_TEST(channelFlow_Heating_FlowRate)

solvers/IncNavierStokesSolver/Tests/ChanFlow_Tri_Lagrange.tst

+<?xml version="1.0" encoding="utf-8" ?>
+<test>
+    <description>Channel Flow P=8</description>
+    <executable>IncNavierStokesSolver</executable>
+    <parameters>ChanFlow_Tri_Lagrange.xml</parameters>
+    <files>
+        <file description="Session File">ChanFlow_Tri_Lagrange.xml</file>
+    </files>
+    <metrics>
+        <metric type="L2" id="1">
+            <value variable="u" tolerance="1e-12">2.44199e-12</value>
+            <value variable="v" tolerance="1e-12">4.18302e-13</value>
+            <value variable="p" tolerance="1e-8">3.31598e-11</value>
+        </metric>
+        <metric type="Linf" id="2">
+            <value variable="u" tolerance="1e-12">1.24133e-11</value>
+            <value variable="v" tolerance="1e-12">2.05103e-12</value>
+            <value variable="p" tolerance="1e-8">5.98182e-10</value>
+        </metric>
+    </metrics>
+</test>

solvers/IncNavierStokesSolver/Tests/ChanFlow_Tri_Lagrange.xml

+<?xml version="1.0" encoding="utf-8" ?>
+<NEKTAR>
+    <GEOMETRY DIM="2" SPACE="2">
+        <VERTEX COMPRESSED="B64Z-LittleEndian" BITSIZE="64">eJxjYMAPGKF0KBhctUeXZ0KTR1fHjKYeXZ4Fq/6ncHlWrPoR8mxYzf0Al2fH6iuEPAcOc2GAE00e3f1cWOUR5nBjNR9hPw8Od8EAL5o8uv18WOUR9vNjNR9hDwDiLzuR</VERTEX>
+        <EDGE COMPRESSED="B64Z-LittleEndian" BITSIZE="64">eJx1kMsSgjAQBIOCIgoqCigv+f+v9JI+0FXJZaunZmFnQti/LDEPick7SodzMa+QH98psXeWzn4pxn9JfLcSs3eVzv9uYvx1nPTDHY2YO+5i9h7S6eMpnbtbMf6XdHK/pZOzE+PvpdPHECe56OUjppevmL1ROvkm6fQ4i/Ev0sm3Sqf3nxj/Jp3cfzdUA+4A</EDGE>
+        <ELEMENT>
+            <T COMPRESSED="B64Z-LittleEndian" BITSIZE="64">eJx1z8EOgjAQhGFAQEALVVFQUHj/p/Tyz4FJ2suXmTSb3Sw7vhyLRD5hmegrrK1XPmNjc9S32Nk85QtebW7AHgebqxxRdzXW3/Bu+yg/EnupH/Fp+yq/bK7umHDGNwbLH+xN9QuuOFj+ou6M1v9ws3/KO+rOPw99A+wA</T>
+        </ELEMENT>
+        <COMPOSITE>
+            <C ID="0"> T[0-17] </C>
+            <C ID="1"> E[0,13,23] </C>
+            <C ID="2"> E[24,27,30] </C>
+            <C ID="3"> E[11,22,32] </C>
+            <C ID="4"> E[4,8,12] </C>
+        </COMPOSITE>
+        <DOMAIN>
+            <D ID="0"> C[0] </D>
+        </DOMAIN>
+    </GEOMETRY>
+    <EXPANSIONS>
+        <E COMPOSITE="C[0]" NUMMODES="7" FIELDS="u,v,p" TYPE="GLL_LAGRANGE" />
+    </EXPANSIONS>
+    <CONDITIONS>
+        <SOLVERINFO>
+            <I PROPERTY="SolverType" VALUE="VelocityCorrectionScheme" />
+            <I PROPERTY="EQTYPE" VALUE="UnsteadyNavierStokes" />
+            <I PROPERTY="AdvectionForm" VALUE="Convective" />
+            <I PROPERTY="Projection" VALUE="Galerkin" />
+            <I PROPERTY="TimeIntegrationMethod" VALUE="IMEXOrder1" />
+        </SOLVERINFO>
+
+        <PARAMETERS>
+            <P> TimeStep = 0.001     </P>
+            <P> NumSteps = 1000       </P>
+            <P> IO_CheckSteps = 1000       </P>
+            <P> IO_InfoSteps = 1000       </P>
+            <P> Kinvis = 1         </P>
+        </PARAMETERS>
+
+        <VARIABLES>
+            <V ID="0"> u </V>
+            <V ID="1"> v </V>
+            <V ID="2"> p </V>
+        </VARIABLES>
+
+        <BOUNDARYREGIONS>
+            <B ID="0"> C[1,3] </B>
+            <B ID="1"> C[4] </B>
+            <B ID="2"> C[2] </B>
+        </BOUNDARYREGIONS>
+
+        <BOUNDARYCONDITIONS>
+            <REGION REF="0">
+                <D VAR="u" VALUE="0" />
+                <D VAR="v" VALUE="0" />
+                <N VAR="p" USERDEFINEDTYPE="H" VALUE="0" />
+            </REGION>
+            <REGION REF="1">
+                <D VAR="u" VALUE="y*(1-y)" />
+                <D VAR="v" VALUE="0" />
+                <N VAR="p" USERDEFINEDTYPE="H" VALUE="0" />
+            </REGION>
+            <REGION REF="2">
+                <N VAR="u" VALUE="0" />
+                <N VAR="v" VALUE="0" />
+                <D VAR="p" VALUE="0" />
+            </REGION>
+        </BOUNDARYCONDITIONS>
+
+        <FUNCTION NAME="InitialConditions">
+            <E VAR="u" VALUE="0" />
+            <E VAR="v" VALUE="0" />
+            <E VAR="p" VALUE="0" />
+        </FUNCTION>
+
+        <FUNCTION NAME="ExactSolution">
+            <E VAR="u" VALUE="y*(1-y)" />
+            <E VAR="v" VALUE="0" />
+            <E VAR="p" VALUE="-2*Kinvis*(x-1)" />
+        </FUNCTION>
+    </CONDITIONS>
+    <Metadata>
+        <Provenance>
+            <Hostname>aegir.lan</Hostname>
+            <NektarVersion>5.7.0</NektarVersion>
+            <Timestamp>15-Nov-2024 19:13:25</Timestamp>
+        </Provenance>
+        <NekMeshCommandLine>square-tris.msh ChanFlow_Lagrange_Tri.xml </NekMeshCommandLine>
+    </Metadata>
+</NEKTAR>