Implemented alignment to PIV field. It compiles but it needs to be tested.

Author: sknepnek

CMakeBoostSetup.txt

@@ -44,7 +44,7 @@ set(Boost_USE_MULTITHREAD "OFF")
 set(Boost_ADDITIONAL_VERSIONS "1.55.0" "1.55" "1.54" "1.54.0" "1.53" "1.53.0" "1.52" "1.52.0" "1.51" "1.51.0" "1.50" "1.50.0" "1.49" "1.49.0" "1.48" "1.48.0" "1.48.0.2" "1.47" "1.47.0" "1.46.1" "1.46" "1.46.0" "1.45" "1.45.0" "1.44" "1.44.0" "1.42" "1.42.0" "1.41.0" "1.41" "1.40.0" "1.40" "1.39.0" "1.39" "1.38.0")
 
 # first, see if we can get any supported version of Boost
-find_package(Boost COMPONENTS regex iostreams REQUIRED)
+find_package(Boost COMPONENTS regex iostreams filesystem REQUIRED)
 
 
 # if we get boost 1.35 or greater, we need to get the system library too

CMakeLists.txt

@@ -63,7 +63,7 @@ include (CMakeCGALSetup.txt)
 
 ################################
 ## Define common libraries used by every target in MEMBRANE
-set(BOOST_LIBS 	${Boost_REGEX_LIBRARY} ${Boost_IOSTREAMS_LIBRARY} )
+set(BOOST_LIBS 	${Boost_REGEX_LIBRARY} ${Boost_IOSTREAMS_LIBRARY} ${Boost_FILESYSTEM_LIBRARY} )
 set(GSL_LIBS ${GSL_LIBRARIES})
 set(MATH_LIB -lm)
 set(SAMoS_LIBS ${GSL_LIBS} ${MATH_LIB} ${BOOST_LIBS} ${VTK_LIBS})

src/aligner/external/external_piv_aligner.cpp

@@ -33,39 +33,84 @@
 void ExternalPIVAlign::compute()
 {
   int N = m_system->size();
-  double hx = m_hx, hy = m_hy, hz = m_hz;  
-  
-  for  (int i = 0; i < N; i++)
+  double gamma = m_gamma;
+  double dx = (m_xhi - m_xlow)/(m_nx-1), dy = (m_yhi - m_ylow)/(m_ny - 1);
+
+  double t = static_cast<double>(m_step_counter) / static_cast<double>(m_n_piv_steps);
+  for (int i = 0; i < N; i++)
   {
     Particle& pi = m_system->get_particle(i);
     if (m_has_params)
+      gamma = m_type_params[pi.get_type()-1].gamma;
+    if (pi.x > m_xhi || pi.x < m_xlow || pi.y > m_yhi || pi.y < m_ylow)
+      throw runtime_error("One of the cells is outside the PIV mesh range. Please confien the system or make a larger PIV mesh");
+
+    int I = int((pi.x - m_xlow) / dx);
+    int J = int((pi.y - m_ylow) / dy);
+
+    double u_c = m_u[m_piv_frame_counter].piv_data[I][J], u_n = m_u[m_piv_frame_counter + 1].piv_data[I][J];
+    double v_c = m_v[m_piv_frame_counter].piv_data[I][J], v_n = m_v[m_piv_frame_counter + 1].piv_data[I][J];
+
+    double u = t * u_c + (1.0 - t) * u_n;
+    double v = t * v_c + (1.0 - t) * v_n;
+
+    double vel = sqrt(u * u + v * v);
+    if (vel > 0.0)
     {
-      hx = m_type_params[pi.get_type()-1].hx;
-      hy = m_type_params[pi.get_type()-1].hy;
-      hz = m_type_params[pi.get_type()-1].hz;
+      u *= gamma/vel;
+      v *= gamma/vel;
     }
-    pi.tau_x += hy*pi.nz - hz*pi.ny;
-    pi.tau_y += hz*pi.nx - hx*pi.nz;
-    pi.tau_z += hx*pi.ny - hy*pi.nx;
+    pi.tau_z += u*pi.ny - v*pi.nx;
+  
   }
+  m_step_counter++;
+  if (m_step_counter % m_n_piv_steps == 0)
+  {
+    m_piv_frame_counter++;
+    m_step_counter = 0;
+  }
+  if (m_piv_frame_counter >= m_u.size())
+    throw runtime_error("Simulation time exceeded avaliable PIV snapshots.");
 }
 
 // Private member function
-void ExternalPIVAlign::__read_csv(const string & fame, PIVData & data)
+void ExternalPIVAlign::__read_csv(const string & fname, PIVData & data)
 {
-  ifstream f(fname.c_str());
-  while (f.good())
+  std::ifstream f(fname.c_str());
+  if (f)
   {
-    string line;
-    getline(f, line);
-    if (line.length() > 0)
+    while (f.good())
     {
-      istringstream buffer(line);
-      string sval;
-      vector<double> fline;
-      while (getline(buffer, sval, ','))
-        fline.push_back(stod(sval));
-      data.piv_data.push_back(fline);
+      string line;
+      std::getline(f, line);
+      if (line.length() > 0)
+      {
+        std::istringstream buffer(line);
+        string sval;
+        vector<double> fline;
+        while (getline(buffer, sval, ','))
+          fline.push_back(stod(sval));
+        data.piv_data.push_back(fline);
+      }
     }
   }
+  else
+  {
+    throw runtime_error("Could not open file " + fname + ".");
+  }
+}
+
+// Get a list of all files matching base name
+void ExternalPIVAlign::__get_file_list(const string & base_name, vector<string> & all_matching_files)
+{
+  const boost::regex filter(base_name + ".*\\."+m_ext);
+  directory_iterator end_itr;
+  for(directory_iterator i(m_path); i != end_itr; ++i)
+  { 
+     if( !is_regular_file(i->status()) ) continue;
+     boost::smatch what;
+     if( !boost::regex_match(i->path().filename().string(), what, filter) ) continue;
+     all_matching_files.push_back(m_path+i->path().filename().string());
+  } 
+  sort(all_matching_files.begin(), all_matching_files.end());
 }