From 7d730126a48fe759c66addd4a33a0080ff21f0c0 Mon Sep 17 00:00:00 2001
From: Mauro Perego <mperego@sandia.gov>
Date: Thu, 3 Sep 2015 09:50:11 -0600
Subject: [PATCH] add function procsSharingVertex to dycore interface

This function computes the ranks of processes that share a vertex.
This is called by Albany-Felix starting with Albany hash:
39ed4d5f0c4fb790622cfeaa7943806afd0ef9e1 (Sept. 3, 2015)
to optimize the construction of the FE grid.

(Albany builds newer than that commit will require this change to MPAS.
---
 .../Interface_velocity_solver.cpp             | 33 ++++++++++++++++---
 .../Interface_velocity_solver.hpp             |  1 +
 2 files changed, 30 insertions(+), 4 deletions(-)

diff --git a/src/core_landice/mode_forward/Interface_velocity_solver.cpp b/src/core_landice/mode_forward/Interface_velocity_solver.cpp
index d086bde09ef1..ef1dff7b551d 100644
--- a/src/core_landice/mode_forward/Interface_velocity_solver.cpp
+++ b/src/core_landice/mode_forward/Interface_velocity_solver.cpp
@@ -29,7 +29,7 @@ int nVertices, nEdges, nTriangles, nGlobalVertices, nGlobalEdges,
 int maxNEdgesOnCell_F;
 int const *cellsOnEdge_F, *cellsOnVertex_F, *verticesOnCell_F,
     *verticesOnEdge_F, *edgesOnCell_F, *indexToCellID_F, *nEdgesOnCells_F,
-    *dirichletCellsMask_F, *floatingEdgesMask_F;
+    *dirichletCellsMask_F, *floatingEdgesMask_F, *verticesMask_F;
 std::vector<double> layersRatio, levelsNormalizedThickness;
 int nLayers;
 double const *xCell_F, *yCell_F, *zCell_F, *xVertex_F,  *yVertex_F, *zVertex_F, *areaTriangle_F;
@@ -47,7 +47,8 @@ int ice_present_bit_value;
 
 //void *phgGrid = 0;
 std::vector<int> edgesToReceive, fCellsToReceive, indexToTriangleID,
-    verticesOnTria, trianglesOnEdge, trianglesPositionsOnEdge, verticesOnEdge;
+    verticesOnTria, trianglesOnEdge, trianglesPositionsOnEdge, verticesOnEdge,
+    trianglesProcIds,  reduced_ranks;
 std::vector<int> indexToVertexID, vertexToFCell, indexToEdgeID, edgeToFEdge,
     mask, fVertexToTriangleID, fCellToVertex, floatingEdgesIds, dirichletNodesIDs;
 std::vector<double> temperatureOnTetra, velocityOnVertices, velocityOnCells,
@@ -160,6 +161,9 @@ void velocity_solver_set_grid_data(int const* _nCells_F, int const* _nEdges_F,
   recvVerticesList_F = new exchangeList_Type(
       unpackMpiArray(recvVerticesArray_F));
 
+  trianglesProcIds.resize(nVertices_F);
+  getProcIds(trianglesProcIds, recvVerticesList_F);
+
   if (radius > 10) {
     xCellProjected.resize(nCells_F);
     yCellProjected.resize(nCells_F);
@@ -491,9 +495,10 @@ void velocity_solver_finalize() {
  *
  */
 
-void velocity_solver_compute_2d_grid(int const* verticesMask_F, int const* _dirichletCellsMask_F, int const* _floatingEdgesMask_F) {
+void velocity_solver_compute_2d_grid(int const* _verticesMask_F, int const* _dirichletCellsMask_F, int const* _floatingEdgesMask_F) {
   int numProcs, me;
 
+  verticesMask_F = _verticesMask_F;
   dirichletCellsMask_F = _dirichletCellsMask_F;
   floatingEdgesMask_F = _floatingEdgesMask_F;
 
@@ -1379,9 +1384,12 @@ void createReducedMPI(int nLocalEntities, MPI_Comm& reduced_comm_id) {
   int nonEmpty = int(nLocalEntities > 0);
   MPI_Allgather(&nonEmpty, 1, MPI_INT, &haveElements[0], 1, MPI_INT, comm);
   std::vector<int> ranks;
+  reduced_ranks.resize(numProcs,0);
   for (int i = 0; i < numProcs; i++) {
-    if (haveElements[i])
+    if (haveElements[i]) {
+      reduced_ranks[i] = ranks.size();
       ranks.push_back(i);
+    }
   }
 
   MPI_Comm_group(comm, &world_group_id);
@@ -1783,3 +1791,20 @@ int prismType(long long int const* prismVertexMpasIds, int& minIndex)
     }
   }
 
+  void procsSharingVertex(const int vertex, std::vector<int>& procIds) {
+    int fCell = vertexToFCell[vertex];
+    procIds.clear();
+    int nEdg = nEdgesOnCells_F[fCell];
+    int me;
+    MPI_Comm_rank(comm, &me);
+    procIds.reserve(nEdg);
+    for(int i=0; i<nEdg; ++i) {
+      int fVertex = verticesOnCell_F[maxNEdgesOnCell_F * fCell + i]-1;
+      if (verticesMask_F[fVertex] & dynamic_ice_bit_value) {
+        int proc = trianglesProcIds[fVertex];
+        if(proc != me)
+          procIds.push_back(reduced_ranks[proc]);
+      }
+    }
+  }
+
diff --git a/src/core_landice/mode_forward/Interface_velocity_solver.hpp b/src/core_landice/mode_forward/Interface_velocity_solver.hpp
index b3ae404bdd0b..5441d3130212 100644
--- a/src/core_landice/mode_forward/Interface_velocity_solver.hpp
+++ b/src/core_landice/mode_forward/Interface_velocity_solver.hpp
@@ -282,6 +282,7 @@ void computeMap();
 
 void setBdFacesOnPrism (const std::vector<std::vector<std::vector<int> > >& prismStruct, const std::vector<int>& prismFaceIds, std::vector<int>& tetraPos, std::vector<int>& facePos);
 void tetrasFromPrismStructured (int const* prismVertexMpasIds, int const* prismVertexGIds, int tetrasIdsOnPrism[][4]);
+void procsSharingVertex(const int vertex, std::vector<int>& procIds);
 
 bool belongToTria(double const* x, double const* t, double bcoords[3], double eps = 1e-3);