From 20ff8f986932d53d29b617aa94124d8da7673e1f Mon Sep 17 00:00:00 2001
From: Joachim Schoeberl <joachim.schoeberl@tuwien.ac.at>
Date: Mon, 3 Sep 2012 09:49:43 +0000
Subject: [PATCH] parallel curved elements

---
 libsrc/meshing/curvedelems.cpp  | 82 +++++++++++++++++++++++++++++++--
 libsrc/meshing/parallelmesh.cpp |  2 +-
 libsrc/meshing/paralleltop.cpp  |  2 +-
 3 files changed, 79 insertions(+), 7 deletions(-)

diff --git a/libsrc/meshing/curvedelems.cpp b/libsrc/meshing/curvedelems.cpp
index b938c976..a7c2acaf 100644
--- a/libsrc/meshing/curvedelems.cpp
+++ b/libsrc/meshing/curvedelems.cpp
@@ -11,6 +11,74 @@ namespace netgen
   //   bool rational = true;
 
 
+  template <typename T>
+  inline void MyMPI_ExchangeTableTmp (TABLE<T> & send_data, 
+				   TABLE<T> & recv_data, int tag,
+				   MPI_Comm comm = MPI_COMM_WORLD)
+  {
+    int ntasks, rank;
+    MPI_Comm_size(comm, &ntasks);
+    MPI_Comm_rank(comm, &rank);
+
+    Array<int> send_sizes(ntasks);
+    Array<int> recv_sizes(ntasks);
+    for (int i = 0; i < ntasks; i++)
+      send_sizes[i] = send_data[i].Size();
+    
+    cout << "id = " << rank << ", sendsize = " << send_sizes << endl;
+    MPI_Alltoall (&send_sizes[0], 1, MPI_INT, 
+		  &recv_sizes[0], 1, MPI_INT, comm);
+    /*
+      // in-place is buggy !
+    MPI_Alltoall (MPI_IN_PLACE, 1, MPI_INT, 
+		  &recv_sizes[0], 1, MPI_INT, comm);
+    */
+    MPI_Barrier(comm);
+    cout << "id = " << rank << ", recvsize = " << recv_sizes << endl;
+    MPI_Barrier(comm);
+
+    cout << "recv_table = " << recv_data << endl;
+
+    for (int i = 0; i < ntasks; i++)
+      if (recv_sizes[i])
+	recv_data.SetEntrySize (i, recv_sizes[i], sizeof(T));
+
+    cout << "recv_table2 = " << recv_data << endl;
+
+
+    MPI_Barrier(comm);
+    for (int i = 0; i < ntasks; i++)
+      cout << "id = " << rank << ", recvtagsize(" << i << ") = " << recv_data[i].Size() << endl;
+    MPI_Barrier(comm);
+
+    for (int i = 0; i < ntasks; i++)
+      if (recv_data[i].Size())
+	recv_data[i] = 0;
+    
+    Array<MPI_Request> requests(0);
+    for (int dest = 0; dest < ntasks; dest++)
+      if ( (dest != rank) && send_data[dest].Size())
+	requests.Append (MyMPI_ISend (send_data[dest], dest, tag, comm));
+
+    for (int dest = 0; dest < ntasks; dest++)
+      if ( (dest != rank) && recv_data[dest].Size())
+	requests.Append (MyMPI_IRecv (recv_data[dest], dest, tag, comm));
+    
+    MPI_Barrier(comm);
+    cout << "comm started, id = " << rank << ", num requ = " << requests.Size() << endl;
+    MPI_Barrier(comm);
+    
+    if (requests.Size())
+      MPI_Waitall (requests.Size(), &requests[0], MPI_STATUS_IGNORE);
+    MPI_Barrier(comm);
+    cout << "table done, id = " << rank << endl;
+    sleep(1);
+    MPI_Barrier(comm);
+  }
+
+
+
+
   
   static void ComputeGaussRule (int n, Array<double> & xi, Array<double> & wi)
   {
@@ -599,6 +667,7 @@ namespace netgen
 	  }
       }
 
+    
     MyMPI_ExchangeTable (send_orders, recv_orders, MPI_TAG_CURVE, curve_comm);
 
     if (ntasks > 1 && working)
@@ -667,8 +736,6 @@ namespace netgen
 	  jacpols2[i] = new JacobiRecPol (100, i, 2);
       }
 
-
-
     PrintMessage (3, "Curving edges");
 
     if (mesh.GetDimension() == 3 || rational)
@@ -689,7 +756,7 @@ namespace netgen
 		{
 		  PointIndex pi1 = el[edges[i2][0]];
 		  PointIndex pi2 = el[edges[i2][1]];
-		
+
 		  bool swap = pi1 > pi2;
 		
 		  Point<3> p1 = mesh[pi1];
@@ -704,6 +771,7 @@ namespace netgen
 		}
 	    }
 
+
 #ifdef PARALLEL
 	if (ntasks > 1)
 	  {
@@ -727,8 +795,10 @@ namespace netgen
 			}
 		    }
 		}
+	    
 	    MyMPI_ExchangeTable (senddata, recvdata, MPI_TAG_CURVE, curve_comm);
 	    
+
 	    Array<int> cnt(ntasks);
 	    cnt = 0;
 	    if (working)
@@ -754,6 +824,7 @@ namespace netgen
 	  }
 #endif    
 
+
 	if (working)
 	  for (int e = 0; e < surfnr.Size(); e++)
 	    {
@@ -943,9 +1014,10 @@ namespace netgen
 	      partop.GetDistantEdgeNums (e+1, procs);
 	      for (int j = 0; j < procs.Size(); j++)
 		{
-		  use_edge[e] = int(recvdata[procs[j]][cnt[procs[j]]++]);
-		  if (use_edge[e])
+		  int get_edge = int(recvdata[procs[j]][cnt[procs[j]]++]);
+		  if (get_edge)
 		    {
+		      use_edge[e] = 1;
 		      edge_surfnr1[e] = int (recvdata[procs[j]][cnt[procs[j]]++]);
 		      edge_surfnr2[e] = int (recvdata[procs[j]][cnt[procs[j]]++]);
 		      edge_gi0[e].edgenr = int (recvdata[procs[j]][cnt[procs[j]]++]);
diff --git a/libsrc/meshing/parallelmesh.cpp b/libsrc/meshing/parallelmesh.cpp
index 55543550..43644f8b 100644
--- a/libsrc/meshing/parallelmesh.cpp
+++ b/libsrc/meshing/parallelmesh.cpp
@@ -706,7 +706,7 @@ namespace netgen
 
     NgProfiler::RegionTimer regloc(timerloc);
     stringstream str;
-    str << "Got " << GetNE() << " elements and " 
+    str << "p" << id << ": got " << GetNE() << " elements and " 
 	 << GetNSE() << " surface elements";
     cout << str.str() << endl;
     // PrintMessage (2, "Got ", GetNE(), " elements and ", GetNSE(), " surface elements");
diff --git a/libsrc/meshing/paralleltop.cpp b/libsrc/meshing/paralleltop.cpp
index 3048fa2c..3f5ac9f2 100644
--- a/libsrc/meshing/paralleltop.cpp
+++ b/libsrc/meshing/paralleltop.cpp
@@ -211,7 +211,7 @@ namespace netgen
 
     (*testout) << "UPDATE COARSE GRID PARALLEL TOPOLOGY " << endl;
     if (id == 0)
-      PrintMessage (1, "UPDATE COARSE GRID PARALLEL TOPOLOGY ");
+      PrintMessage (1, "update parallel topology");
 
 
     // UpdateCoarseGridGlobal();