Fixed history point filter to select correct element in parallel manifold case.

library/SolverUtils/Filters/FilterHistoryPoints.cpp

@@ -155,23 +155,50 @@ namespace Nektar
             int vRank = vComm->GetRank();
             Array<OneD, int>  procList(m_historyPoints.size(), -1);
             Array<OneD, int> idList(m_historyPoints.size());
+            Array<OneD, NekDouble>  locCoords(3);
             std::vector<Array<OneD, NekDouble> > LocCoords; 
             
             for (i = 0; i < m_historyPoints.size(); ++i)
             {
-                Array<OneD, NekDouble>  locCoords(3);
 
+                // Determine the expansion and local coordinates
                 m_historyPoints[i]->GetCoords(  gloCoord[0],
                                                 gloCoord[1],
                                                 gloCoord[2]);
 
-                idList[i] = pFields[0]->GetExpIndex(gloCoord,locCoords,NekConstants::kGeomFactorsTol);
+                idList[i] = pFields[0]->GetExpIndex(gloCoord,locCoords,
+                                                NekConstants::kGeomFactorsTol);
+
+                // Check if the reverse mapping of the local coordinates gives
+                // the correct coordinates of the history point. This ensures
+                // that the correct element is chosen in the manifold case.
+                if (idList[i] != -1)
+                {
+                    SpatialDomains::GeometrySharedPtr g =
+                                    pFields[0]->GetExp(idList[i])->GetGeom();
+                    StdRegions::StdExpansionSharedPtr e = g->GetXmap();
+                    Array<OneD, NekDouble> xvals(e->GetTotPoints());
+                    Array<OneD, NekDouble> yvals(e->GetTotPoints());
+                    Array<OneD, NekDouble> zvals(e->GetTotPoints());
+                    e->BwdTrans(g->GetCoeffs(0), xvals);
+                    e->BwdTrans(g->GetCoeffs(1), yvals);
+                    e->BwdTrans(g->GetCoeffs(2), zvals);
+                    double x = e->PhysEvaluate(locCoords, xvals) - gloCoord[0];
+                    double y = e->PhysEvaluate(locCoords, yvals) - gloCoord[1];
+                    double z = e->PhysEvaluate(locCoords, zvals) - gloCoord[2];
+                    if (x*x + y*y + z*z > NekConstants::kGeomFactorsTol)
+                    {
+                        idList[i] = -1;
+                    }
+                }
 
                 // Save Local coordinates for later
                 LocCoords.push_back(locCoords);
-                // Set element id to Vid of m_historyPoints;
+                // Set element id to Vid of m_history point for later use
                 m_historyPoints[i]->SetVid(idList[i]);
                 
+                // If a matching element is found on this process, note the
+                // process ID
                 if (idList[i] != -1) 
                 {
                     if(m_isHomogeneous1D)
@@ -198,6 +225,9 @@ namespace Nektar
                     }
                 }
             }
+
+            // Reduce process IDs for all history points. The process with
+            // largest rank will handle the history point
             vComm->AllReduce(procList, LibUtilities::ReduceMax);
 
             // Determine the element in which each history point resides.
@@ -205,11 +235,14 @@ namespace Nektar
             for (i = 0; i < m_historyPoints.size(); ++i)
             {
                 // If point lies on partition boundary, only the proc with max
-                // rank retains posession.
+                // rank retains possession.
                 if (procList[i] != vRank)
                 {
                     idList[i] = -1;
                 }
+
+                // If the current process owns this history point, add it to its
+                // local list of history points.
                 if (idList[i] != -1)
                 {
                     m_historyLocalPointMap[m_historyList.size()] = i;