diff --git a/src/writer.cpp b/src/writer.cpp index 946db78..da53ae4 100644 --- a/src/writer.cpp +++ b/src/writer.cpp @@ -80,6 +80,94 @@ void dump_particles(std::string const & name, std::vector const ofs << "\n"; } +void dump_particles(std::string const & dir, size_t count, std::vector const & x, + std::vector const & theta, + std::vector const & v, + std::vector const & a, + std::vector const & omega, + std::vector const & alpha, + double r_part) { + std::stringstream out_file_name; + out_file_name << dir << "/particles_" << count << ".vtk"; + std::ofstream ofs(out_file_name.str()); + + if (!ofs.good()) { + std::cerr << "Unable to create a dump file at " << out_file_name.str() << std::endl; + exit(EXIT_FAILURE); + } + + ofs << "# vtk DataFile Version 4.0" << "\n"; + ofs << "Generated by libFractalCommon" << "\n"; + ofs << "ASCII" << "\n"; + ofs << "DATASET POLYDATA" << "\n"; + ofs << "POINTS " << x.size() << " FLOAT" << "\n"; + + // Coordinates are normalized to avoid issues with ParaView and small particles + for (auto const & p : x) { + ofs << p[0] / r_part << " " << p[1] / r_part << " " << p[2] / r_part << " "; + } + + ofs << "\n\n"; + ofs << "POINT_DATA " << x.size() << "\n"; + ofs << "FIELD FieldData 7" << "\n"; + ofs << "v 1 " << x.size() << " double\n"; + for (auto const & p : v) { + ofs << p.norm() << " "; + } + ofs << "\n"; + ofs << "a 1 " << x.size() << " double\n"; + for (auto const & p : a) { + ofs << p.norm() << " "; + } + ofs << "\n"; + ofs << "omega 1 " << x.size() << " double\n"; + for (auto const & p : omega) { + ofs << p.norm() << " "; + } + ofs << "\n"; + ofs << "alpha 1 " << x.size() << " double\n"; + for (auto const & p : alpha) { + ofs << p.norm() << " "; + } + ofs << "\n"; + ofs << "thetax 3 " << x.size() << " double\n"; + for (auto const & t : theta) { + Eigen::Matrix3d m; + m = Eigen::AngleAxis(t[0], Eigen::Vector3d::UnitX()) + * Eigen::AngleAxis(t[1], Eigen::Vector3d::UnitY()) + * Eigen::AngleAxis(t[2], Eigen::Vector3d::UnitZ()); + + auto unit = Eigen::Vector3d::UnitX(); + auto orient = m*unit; + + ofs << orient[0] << " " << orient[1] << " " << orient[2] << " "; + } + ofs << "\n" << "thetay 3 " << x.size() << " double" << "\n"; + for (auto const & t : theta) { + Eigen::Matrix3d m; + m = Eigen::AngleAxis(t[0], Eigen::Vector3d::UnitX()) + * Eigen::AngleAxis(t[1], Eigen::Vector3d::UnitY()) + * Eigen::AngleAxis(t[2], Eigen::Vector3d::UnitZ()); + + auto unit = Eigen::Vector3d::UnitY(); + auto orient = m*unit; + + ofs << orient[0] << " " << orient[1] << " " << orient[2] << " "; + } + ofs << "\n" << "thetaz 3 " << x.size() << " double" << "\n"; + for (auto const & t : theta) { + Eigen::Matrix3d m; + m = Eigen::AngleAxis(t[0], Eigen::Vector3d::UnitX()) + * Eigen::AngleAxis(t[1], Eigen::Vector3d::UnitY()) + * Eigen::AngleAxis(t[2], Eigen::Vector3d::UnitZ()); + + auto unit = Eigen::Vector3d::UnitZ(); + auto orient = m*unit; + + ofs << orient[0] << " " << orient[1] << " " << orient[2] << " "; + } +} + void dump_particles(std::string const & dir, size_t count, std::vector const & x, std::vector const & v, std::vector const & a, diff --git a/src/writer.h b/src/writer.h index 3755d93..5f94202 100644 --- a/src/writer.h +++ b/src/writer.h @@ -70,4 +70,12 @@ void dump_particles(std::string const & dir, size_t count, std::vector const & alpha, double r_part); +void dump_particles(std::string const & dir, size_t count, std::vector const & x, + std::vector const & theta, + std::vector const & v, + std::vector const & a, + std::vector const & omega, + std::vector const & alpha, + double r_part); + #endif //SOOT_AFM_WRITER_H diff --git a/validation/01_particles_colliding.cpp b/validation/01_particles_colliding.cpp index e7b5cec..d856fa9 100644 --- a/validation/01_particles_colliding.cpp +++ b/validation/01_particles_colliding.cpp @@ -12,6 +12,7 @@ #include #include "../src/energy.h" +#include "../src/writer.h" static const std::filesystem::path OUTPUT_FILE = "/dev/stdout"; @@ -28,6 +29,8 @@ int main() { const double dt = 1e-13; const double t_tot = 1.0e-10; const auto n_steps = size_t(t_tot / dt); + const size_t n_dumps = 300; + const size_t dump_period = n_steps / n_dumps; // General parameters const double rho = 1700.0; @@ -83,6 +86,12 @@ int main() { for (size_t n = 0; n < n_steps; n ++) { +// if (n % dump_period == 0) { +// std::cout << "Dump #" << n / dump_period << std::endl; +// dump_particles("run", n / dump_period, system.get_x(), system.get_theta(), +// system.get_v(), system.get_a(), +// system.get_omega(), system.get_alpha(), r_part); +// } if (double(n) * dt * 1e6 <= 4.0e-6) { t_span.emplace_back(double(n) * dt * 1e6); ke_trs_span.emplace_back(compute_ke_trs(system.get_v(), mass)); diff --git a/validation/02_particles_colliding.cpp b/validation/02_particles_colliding.cpp index e982217..5b95d6e 100644 --- a/validation/02_particles_colliding.cpp +++ b/validation/02_particles_colliding.cpp @@ -12,6 +12,7 @@ #include #include "../src/energy.h" +#include "../src/writer.h" static const std::filesystem::path OUTPUT_FILE = "/dev/stdout"; @@ -28,7 +29,7 @@ int main() { const double dt = 1e-13; const double t_tot = 1.0e-10; const auto n_steps = size_t(t_tot / dt); - const size_t n_dumps = 1000; + const size_t n_dumps = 300; const size_t dump_period = n_steps / n_dumps; // General parameters @@ -85,6 +86,12 @@ int main() { for (size_t n = 0; n < n_steps; n ++) { +// if (n % dump_period == 0) { +// std::cout << "Dump #" << n / dump_period << std::endl; +// dump_particles("run", n / dump_period, system.get_x(), system.get_theta(), +// system.get_v(), system.get_a(), +// system.get_omega(), system.get_alpha(), r_part); +// } if (double(n) * dt * 1e6 <= 4.0e-6) { t_span.emplace_back(double(n) * dt * 1e6); ke_trs_span.emplace_back(compute_ke_trs(system.get_v(), mass)); diff --git a/validation/03_particles_colliding.cpp b/validation/03_particles_colliding.cpp index 7dd634b..a377548 100644 --- a/validation/03_particles_colliding.cpp +++ b/validation/03_particles_colliding.cpp @@ -12,6 +12,7 @@ #include #include "../src/energy.h" +#include "../src/writer.h" static const std::filesystem::path OUTPUT_FILE = "/dev/stdout"; @@ -28,6 +29,8 @@ int main() { const double dt = 1e-13; const double t_tot = 1.0e-10; const auto n_steps = size_t(t_tot / dt); + const size_t n_dumps = 300; + const size_t dump_period = n_steps / n_dumps; // General parameters const double rho = 1700.0; @@ -83,6 +86,12 @@ int main() { for (size_t n = 0; n < n_steps; n ++) { +// if (n % dump_period == 0) { +// std::cout << "Dump #" << n / dump_period << std::endl; +// dump_particles("run", n / dump_period, system.get_x(), system.get_theta(), +// system.get_v(), system.get_a(), +// system.get_omega(), system.get_alpha(), r_part); +// } if (double(n) * dt * 1e6 <= 4.0e-6) { t_span.emplace_back(double(n) * dt * 1e6); ke_trs_span.emplace_back(compute_ke_trs(system.get_v(), mass)); diff --git a/validation/04_particles_colliding.cpp b/validation/04_particles_colliding.cpp index db9f03b..91005bd 100644 --- a/validation/04_particles_colliding.cpp +++ b/validation/04_particles_colliding.cpp @@ -12,6 +12,7 @@ #include #include "../src/energy.h" +#include "../src/writer.h" static const std::filesystem::path OUTPUT_FILE = "/dev/stdout"; @@ -28,6 +29,8 @@ int main() { const double dt = 1e-13; const double t_tot = 1.0e-10; const auto n_steps = size_t(t_tot / dt); + const size_t n_dumps = 300; + const size_t dump_period = n_steps / n_dumps; // General parameters const double rho = 1700.0; @@ -84,6 +87,12 @@ int main() { for (size_t n = 0; n < n_steps; n ++) { +// if (n % dump_period == 0) { +// std::cout << "Dump #" << n / dump_period << std::endl; +// dump_particles("run", n / dump_period, system.get_x(), system.get_theta(), +// system.get_v(), system.get_a(), +// system.get_omega(), system.get_alpha(), r_part); +// } if (double(n) * dt * 1e6 <= 4.0e-6) { t_span.emplace_back(double(n) * dt * 1e6); ke_trs_span.emplace_back(compute_ke_trs(system.get_v(), mass)); diff --git a/validation/05_hamaker.cpp b/validation/05_hamaker.cpp index 8f4676f..e209aa0 100644 --- a/validation/05_hamaker.cpp +++ b/validation/05_hamaker.cpp @@ -13,6 +13,7 @@ #include #include "../src/energy.h" +#include "../src/writer.h" #include "fixed_step_handler.h" static const std::filesystem::path OUTPUT_FILE = "/dev/stdout"; @@ -116,6 +117,12 @@ int main() { ofs << "D\tKE\tU\n"; for (size_t n = 0; n < n_steps; n ++) { +// if (n % dump_period == 0) { +// std::cout << "Dump #" << n / dump_period << std::endl; +// dump_particles("run", n / dump_period, system.get_x(), system.get_theta(), +// system.get_v(), system.get_a(), +// system.get_omega(), system.get_alpha(), r_part); +// } double separation = (system.get_x()[1] - system.get_x()[0]).norm() - 2.0 * r_part; if (separation <= h0) { has_collided = true; diff --git a/validation/06_sintering.cpp b/validation/06_sintering.cpp index ffd91a9..7b3e9c3 100644 --- a/validation/06_sintering.cpp +++ b/validation/06_sintering.cpp @@ -13,6 +13,7 @@ #include #include "../src/energy.h" +#include "../src/writer.h" static const std::filesystem::path OUTPUT_FILE = "/dev/stdout"; @@ -104,7 +105,11 @@ int main() { for (size_t n = 0; n < n_steps; n ++) { // if (n % dump_period == 0) { // std::cout << "Dump #" << n / dump_period << std::endl; -// write_particles("run", system.get_x(), system.get_theta(), r_part); +// dump_particles("run", n / dump_period, system.get_x(), system.get_theta(), +// system.get_v(), system.get_a(), +// system.get_omega(), system.get_alpha(), r_part); +// dump_necks("run", n / dump_period, system.get_x(), std::vector{false, true, false, true, false, true, false, false, +// true}, r_part); // } if (n % thermo_dump_period == 0) { auto ke_trs = compute_ke_trs(system.get_v(), mass); diff --git a/validation/07_sintering.cpp b/validation/07_sintering.cpp index 6699ce8..2dbd12d 100644 --- a/validation/07_sintering.cpp +++ b/validation/07_sintering.cpp @@ -13,6 +13,7 @@ #include #include "../src/energy.h" +#include "../src/writer.h" static const std::filesystem::path OUTPUT_FILE = "/dev/stdout"; @@ -104,7 +105,11 @@ int main() { for (size_t n = 0; n < n_steps; n ++) { // if (n % dump_period == 0) { // std::cout << "Dump #" << n / dump_period << std::endl; -// write_particles("run", system.get_x(), system.get_theta(), r_part); +// dump_particles("run", n / dump_period, system.get_x(), system.get_theta(), +// system.get_v(), system.get_a(), +// system.get_omega(), system.get_alpha(), r_part); +// dump_necks("run", n / dump_period, system.get_x(), std::vector{false, true, false, true, false, true, false, false, +// true}, r_part); // } if (n % thermo_dump_period == 0) { auto ke_trs = compute_ke_trs(system.get_v(), mass); diff --git a/validation/08_sintering.cpp b/validation/08_sintering.cpp index 55451f5..8e90fdd 100644 --- a/validation/08_sintering.cpp +++ b/validation/08_sintering.cpp @@ -13,6 +13,7 @@ #include #include "../src/energy.h" +#include "../src/writer.h" static const std::filesystem::path OUTPUT_FILE = "/dev/stdout"; @@ -104,7 +105,11 @@ int main() { for (size_t n = 0; n < n_steps; n ++) { // if (n % dump_period == 0) { // std::cout << "Dump #" << n / dump_period << std::endl; -// write_particles("run", system.get_x(), system.get_theta(), r_part); +// dump_particles("run", n / dump_period, system.get_x(), system.get_theta(), +// system.get_v(), system.get_a(), +// system.get_omega(), system.get_alpha(), r_part); +// dump_necks("run", n / dump_period, system.get_x(), std::vector{false, true, false, true, false, true, false, false, +// true}, r_part); // } if (n % thermo_dump_period == 0) { auto ke_trs = compute_ke_trs(system.get_v(), mass); diff --git a/validation/09_sintering.cpp b/validation/09_sintering.cpp index ddc93c3..f508945 100644 --- a/validation/09_sintering.cpp +++ b/validation/09_sintering.cpp @@ -13,6 +13,7 @@ #include #include "../src/energy.h" +#include "../src/writer.h" static const std::filesystem::path OUTPUT_FILE = "/dev/stdout"; @@ -104,7 +105,11 @@ int main() { for (size_t n = 0; n < n_steps; n ++) { // if (n % dump_period == 0) { // std::cout << "Dump #" << n / dump_period << std::endl; -// write_particles("run", system.get_x(), system.get_theta(), r_part); +// dump_particles("run", n / dump_period, system.get_x(), system.get_theta(), +// system.get_v(), system.get_a(), +// system.get_omega(), system.get_alpha(), r_part); +// dump_necks("run", n / dump_period, system.get_x(), std::vector{false, true, false, true, false, true, false, false, +// true}, r_part); // } if (n % thermo_dump_period == 0) { auto ke_trs = compute_ke_trs(system.get_v(), mass);