Implement postprocess

CMakeLists.txt

@@ -98,6 +98,7 @@ target_link_libraries(mirco_core PUBLIC Trilinos::all_selected_libs OpenMP::Open
 
 add_library(mirco_io
   src/mirco_setparameters.cpp
+  src/mirco_postprocess.cpp
   )
 target_link_libraries(mirco_io PUBLIC Trilinos::all_selected_libs)
 

TestingFramework.cmake

@@ -2,7 +2,7 @@
 macro(mirco_framework_test name_of_input_file)
   set (input_location ${PROJECT_SOURCE_DIR}/Input/${name_of_input_file})
   add_test(NAME ${name_of_input_file}
-  COMMAND ./mirco ${input_location})
+  COMMAND ./mirco ${input_location} nametest)
 
 endmacro(mirco_framework_test)
 

src/main.cpp

@@ -4,19 +4,22 @@
 #include <string>
 
 #include "mirco_evaluate.h"
+#include "mirco_filesystem_utils.h"
+#include "mirco_postprocess.h"
 #include "mirco_setparameters.h"
 #include "mirco_topology.h"
 #include "mirco_topologyutilities.h"
 
-
 int main(int argc, char* argv[])
 {
-  TEUCHOS_TEST_FOR_EXCEPTION(
-      argc != 2, std::invalid_argument, "The code expects (only) an input file as argument");
+  TEUCHOS_TEST_FOR_EXCEPTION(argc != 3, std::invalid_argument,
+      "The code expects input file with a full path and a prefix for output files.");
   // reading the input file name from the command line
   std::string inputFileName = argv[1];
+  std::string outputFileName = argv[2];
 
-  const auto start = std::chrono::high_resolution_clock::now();
+  // following function generates the actual path of the output file.
+  UTILS::ChangeRelativePath(outputFileName, inputFileName);
 
   bool WarmStartingFlag = false;
   int Resolution = 0;
@@ -55,17 +58,25 @@ int main(int argc, char* argv[])
 
   auto max_and_mean = MIRCO::ComputeMaxAndMean(topology);
 
-  // Initialise Pressure
-  double pressure = 0.0;
+  // Number of contact points
+  int nf = 0;
+  // Coordinates of the points in contact in the last iteration.
+  std::vector<double> xvf, yvf;
+  // Contact force at (xvf,yvf) calculated in the last iteration.
+  std::vector<double> pf;
 
-  MIRCO::Evaluate(pressure, Delta, LateralLength, GridSize, Tolerance, MaxIteration,
-      CompositeYoungs, WarmStartingFlag, ElasticComplianceCorrection, topology, max_and_mean.max_,
-      meshgrid);
+  // Set start time
+  const auto start = std::chrono::high_resolution_clock::now();
 
-  std::cout << "Mean pressure is: " << std::to_string(pressure) << std::endl;
+  MIRCO::Evaluate(Delta, LateralLength, GridSize, Tolerance, MaxIteration, CompositeYoungs,
+      WarmStartingFlag, ElasticComplianceCorrection, topology, max_and_mean.max_, meshgrid, xvf,
+      yvf, pf, nf);
 
+  // Total MIRCO simulation time
   const auto finish = std::chrono::high_resolution_clock::now();
   const double elapsedTime =
       std::chrono::duration_cast<std::chrono::seconds>(finish - start).count();
-  std::cout << "Elapsed time is: " + std::to_string(elapsedTime) + "s." << std::endl;
+
+  MIRCO::PostProcess(
+      xvf, yvf, pf, nf, GridSize, ngrid, meshgrid, LateralLength, elapsedTime, outputFileName);
 }

src/mirco_evaluate.cpp

@@ -21,10 +21,11 @@
 #include "mirco_nonlinearsolver.h"
 
 
