From 8850969187b18ca759201631190ca9f001d4edc4 Mon Sep 17 00:00:00 2001
From: Guilhem Lavaux <guilhem.lavaux@gmail.com>
Date: Sun, 15 Jan 2012 16:39:32 -0500
Subject: [PATCH] Hack the voz1b1 mesh

---
 mytools/hack_mesh.cpp | 88 +++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 88 insertions(+)
 create mode 100644 mytools/hack_mesh.cpp

diff --git a/mytools/hack_mesh.cpp b/mytools/hack_mesh.cpp
new file mode 100644
index 0000000..6cec2e1
--- /dev/null
+++ b/mytools/hack_mesh.cpp
@@ -0,0 +1,88 @@
+#include <iostream>
+#include <fstream>
+#include <netcdfcpp.h>
+#include <vector>
+#include <set>
+
+using namespace std;
+
+int main(int argc, char **argv)
+{
+  char *n =  argv[1];
+  char *n2 = argv[2];
+  ifstream f(n);
+  ofstream f2(n2);
+  int np, nvp;
+  int *ids;
+  float *vols;
+  set<int> mask_id;
+
+  NcFile fp("params.nc");
+  NcVar *v_id = fp.get_var("particle_ids");
+  long int *e_N = v_id->edges();
+  int *pid = new int[e_N[0]];
+
+  v_id->get(pid, e_N[0]);
+
+  for (int i = 0; i < e_N[0]; i++)
+    mask_id.insert(pid[i]);
+  delete[] pid;
+  delete[] e_N;
+
+  f.read((char*)&np, sizeof(int));
+  f.read((char*)&nvp,sizeof(int));
+  f2.write((char*)&np, sizeof(int));
+
+  ids = new int[nvp];
+  vols = new float[nvp];
+  f.read((char*)ids, sizeof(int)*nvp);
+  f.read((char *)vols, sizeof(float)*nvp);
+
+  vector< vector<int> > adjs;
+  vector<int> adj0, new_ids;
+  vector<float>  new_vols;
+  for (int i = 0; i < nvp; i++)
+    {
+      int nadj;
+      bool  masked_particle = false;
+
+
+      if (mask_id.find(ids[i]) != mask_id.end())
+	{
+	  masked_particle = true;
+	}
+
+      f.read((char*)&nadj, sizeof(int));
+      adj0.resize(nadj);
+      
+      if (nadj != 0)
+	{
+	  f.read((char*)&adj0[0], sizeof(float)*nadj);
+	  for (int j = 0; j < nadj; j++)
+	    if (mask_id.find(adj0[j]) != mask_id.end())
+	      masked_particle = true;
+	}
+      
+      if (!masked_particle)
+	{
+	  new_vols.push_back(vols[i]);
+	  new_ids.push_back(ids[i]);
+	  adjs.push_back(adj0);
+	}
+    }
+
+  nvp = new_ids.size();
+  f2.write((char*)&nvp, sizeof(int));
+  f2.write((char*)&new_ids[0], sizeof(int)*nvp);
+  f2.write((char*)&new_vols[0], sizeof(float)*nvp);
+  
+  for (int i = 0; i < nvp; i++)
+    {
+      int nadj = adjs[i].size();
+      f2.write((char*)&nadj, sizeof(int));
+      if (nadj != 0)
+	f2.write((char*)&adjs[i], sizeof(int)*nadj);
+    }
+
+  return 0;
+}