MPS-Basic
Loading...
Searching...
No Matches
particles_exporter.cpp
Go to the documentation of this file.
1#include "particles_exporter.hpp"
2
3#include <fstream>
4#include <iomanip>
5#include <iostream>
6#include <sstream>
7
8using std::cerr;
9using std::endl;
10namespace fs = std::filesystem;
11
12bool ParticlesExporter::isBigEndian() const {
13 uint32_t i = 1; // 0x00000001
14 uint8_t first_byte = *reinterpret_cast<uint8_t*>(&i);
15 return first_byte == 0;
16}
17
18std::string ParticlesExporter::dataArrayBegin(
19 const std::string& type, const std::string& name, int numberOfComponents, const std::string& format, size_t offset
20) const {
21 // --- assertion ---
22 // allowed types and format https://docs.vtk.org/en/latest/design_documents/VTKFileFormats.html
23 // -----------------
24 assert(
25 type == "Int8" || type == "UInt8" || type == "Int16" || type == "UInt16" || type == "UInt32" ||
26 type == "Int32" || type == "UInt64" || type == "Int64" || type == "Float32" || type == "Float64"
27 ); // allowed data types
28 assert(format == "ascii" || format == "binary" || format == "appended"); // allowed formats
29 assert(
30 format == "appended" && offset != SIZE_MAX || format != "appended" && offset == SIZE_MAX
31 ); // offset is required for appended format
32
33 std::string ret;
34 ret += "<DataArray type=\"" + type + "\" Name=\"" + name + "\" NumberOfComponents=\"" +
35 std::to_string(numberOfComponents) + "\" format=\"" + format + "\"";
36 if (format == "appended") {
37 ret += " offset=\"" + std::to_string(offset) + "\"";
38 }
39 ret += ">";
40 return ret;
41}
42
43std::string ParticlesExporter::dataArrayEnd() const {
44 return "</DataArray>";
45}
46
47void ParticlesExporter::setParticles(const Particles& particles) {
48 this->particles = particles;
49}
50
51void ParticlesExporter::toProf(const fs::path& path, const double& time) {
52 std::ofstream ofs(path);
53 if (ofs.fail()) {
54 cerr << "cannot write " << path << endl;
55 std::exit(-1);
56 }
57
58 ofs << time << endl;
59 ofs << particles.size() << endl;
60 for (const auto& p : particles) {
61 ofs << static_cast<int>(p.type) << " ";
62 ofs << p.position.x() << " " << p.position.y() << " " << p.position.z() << " ";
63 ofs << p.velocity.x() << " " << p.velocity.y() << " " << p.velocity.z();
64 ofs << endl;
65 }
66}
67void ParticlesExporter::toVtuAscii(const fs::path& path, const double& time, const double& n0ForNumberDensity) {
68 std::ofstream ofs(path);
69 if (ofs.fail()) {
70 cerr << "cannot write " << path << endl;
71 std::exit(-1);
72 }
73
74 // header
75 ofs << "<?xml version='1.0' encoding='UTF-8'?>" << endl;
76 ofs << "<VTKFile type=\"UnstructuredGrid\" version=\"1.0\" header_type=\"UInt64\" byte_order=\"";
77 if (isBigEndian()) {
78 ofs << "BigEndian";
79 } else {
80 ofs << "LittleEndian";
81 }
82 ofs << "\">" << endl;
83 ofs << "<UnstructuredGrid>" << endl
84 << "<Piece NumberOfPoints=\"" << particles.size() << "\" NumberOfCells=\"" << particles.size() << "\">" << endl;
85
89 ofs << "<Points>" << endl;
90 ofs << dataArrayBegin("Float64", "Position", 3, "ascii") << endl;
91 for (const auto& p : particles) {
92 ofs << p.position.x() << " " << p.position.y() << " " << p.position.z() << endl;
93 }
94 ofs << dataArrayEnd() << endl;
95 ofs << "</Points>" << endl;
96
97 // ---------------------
98 // ----- PointData -----
99 // ---------------------
100 ofs << "<PointData>" << endl;
101 // Particle Type
102 ofs << dataArrayBegin("Int32", "Particle Type", 1, "ascii") << endl;
103 for (const auto& p : particles) {
104 ofs << static_cast<int32_t>(p.type) << endl;
105 }
106 ofs << dataArrayEnd() << endl;
107 // Velocity
108 ofs << dataArrayBegin("Float64", "Velocity", 3, "ascii") << endl;
109 for (const auto& p : particles) {
110 ofs << p.velocity.x() << " " << p.velocity.y() << " " << p.velocity.z() << endl;
111 }
112 ofs << dataArrayEnd() << endl;
113 // Pressure
114 ofs << dataArrayBegin("Float64", "Pressure", 1, "ascii") << endl;
115 for (const auto& p : particles) {
116 ofs << p.pressure << endl;
117 }
118 ofs << dataArrayEnd() << endl;
119 // Number Density
120 ofs << dataArrayBegin("Float64", "Number Density", 1, "ascii") << endl;
121 for (const auto& p : particles) {
122 ofs << p.numberDensity << endl;
123 }
124 ofs << dataArrayEnd() << endl;
125 // Number Density Ratio
126 ofs << dataArrayBegin("Float64", "Number Density Ratio", 1, "ascii") << endl;
127 for (const auto& p : particles) {
128 ofs << p.numberDensity / n0ForNumberDensity << endl;
129 }
130 ofs << dataArrayEnd() << endl;
131 // Boundary Condition
132 ofs << dataArrayBegin("Int32", "Boundary Condition", 1, "ascii") << endl;
133 for (const auto& p : particles) {
134 ofs << static_cast<int32_t>(p.boundaryCondition) << endl;
135 }
136 ofs << dataArrayEnd() << endl;
137 // Fluid Type
138 ofs << dataArrayBegin("Int32", "Fluid Type", 1, "ascii") << endl;
139 for (const auto& p : particles) {
140 ofs << p.fluidType << endl;
141 }
142 ofs << dataArrayEnd() << endl;
143 ofs << "</PointData>" << endl;
144
148 ofs << "<Cells>" << endl;
149 // connectivity
150 ofs << dataArrayBegin("Int64", "connectivity", 1, "ascii") << endl;
151 for (size_t i = 0; i < particles.size(); i++) {
152 ofs << i << endl;
153 }
154 ofs << dataArrayEnd() << endl;
155 // offsets
156 ofs << dataArrayBegin("Int64", "offsets", 1, "ascii") << endl;
157 for (size_t i = 0; i < particles.size(); i++) {
158 ofs << i + 1 << endl;
159 }
160 ofs << dataArrayEnd() << endl;
161 // types
162 ofs << dataArrayBegin("UInt8", "types", 1, "ascii") << endl;
163 for (const auto& p : particles) {
164 ofs << 1 << endl; // 1: particle
165 }
166 ofs << dataArrayEnd() << endl;
167 ofs << "</Cells>" << endl;
168 ofs << "</Piece>" << endl;
169 // ---------------------
170 // ---- Field data ----
171 // ---------------------
172 ofs << "<FieldData>" << endl;
173 ofs << "<DataArray type=\"Float64\" Name=\"Time\" NumberOfTuples=\"1\" format=\"ascii\">" << endl;
174 ofs << time << endl;
175 ofs << "</DataArray>" << endl;
176 ofs << "</FieldData>" << endl;
177 ofs << "</UnstructuredGrid>" << endl;
178 ofs << "</VTKFile>" << endl;
179}
180
181void ParticlesExporter::toVtuBinary(const fs::path& path, const double& time, const double& n0ForNumberDensity) {
182 std::ofstream ofs(path);
183 if (ofs.fail()) {
184 cerr << "cannot write " << path << endl;
185 std::exit(-1);
186 }
187
188 std::stringstream binaryData; // binary data to be written
189
190 // header
191 ofs << "<?xml version='1.0' encoding='UTF-8'?>" << endl;
192 ofs << "<VTKFile type=\"UnstructuredGrid\" version=\"1.0\" header_type=\"UInt64\" byte_order=\"";
193 if (isBigEndian()) {
194 ofs << "BigEndian";
195 } else {
196 ofs << "LittleEndian";
197 }
198 ofs << "\">" << endl;
199 ofs << "<UnstructuredGrid>" << endl
200 << "<Piece NumberOfPoints=\"" << particles.size() << "\" NumberOfCells=\"" << particles.size() << "\">" << endl;
204 ofs << "<Points>" << endl;
205 ofs << dataArrayBegin("Float64", "Position", 3, "appended", binaryData.tellp()) << endl;
206 ofs << dataArrayEnd() << endl;
207 uint64_t length_points = particles.size() * 3 * sizeof(double);
208 binaryData.write(reinterpret_cast<char*>(&length_points), sizeof(uint64_t));
209 for (const auto& p : particles) {
210 double x = p.position.x();
211 double y = p.position.y();
212 double z = p.position.z();
213 binaryData.write(reinterpret_cast<char*>(&x), sizeof(double));
214 binaryData.write(reinterpret_cast<char*>(&y), sizeof(double));
215 binaryData.write(reinterpret_cast<char*>(&z), sizeof(double));
216 }
217 ofs << "</Points>" << endl;
218 // ---------------------
219 // ----- PointData -----
220 // ---------------------
221 ofs << "<PointData>" << endl;
222 // Particle Type
223 ofs << dataArrayBegin("Int32", "Particle Type", 1, "appended", binaryData.tellp()) << endl;
224 ofs << dataArrayEnd() << endl;
225 uint64_t length_particle_type = particles.size() * sizeof(int32_t);
226 binaryData.write(reinterpret_cast<char*>(&length_particle_type), sizeof(uint64_t));
227 for (const auto& p : particles) {
228 int32_t type = static_cast<int32_t>(p.type);
229 binaryData.write(reinterpret_cast<char*>(&type), sizeof(int32_t));
230 }
231 // Velocity
232 ofs << dataArrayBegin("Float64", "Velocity", 3, "appended", binaryData.tellp()) << endl;
233 ofs << dataArrayEnd() << endl;
234 uint64_t length_velocity = particles.size() * 3 * sizeof(double);
235 binaryData.write(reinterpret_cast<char*>(&length_velocity), sizeof(uint64_t));
236 for (const auto& p : particles) {
237 double vx = p.velocity.x();
238 double vy = p.velocity.y();
239 double vz = p.velocity.z();
240 binaryData.write(reinterpret_cast<char*>(&vx), sizeof(double));
241 binaryData.write(reinterpret_cast<char*>(&vy), sizeof(double));
242 binaryData.write(reinterpret_cast<char*>(&vz), sizeof(double));
243 }
244 // Pressure
245 ofs << dataArrayBegin("Float64", "Pressure", 1, "appended", binaryData.tellp()) << endl;
246 ofs << dataArrayEnd() << endl;
247 uint64_t length_pressure = particles.size() * sizeof(double);
248 binaryData.write(reinterpret_cast<char*>(&length_pressure), sizeof(uint64_t));
249 for (const auto& p : particles) {
250 double pressure = p.pressure;
251 binaryData.write(reinterpret_cast<char*>(&pressure), sizeof(double));
252 }
253 // Number Density
254 ofs << dataArrayBegin("Float64", "Number Density", 1, "appended", binaryData.tellp()) << endl;
255 ofs << dataArrayEnd() << endl;
256 uint64_t length_number_density = particles.size() * sizeof(double);
257 binaryData.write(reinterpret_cast<char*>(&length_number_density), sizeof(uint64_t));
258 for (const auto& p : particles) {
259 double number_density = p.numberDensity;
260 binaryData.write(reinterpret_cast<char*>(&number_density), sizeof(double));
261 }
262 // Number Density Ratio
263 ofs << dataArrayBegin("Float64", "Number Density Ratio", 1, "appended", binaryData.tellp()) << endl;
264 ofs << dataArrayEnd() << endl;
265 uint64_t length_number_density_ratio = particles.size() * sizeof(double);
266 binaryData.write(reinterpret_cast<char*>(&length_number_density_ratio), sizeof(uint64_t));
267 for (const auto& p : particles) {
268 double number_density_ratio = p.numberDensity / n0ForNumberDensity;
269 binaryData.write(reinterpret_cast<char*>(&number_density_ratio), sizeof(double));
270 }
271 // Boundary Condition
272 ofs << dataArrayBegin("Int32", "Boundary Condition", 1, "appended", binaryData.tellp()) << endl;
273 ofs << dataArrayEnd() << endl;
274 uint64_t length_boundary_condition = particles.size() * sizeof(int32_t);
275 binaryData.write(reinterpret_cast<char*>(&length_boundary_condition), sizeof(uint64_t));
276 for (const auto& p : particles) {
277 int32_t boundary_condition = static_cast<int32_t>(p.boundaryCondition);
278 binaryData.write(reinterpret_cast<char*>(&boundary_condition), sizeof(int32_t));
279 }
280 // Fluid Type
281 ofs << dataArrayBegin("Int32", "Fluid Type", 1, "appended", binaryData.tellp()) << endl;
282 ofs << dataArrayEnd() << endl;
283 uint64_t length_fluid_type = particles.size() * sizeof(int32_t);
284 binaryData.write(reinterpret_cast<char*>(&length_fluid_type), sizeof(uint64_t));
285 for (const auto& p : particles) {
286 int32_t fluid_type = p.fluidType;
287 binaryData.write(reinterpret_cast<char*>(&fluid_type), sizeof(int32_t));
288 }
289 ofs << "</PointData>" << endl;
293 ofs << "<Cells>" << endl;
294 // connectivity
295 ofs << dataArrayBegin("Int64", "connectivity", 1, "appended", binaryData.tellp()) << endl;
296 ofs << dataArrayEnd() << endl;
297 uint64_t length_connectivity = particles.size() * sizeof(int64_t);
298 binaryData.write(reinterpret_cast<char*>(&length_connectivity), sizeof(uint64_t));
299 for (size_t i = 0; i < particles.size(); i++) {
300 int64_t connectivity = i;
301 binaryData.write(reinterpret_cast<char*>(&connectivity), sizeof(int64_t));
302 }
303 // offsets
304 ofs << dataArrayBegin("Int64", "offsets", 1, "appended", binaryData.tellp()) << endl;
305 ofs << dataArrayEnd() << endl;
306 uint64_t length_offset = particles.size() * sizeof(int64_t);
307 binaryData.write(reinterpret_cast<char*>(&length_offset), sizeof(uint64_t));
308 for (size_t i = 0; i < particles.size(); i++) {
309 int64_t offset = i + 1;
310 binaryData.write(reinterpret_cast<char*>(&offset), sizeof(int64_t));
311 }
312 // types
313 ofs << dataArrayBegin("UInt8", "types", 1, "appended", binaryData.tellp()) << endl;
314 ofs << dataArrayEnd() << endl;
315 uint64_t length_types = particles.size() * sizeof(uint8_t);
316 binaryData.write(reinterpret_cast<char*>(&length_types), sizeof(uint64_t));
317 for (const auto& p : particles) {
318 uint8_t type = 1; // 1: particle
319 binaryData.write(reinterpret_cast<char*>(&type), sizeof(uint8_t));
320 }
321 ofs << "</Cells>" << endl;
322 ofs << "</Piece>" << endl;
323 // ---------------------
324 // ---- Field data ----
325 // ---------------------
326 ofs << "<FieldData>" << endl;
327 ofs << "<DataArray type=\"Float64\" Name=\"Time\" NumberOfTuples=\"1\" format=\"appended\" offset=\""
328 << binaryData.tellp() << "\"/>" << endl;
329 // Time
330 uint64_t length_time = sizeof(double);
331 binaryData.write(reinterpret_cast<char*>(&length_time), sizeof(uint64_t));
332 double time_copied = time; // copy time to use reinterpret_cast
333 binaryData.write(reinterpret_cast<char*>(&time_copied), sizeof(double));
334 ofs << "</FieldData>" << endl;
335 ofs << "</UnstructuredGrid>" << endl;
336 // ---------------------
337 // --- Appended Data ---
338 // ---------------------
339 ofs << "<AppendedData encoding=\"raw\">" << endl;
340 ofs << "_" << binaryData.str() << endl;
341 ofs << "</AppendedData>" << endl;
342 ofs << "</VTKFile>" << endl;
343}
344
345void ParticlesExporter::toVtu(
346 const std::filesystem::path& path, const double& time, const double& n0ForNumberDensity, const bool binary
347) {
348 if (binary) {
349 toVtuBinary(path, time, n0ForNumberDensity);
350 } else {
351 toVtuAscii(path, time, n0ForNumberDensity);
352 }
353}
354
355void ParticlesExporter::toCsv(const fs::path& path, const double& time) {
356 std::ofstream ofs(path);
357 if (ofs.fail()) {
358 cerr << "cannot write " << path << endl;
359 std::exit(-1);
360 }
361
362 ofs << time << endl;
363 ofs << particles.size() << endl;
364 ofs << "type,fluidType,x,y,z,vx,vy,vz" << endl;
365 for (const auto& p : particles) {
366 ofs << static_cast<int>(p.type) << ",";
367 ofs << p.fluidType << ",";
368 ofs << p.position.x() << "," << p.position.y() << "," << p.position.z() << ",";
369 ofs << p.velocity.x() << "," << p.velocity.y() << "," << p.velocity.z();
370 ofs << endl;
371 }
372}
ParticlesExporter::toVtuAscii
void toVtuAscii(const std::filesystem::path &path, const double &time, const double &n0ForNumberDensity=1.0)
Export the particles to a file in the VTK zlib format.
Definition particles_exporter.cpp:67
ParticlesExporter::setParticles
void setParticles(const Particles &particles)
Set the particles to export to a file. This method is required before exporting.
Definition particles_exporter.cpp:47
ParticlesExporter::toVtu
void toVtu(const std::filesystem::path &path, const double &time, const double &n0ForNumberDensity=1.0, const bool binary=false)
Export the particles to a file in the VTK zlib format.
Definition particles_exporter.cpp:345
ParticlesExporter::toVtuBinary
void toVtuBinary(const std::filesystem::path &path, const double &time, const double &n0ForNumberDensity=1.0)
Export the particles to a file in the VTK zlib format.
Definition particles_exporter.cpp:181
ParticlesExporter::toProf
void toProf(const std::filesystem::path &path, const double &time)
Export the particles to a file in the Prof format.
Definition particles_exporter.cpp:51
ParticlesExporter::dataArrayBegin
std::string dataArrayBegin(const std::string &type, const std::string &name, int numberOfComponents, const std::string &format, size_t offset=SIZE_MAX) const
write the header of DataArray for the VTK format
Definition particles_exporter.cpp:18
ParticlesExporter::dataArrayEnd
std::string dataArrayEnd() const
write the end of DataArray for the VTK format
Definition particles_exporter.cpp:43
ParticlesExporter::isBigEndian
bool isBigEndian() const
detect the endian of the system
Definition particles_exporter.cpp:12
ParticlesExporter::toCsv
void toCsv(const std::filesystem::path &path, const double &time)
Export the particles to a file in the CSV format.
Definition particles_exporter.cpp:355
Particles
A collection of particles.
Definition particles.hpp:10
Particles::size
int size() const
Get the number of particles.
Definition particles.cpp:21
particles_exporter.hpp