-void MIRCO::Evaluate(double& pressure, double Delta, double LateralLength, double GridSize,
-    double Tolerance, int MaxIteration, double CompositeYoungs, bool WarmStartingFlag,
+void MIRCO::Evaluate(double Delta, double LateralLength, double GridSize, double Tolerance,
+    int MaxIteration, double CompositeYoungs, bool WarmStartingFlag,
     double ElasticComplianceCorrection, Teuchos::SerialDenseMatrix<int, double>& topology,
-    double zmax, std::vector<double>& meshgrid)
+    double zmax, std::vector<double>& meshgrid, std::vector<double>& xvf, std::vector<double>& yvf,
+    std::vector<double>& pf, int& nf)
 {
   omp_set_num_threads(6);  // 6 seems to be optimal
 
@@ -40,21 +41,14 @@ void MIRCO::Evaluate(double& pressure, double Delta, double LateralLength, doubl
 
   // Coordinates of the points predicted to be in contact.
   std::vector<double> xv0, yv0;
-  // Coordinates of the points in contact in the previous iteration.
-  std::vector<double> xvf, yvf;
   // Indentation value of the half space at the predicted points of contact.
   std::vector<double> b0;
-  // Contact force at (xvf,yvf) predicted in the previous iteration.
-  std::vector<double> pf;
 
   // x0 --> contact forces at (xvf,yvf) predicted in the previous iteration but
   // are a part of currect predicted contact set. x0 is calculated in the
   // Warmstart function to be used in the NNLS to accelerate the simulation.
   Teuchos::SerialDenseMatrix<int, double> x0;
 
-  // The number of nodes in contact in the previous iteration.
-  int nf = 0;
-
   // The influence coefficient matrix (Discrete version of Green Function)
   Teuchos::SerialDenseMatrix<int, double> A;
   // Solution containing force
@@ -121,12 +115,4 @@ void MIRCO::Evaluate(double& pressure, double Delta, double LateralLength, doubl
 
   TEUCHOS_TEST_FOR_EXCEPTION(ErrorForce > Tolerance, std::out_of_range,
       "The solution did not converge in the maximum number of iternations defined");
-  // @{
-
-  // Calculate the final force value at the end of the iteration.
-  const double force = force0[k - 1];
-
-  // Mean pressure
-  double sigmaz = force / pow(LateralLength, 2);
-  pressure = sigmaz;
 }

src/mirco_evaluate.h

@@ -9,7 +9,6 @@ namespace MIRCO
   /**
    * @brief Relate the far-field displacement with pressure
    *
-   * @param pressure Pressure
    * @param Delta Far-field displacement (Gap)
    * @param LateralLength Lateral side of the surface [micrometers]
    * @param GridSize Grid size
@@ -21,11 +20,16 @@ namespace MIRCO
    * @param topology Topology matrix containing heights
    * @param zmax Maximum height
    * @param meshgrid Meshgrid vector
+   * @param xvf x-Coordinate vector of points that are in contact, calculated in each iteration
+   * @param yvf y-Coordinate vector of points that are in contact, calculated in each iteration
+   * @param pf Contact force vector at (xvf,yvf) calculated in each iteration.
+   * @param nf Number of contact points
    */
-  void Evaluate(double& pressure, double Delta, double LateralLength, double GridSize,
-      double Tolerance, int MaxIteration, double CompositeYoungs, bool WarmStartingFlag,
+  void Evaluate(double Delta, double LateralLength, double GridSize, double Tolerance,
+      int MaxIteration, double CompositeYoungs, bool WarmStartingFlag,
       double ElasticComplianceCorrection, Teuchos::SerialDenseMatrix<int, double>& topology,
-      double zmax, std::vector<double>& meshgrid);
+      double zmax, std::vector<double>& meshgrid, std::vector<double>& xvf,
+      std::vector<double>& yvf, std::vector<double>& pf, int& n);
 }  // namespace MIRCO
 
 

src/mirco_postprocess.cpp

@@ -0,0 +1,82 @@
+#include "mirco_postprocess.h"
+
+#include <Teuchos_SerialDenseMatrix.hpp>
+#include <cmath>
+#include <fstream>
+#include <iostream>
+#include <numeric>
+#include <vector>
+
+long findClosestIndex(const std::vector<double>& sorted_array, double x)
+{
+  auto iter_geq = std::lower_bound(sorted_array.begin(), sorted_array.end(), x);
+
+  if (iter_geq == sorted_array.begin())
+  {
+    return 0;
+  }
+
+  double a = *(iter_geq - 1);
+  double b = *(iter_geq);
+
+  if (fabs(x - a) < fabs(x - b))
+  {
+    return iter_geq - sorted_array.begin() - 1;
+  }
+  return iter_geq - sorted_array.begin();
+}
+
+void MIRCO::PostProcess(std::vector<double> xvf, std::vector<double> yvf, std::vector<double> pf,
+    int nf, double GridSize, int ngrid, std::vector<double>& meshgrid, double LateralLength,
+    double elapsedTime, std::string outputFileName)
+{
+  double forceArray[ngrid][ngrid] = {0};
+  int ix;
+  int iy;
+  double areaContact;
+
+  for (int i = 0; i < nf; i++)
+  {
+    ix = findClosestIndex(meshgrid, xvf[i]);
+    iy = findClosestIndex(meshgrid, yvf[i]);
+    forceArray[ix][ngrid-1-iy] = pf[i];
+    // std::cout << ix << "\t" << iy  << std::endl;
+    // std::cout << xvf[i] << "\t" << yvf[i] << "\t" << pf[i] << std::endl;
+  }
+
+  // Store force at every point (including contact and non-contact points)
+  std::ofstream outputForceFile(outputFileName + "_force" + ".dat");
+
+  for (int i = 0; i < ngrid; i++)
+  {
+    for (int j = 0; j < ngrid; j++)
+    {
+      outputForceFile << forceArray[j][ngrid - 1 - i] << ' ';
+    }
+    outputForceFile << "\n";
+  }
+  outputForceFile.close();
+
+  // Calculate average pressure via total force
+  auto totalForce = std::accumulate(pf.begin(), pf.end(), 0.0);
+  double pressure = totalForce / pow(LateralLength, 2);
+
+  // Total contact area
+  areaContact = nf * (pow(GridSize, 2) / pow(LateralLength, 2)) * 100;
+
+  // Store overall information
+  std::ofstream infoFile(outputFileName + "_info" + ".csv");
+  infoFile << "total_force"
+           << "\t"
+           << "total_area"
+           << "\t"
+           << "ave_pressure"
+           << "\t"
+           << "n_contact"
+           << "\t"
+           << "time"
+           << "\n";
+  infoFile << totalForce << "\t" << areaContact << "\t" << pressure << "\t" << nf << "\t"
+           << std::to_string(elapsedTime);
+  infoFile.close();
+}

src/mirco_postprocess.h

@@ -0,0 +1,29 @@
+#ifndef SRC_POSTPROCESS_H_
+#define SRC_POSTPROCESS_H_
+
+#include <Teuchos_SerialDenseMatrix.hpp>
+#include <string>
+
+namespace MIRCO
+{
+  /**
+   * @brief Postprocessing of contact force and contact area
+   *
+   * @param xvf x-Coordinate vector of points that are in contact, calculated in each iteration
+   * @param yvf y-Coordinate vector of points that are in contact, calculated in each iteration
+   * @param pf Contact force vector at (xvf,yvf) calculated in each iteration.
+   * @param nf Number of contact points
+   * @param GridSize Grid size
+   * @param ngrid Number of grids
+   * @param meshgrid Meshgrid vector
+   * @param LateralLength Lateral side of the surface [micrometers]
+   * @param elapsedTimeString Number of contact points
+   * @param outputFileName Output file to store results
+   */
+  void PostProcess(std::vector<double> xvf, std::vector<double> yvf, std::vector<double> pf, int nf,
+      double GridSize, int ngrid, std::vector<double>& meshgrid, double LateralLength,
+      double elapsedTime, std::string outputFileName);
+}  // namespace MIRCO
+
+
+#endif  // SRC_POSTPROCESS_H_