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