1/*
2 * VTKWriter_grids_st.hpp
3 *
4 * Created on: Sep 3, 2015
5 * Author: Pietro Incardona
6 */
7
8#ifndef SRC_VTKWRITER_GRIDS_ST_HPP_
9#define SRC_VTKWRITER_GRIDS_ST_HPP_
10
11
12#include <boost/mpl/pair.hpp>
13#include "VTKWriter_grids_util.hpp"
14#include "util/util_debug.hpp"
15#include "util/convert.hpp"
16
17/*! \brief for each combination in the cell grid you can have different grids
18 *
19 * \tparam Grid type of grid
20 *
21 */
22template <typename Grid>
23struct cell_grid
24{
25 //! vector of fused grids
26 openfpm::vector<const Grid *> grids;
27
28 //! combination
29 //! (used to calculate the grid shift from the starting point of the cell)
30 comb<Grid::dims> cmb;
31
32 cell_grid() {}
33
34 /*! \brief construct a cell grid
35 *
36 * \param cmb in which position this grid live
37 *
38 */
39 cell_grid(const comb<Grid::dims> & cmb)
40 :cmb(cmb)
41 {}
42
43 /*! \brief copy contructor
44 *
45 * \param cg element to copy
46 *
47 */
48 cell_grid(const cell_grid<Grid> & cg)
49 {
50 this->operator=(cg);
51 }
52
53 /*! \brief copy constructor
54 *
55 * \param cg element to copy
56 *
57 */
58 cell_grid(cell_grid<Grid> && cg)
59 {
60 this->operator=(cg);
61 }
62
63 /*! \brief Copy the cell grid
64 *
65 * \param cg cell_grid to copy
66 *
67 * \return itself
68 *
69 */
70 cell_grid<Grid> & operator=(const cell_grid<Grid> & cg)
71 {
72 cmb = cg.cmb;
73 grids = cg.grids;
74
75 return *this;
76 }
77
78 /*! \brief Copy the cell grid
79 *
80 * \param cg cell_grid to copy
81 *
82 * \return itself
83 *
84 */
85 cell_grid<Grid> & operator=(cell_grid<Grid> && cg)
86 {
87 cmb = cg.cmb;
88 grids = cg.grids;
89
90 return *this;
91 }
92};
93
94/*! \brief convert a staggered element into a string for vtk write
95 *
96 * \tparam Grid type of the grid
97 * \tparam St space type
98 *
99 */
100template <typename Grid, typename St>
101class ele_g_st
102{
103public:
104
105 //! grid type
106 typedef Grid value_type;
107
108 //! constructor
109 ele_g_st(){};
110
111 /*! \brief convert a staggered grid property into a string
112 *
113 * \param offset shift of the staggered element
114 * \param spacing of the grid
115 * \param dom Part of the grid that is real domain
116 *
117 */
118 ele_g_st(const Point<Grid::dims,St> & offset,
119 const Point<Grid::dims,St> & spacing,
120 const Box<Grid::dims,St> & dom)
121 :offset(offset),spacing(spacing),dom(dom)
122 {}
123
124 //! output string
125 std::string dataset;
126 //! fused grids
127 openfpm::vector<cell_grid<Grid>> g;
128 //! offset where it start the grid
129 Point<Grid::dims,St> offset;
130 //! spacing of the grid
131 Point<Grid::dims,St> spacing;
132 //! Part of the grid that is real domain
133 Box<Grid::dims,size_t> dom;
134
135 /*! \brief Copy constructor
136 *
137 * \param ele element to copy
138 *
139 */
140 inline ele_g_st(const ele_g_st & ele)
141 {
142 this->operator=(ele);
143 }
144
145 /*! \brief Copy constructor
146 *
147 * \param ele element to copy
148 *
149 */
150 inline ele_g_st(ele_g_st && ele)
151 {
152 this->operator=(ele);
153 }
154
155 /*! \brief Copy the object
156 *
157 * \param ele ele_g_st to copy
158 *
159 * \return itself
160 *
161 */
162 ele_g_st<Grid,St> & operator=(const ele_g_st & ele)
163 {
164 dataset = ele.dataset;
165 g = ele.g;
166 offset = ele.offset;
167 spacing = ele.spacing;
168 dom = ele.dom;
169
170 return *this;
171 }
172
173 /*! \brief Copy the object
174 *
175 * \param ele ele_g_st to copy
176 *
177 * \return itself
178 *
179 */
180 ele_g_st<Grid,St> & operator=(ele_g_st && ele)
181 {
182 dataset = ele.dataset;
183 g = ele.g;
184 offset = ele.offset;
185 spacing = ele.spacing;
186 dom = ele.dom;
187
188 return *this;
189 }
190};
191
192/*!
193 *
194 * It write a VTK format file in case of grids defined on a space
195 *
196 * \tparam boost::mpl::pair<G,S>
197 *
198 * where G is the type of grid S is the type of space, float, double ...
199 *
200 */
201template <typename pair>
202class VTKWriter<pair,VECTOR_ST_GRIDS>
203{
204 //! Vector of grids
205 openfpm::vector< ele_g_st<typename pair::first,typename pair::second> > vg;
206
207 /*! \brief Get the total number of points
208 *
209 * \return the total number
210 *
211 */
212 size_t get_total()
213 {
214 size_t tot = 0;
215
216 //! Calculate the full number of vertices
217 for (size_t i = 0 ; i < vg.size() ; i++)
218 {
219 for (size_t j = 0 ; j < vg.get(i).g.size() ; j++)
220 {
221 if (vg.get(i).g.get(j).grids.size() != 0)
222 tot += vg.get(i).g.get(j).grids.get(0)->size();
223 }
224 }
225 return tot;
226 }
227
228 /*! \brief It get the vertex properties list
229 *
230 * It get the vertex properties list of the vertex defined as VTK header
231 *
232 * \return a string that define the vertex properties in graphML format
233 *
234 */
235
236 std::string get_vertex_properties_list()
237 {
238 //! vertex property output string
239 std::string v_out;
240
241 // write the number of vertex
242 v_out += "VERTICES " + std::to_string(get_total()) + " " + std::to_string(get_total() * 2) + "\n";
243
244 // return the vertex properties string
245 return v_out;
246 }
247
248 /*! \brief It get the vertex properties list
249 *
250 * It get the vertex properties list of the vertex defined as a VTK header
251 *
252 * \return a string that define the vertex properties in graphML format
253 *
254 */
255
256 std::string get_point_properties_list()
257 {
258 //! vertex property output string
259 std::string v_out;
260
261 // write the number of vertex
262 v_out += "POINTS " + std::to_string(get_total()) + " float" + "\n";
263
264 // return the vertex properties string
265 return v_out;
266 }
267
268 /*! \brief Create the VTK point definition
269 *
270 * \return the list of points
271 *
272 */
273 std::string get_point_list()
274 {
275 //! vertex node output string
276 std::stringstream v_out;
277
278 //! For each sub-domain
279 for (size_t i = 0 ; i < vg.size() ; i++)
280 {
281 // For each position in the cell
282 for (size_t j = 0 ; j < vg.get(i).g.size() ; j++)
283 {
284 // If there are no grid in this position
285 if (vg.get(i).g.get(j).grids.size() == 0)
286 continue;
287
288 //! Get the iterator
289 auto it = vg.get(i).g.get(j).grids.get(0)->getIterator();
290
291 //! Where the grid is defined
292 Box<pair::first::dims,typename pair::second> dom;
293
294 // Calculate the offset of the grid considering its cell position
295 Point<pair::first::dims,typename pair::second> middle = vg.get(i).spacing / 2;
296 Point<pair::first::dims,typename pair::second> one;
297 one.one();
298 one = one + toPoint<pair::first::dims,typename pair::second>::convert(vg.get(i).g.get(j).cmb);
299 Point<pair::first::dims,typename pair::second> offset = pmul(middle,one) + vg.get(i).offset;
300
301 // if there is the next element
302 while (it.isNext())
303 {
304 Point<pair::first::dims,typename pair::second> p;
305 p = it.get().toPoint();
306 p = pmul(p,vg.get(i).spacing) + offset;
307
308 if (pair::first::dims == 2)
309 v_out << p.toString() << " 0.0" << "\n";
310 else
311 v_out << p.toString() << "\n";
312
313 // increment the iterator
314 ++it;
315 }
316 }
317 }
318
319 // return the vertex list
320 return v_out.str();
321 }
322
323 /*! \brief It generate a name for the property cell component
324 *
325 * \param k component in the cell
326 *
327 * \return property name
328 *
329 */
330 std::string get_prop_components(size_t k)
331 {
332 std::stringstream v_out;
333
334 //! For each sub-domain
335 for (size_t i = 0 ; i < vg.size() ; i++)
336 {
337 // For each position in the cell
338 for (size_t j = 0 ; j < vg.get(i).g.size() ; j++)
339 {
340 if (k < vg.get(i).g.get(j).grids.size())
341 {
342 // get the combination string
343 v_out << vg.get(i).g.get(j).cmb.to_string();
344 }
345 }
346 }
347
348 return v_out.str();
349 }
350
351 /*! \brief Create the VTK properties output
352 *
353 * \param k component
354 * \param prop_name property name
355 *
356 * \return the property output string for the grid
357 *
358 */
359 std::string get_properties_output(size_t k, std::string prop_name)
360 {
361 //! vertex node output string
362 std::stringstream v_out;
363
364 // Check if T is a supported format
365 // for now we support only scalar of native type
366
367 typedef typename boost::mpl::at<typename pair::first::value_type::type,boost::mpl::int_<0>>::type ctype;
368
369 std::string type = getType<ctype>();
370
371 // if the type is not supported return
372 if (type.size() == 0)
373 {
374#ifndef DISABLE_ALL_RTTI
375 std::cerr << "Error " << __FILE__ << ":" << __LINE__ << " the type " << demangle(typeid(ctype).name()) << " is not supported by vtk\n";
376#endif
377 return "";
378 }
379
380 std::string prp_cp = get_prop_components(k);
381
382 // Create point data properties
383 v_out << "SCALARS " << prop_name << "_" << prp_cp << " " << type + "\n";
384
385 // Default lookup table
386 v_out << "LOOKUP_TABLE default\n";
387
388 //! For each sub-domain
389 for (size_t i = 0 ; i < vg.size() ; i++)
390 {
391 // For each position in the cell
392 for (size_t j = 0 ; j < vg.get(i).g.size() ; j++)
393 {
394 // If there are no grid in this position
395 if (vg.get(i).g.get(j).grids.size() == 0)
396 continue;
397
398 if (k < vg.get(i).g.get(j).grids.size())
399 {
400 // Grid source
401 auto & g_src = *vg.get(i).g.get(j).grids.get(k);
402
403 //! Get the iterator
404 auto it = g_src.getIterator();
405
406 //! Where the grid is defined
407 Box<pair::first::dims,typename pair::second> dom;
408
409 // if there is the next element
410 while (it.isNext())
411 {
412 auto key = it.get();
413
414 v_out << std::to_string(g_src.template get<0>(key)) << "\n";
415
416 // increment the iterator
417 ++it;
418 }
419 }
420 else
421 {
422 // Grid source
423 auto & g_src = *vg.get(i).g.get(j).grids.get(0);
424
425 //! Get the iterator
426 auto it = g_src.getIterator();
427
428 //! Where the grid is defined
429 Box<pair::first::dims,typename pair::second> dom;
430
431 // if there is the next element
432 while (it.isNext())
433 {
434 v_out << "0\n";
435
436 // increment the iterator
437 ++it;
438 }
439 }
440 }
441 }
442
443 // return the vertex list
444 return v_out.str();
445 }
446
447 /*! \brief Return the output of the domain property
448 *
449 * \return vtk output
450 *
451 */
452 std::string lastProp()
453 {
454 //! vertex node output string
455 std::stringstream v_out;
456
457 // Create point data properties
458 v_out << "SCALARS domain float\n";
459
460 // Default lookup table
461 v_out << "LOOKUP_TABLE default\n";
462
463 //! For each sub-domain
464 for (size_t i = 0 ; i < vg.size() ; i++)
465 {
466 // For each position in the cell
467 for (size_t j = 0 ; j < vg.get(i).g.size() ; j++)
468 {
469 // If there are no grid in this position
470 if (vg.get(i).g.get(j).grids.size() == 0)
471 continue;
472
473 //! Get the iterator
474 auto it = vg.get(i).g.get(j).grids.get(0)->getIterator();
475
476 // if there is the next element
477 while (it.isNext())
478 {
479 if (vg.get(i).dom.isInside(it.get().toPoint()) == true)
480 v_out << "1.0\n";
481 else
482 v_out << "0.0\n";
483
484 // increment the iterator and counter
485 ++it;
486 }
487 }
488 }
489
490 return v_out.str();
491 }
492
493 /*! \brief Get the maximum number of fused grid
494 *
495 * \return the maximum number of fused grids
496 *
497 */
498 size_t getMaxFused()
499 {
500 size_t max = 0;
501
502 //! For each sub-domain
503 for (size_t i = 0 ; i < vg.size() ; i++)
504 {
505 // For each position in the cell
506 for (size_t j = 0 ; j < vg.get(i).g.size() ; j++)
507 {
508 // If there are no grid in this position
509 if (vg.get(i).g.get(j).grids.size() > max)
510 max = vg.get(i).g.get(j).grids.size();
511 }
512 }
513
514 return max;
515 }
516
517 /*! \brief Create the VTK vertex definition
518 *
519 * \return the string with the vertices as string
520 *
521 */
522 std::string get_vertex_list()
523 {
524 //! vertex node output string
525 std::string v_out;
526
527 size_t k = 0;
528
529 //! For each sub-domain
530 for (size_t i = 0 ; i < vg.size() ; i++)
531 {
532 // For each position in the cell
533 for (size_t j = 0 ; j < vg.get(i).g.size() ; j++)
534 {
535 // If there are no grid in this position
536 if (vg.get(i).g.get(j).grids.size() == 0)
537 continue;
538 //! For each grid point create a vertex
539 auto it = vg.get(i).g.get(j).grids.get(0)->getIterator();
540
541 while (it.isNext())
542 {
543 v_out += "1 " + std::to_string(k) + "\n";
544
545 ++k;
546 ++it;
547 }
548 }
549 }
550 // return the vertex list
551 return v_out;
552 }
553
554 /*! \brief Get the point data header
555 *
556 * \return a string with the point data header for VTK format
557 *
558 */
559
560 std::string get_point_data_header()
561 {
562 std::string v_out;
563
564 v_out += "POINT_DATA " + std::to_string(get_total()) + "\n";
565
566 return v_out;
567 }
568
569 /*! \brief Append the grid to the sub-domain, if for a sub-domain we have a grid that is overlapping
570 * fuse them, otherwise create a new combination and grid
571 *
572 * \param id sub-domain id
573 * \param g grid to output
574 * \param cmb position of the grid
575 *
576 * \return a valid slot, if does not exist it append the grid at the end with the new combination
577 *
578 */
579 void append_grid(size_t id, const typename pair::first & g, const comb<pair::first::dims> & cmb)
580 {
581 for(size_t i = 0 ; i < vg.get(id).g.size() ; i++)
582 {
583 // for each defined grid if exist the combination fuse
584 if (cmb == vg.get(id).g.get(i).cmb)
585 {
586 vg.get(id).g.get(i).grids.add(&g);
587 return;
588 }
589 }
590
591 // if the combination does not exist add the grid
592 cell_grid< typename pair::first> cg(cmb);
593 vg.get(id).g.add(cg);
594 vg.get(id).g.last().grids.add(&g);
595 }
596
597public:
598
599 /*!
600 *
601 * VTKWriter constructor
602 *
603 */
604 VTKWriter()
605 {}
606
607 /*! \brief Add grid dataset
608 *
609 * \param i sub-domain id
610 * \param g Grid to add
611 * \param offset grid offset
612 * \param spacing spacing of the grid
613 * \param dom part of the spacethat is the domain
614 * \param cmb position of the grid
615 *
616 */
617 void add(size_t i,
618 const typename pair::first & g,
619 const Point<pair::first::dims,typename pair::second> & offset,
620 const Point<pair::first::dims,typename pair::second> & spacing,
621 const Box<pair::first::dims,typename pair::second> & dom,
622 const comb<pair::first::dims> & cmb)
623 {
624 //! Increase the size
625 if (i >= vg.size())
626 vg.resize(i+1);
627
628 vg.get(i).offset = offset;
629 vg.get(i).spacing = spacing;
630 vg.get(i).dom = dom;
631
632 // append the grid
633 append_grid(i,g,cmb);
634 }
635
636 /*! \brief It write a VTK file from a graph
637 *
638 * \tparam prp_out which properties to output [default = -1 (all)]
639 *
640 * \param file path where to write
641 * \param g_name of the set of grids
642 * \param ft specify if it is a VTK BINARY or ASCII file [default = ASCII]
643 *
644 * \return true if the file is succeful written
645 *
646 */
647
648 template<int prp = -1> bool write(std::string file,
649 std::string g_name = "grids",
650 file_type ft = file_type::ASCII)
651 {
652 // Header for the vtk
653 std::string vtk_header;
654 // Point list of the VTK
655 std::string point_list;
656 // Vertex list of the VTK
657 std::string vertex_list;
658 // Graph header
659 std::string vtk_binary_or_ascii;
660 // vertex properties header
661 std::string point_prop_header;
662 // edge properties header
663 std::string vertex_prop_header;
664 // Data point header
665 std::string point_data_header;
666 // Data point
667 std::string point_data;
668
669 // VTK header
670 vtk_header = "# vtk DataFile Version 3.0\n"
671 + g_name + "\n";
672
673 // Choose if binary or ASCII
674 if (ft == file_type::ASCII)
675 {vtk_header += "ASCII\n";}
676 else
677 {vtk_header += "BINARY\n";}
678
679 // Data type for graph is DATASET POLYDATA
680 vtk_header += "DATASET POLYDATA\n";
681
682 // point properties header
683 point_prop_header = get_point_properties_list();
684
685 // Get point list
686 point_list = get_point_list();
687
688 // vertex properties header
689 vertex_prop_header = get_vertex_properties_list();
690
691 // Get vertex list
692 vertex_list = get_vertex_list();
693
694 // Get the point data header
695 point_data_header = get_point_data_header();
696
697 // Get the maximum number of fused grids
698 size_t mf = getMaxFused();
699
700 // For each property in the vertex type produce a point data
701 for (size_t i = 0 ; i < mf ; i++)
702 point_data += get_properties_output(i,g_name);
703
704 lastProp();
705
706
707 // write the file
708 std::ofstream ofs(file);
709
710 // Check if the file is open
711 if (ofs.is_open() == false)
712 {std::cerr << "Error cannot create the VTK file: " + file + "\n";}
713
714 ofs << vtk_header << point_prop_header << point_list <<
715 vertex_prop_header << vertex_list << point_data_header << point_data;
716
717 // Close the file
718
719 ofs.close();
720
721 // Completed succefully
722 return true;
723 }
724};
725
726
727#endif /* SRC_VTKWRITER_GRIDS_ST_HPP_ */
728