diff --git a/CMakeLists.txt b/CMakeLists.txt index 56b29195..8c54679a 100644 --- a/CMakeLists.txt +++ b/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) diff --git a/TestingFramework.cmake b/TestingFramework.cmake index e4130c33..090a779b 100644 --- a/TestingFramework.cmake +++ b/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) diff --git a/src/main.cpp b/src/main.cpp index 015e8665..4b36af79 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -4,19 +4,22 @@ #include #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 xvf, yvf; + // Contact force at (xvf,yvf) calculated in the last iteration. + std::vector 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(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); } diff --git a/src/mirco_evaluate.cpp b/src/mirco_evaluate.cpp index 6bfa7100..b0210891 100644 --- a/src/mirco_evaluate.cpp +++ b/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& topology, - double zmax, std::vector& meshgrid) + double zmax, std::vector& meshgrid, std::vector& xvf, std::vector& yvf, + std::vector& 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 xv0, yv0; - // Coordinates of the points in contact in the previous iteration. - std::vector xvf, yvf; // Indentation value of the half space at the predicted points of contact. std::vector b0; - // Contact force at (xvf,yvf) predicted in the previous iteration. - std::vector 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 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 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; } diff --git a/src/mirco_evaluate.h b/src/mirco_evaluate.h index ce120779..089e88f4 100644 --- a/src/mirco_evaluate.h +++ b/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& topology, - double zmax, std::vector& meshgrid); + double zmax, std::vector& meshgrid, std::vector& xvf, + std::vector& yvf, std::vector& pf, int& n); } // namespace MIRCO diff --git a/src/mirco_postprocess.cpp b/src/mirco_postprocess.cpp new file mode 100644 index 00000000..6f2356cd --- /dev/null +++ b/src/mirco_postprocess.cpp @@ -0,0 +1,82 @@ +#include "mirco_postprocess.h" + +#include +#include +#include +#include +#include +#include + +long findClosestIndex(const std::vector& 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 xvf, std::vector yvf, std::vector pf, + int nf, double GridSize, int ngrid, std::vector& 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(); +} \ No newline at end of file diff --git a/src/mirco_postprocess.h b/src/mirco_postprocess.h new file mode 100644 index 00000000..3183ea1a --- /dev/null +++ b/src/mirco_postprocess.h @@ -0,0 +1,29 @@ +#ifndef SRC_POSTPROCESS_H_ +#define SRC_POSTPROCESS_H_ + +#include +#include + +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 xvf, std::vector yvf, std::vector pf, int nf, + double GridSize, int ngrid, std::vector& meshgrid, double LateralLength, + double elapsedTime, std::string outputFileName); +} // namespace MIRCO + + +#endif // SRC_POSTPROCESS_H_