LehrFEM++ 1.0.0
A simple Finite Element Library for teaching
Loading...
Searching...
No Matches
vtk_writer.cc
1
9#include "vtk_writer.h"
10
11#include <lf/base/base.h>
12
13#include <Eigen/Eigen>
14#include <boost/fusion/include/adapt_struct.hpp>
15#include <boost/phoenix/phoenix.hpp>
16#include <boost/phoenix/scope/let.hpp>
17#include <boost/spirit/include/karma.hpp>
18#include <fstream>
19#include <unsupported/Eigen/KroneckerProduct>
20
21#include "eigen_fusion_adapter.h"
22
23template <class RESULT_TYPE, class LAMBDA>
25 public:
26 using result_type = RESULT_TYPE;
27
28 template <class T>
30 return lambda_(std::forward<T>(arg));
31 }
32
33 explicit VariantVisitor(LAMBDA lambda) : lambda_(lambda) {}
34
35 private:
36 LAMBDA lambda_;
37};
38
39template <class RESULT_TYPE, class LAMBDA>
40VariantVisitor<RESULT_TYPE, LAMBDA> make_VariantVisitor(LAMBDA l) {
42}
43
44// clang-format off
45BOOST_FUSION_ADAPT_STRUCT(lf::io::VtkFile::UnstructuredGrid,
46 (std::vector<Eigen::Vector3f>, points)
47 (std::vector<std::vector<unsigned int>>, cells)
48 (std::vector<lf::io::VtkFile::CellType>, cell_types)
49);
50
51BOOST_FUSION_ADAPT_TPL_STRUCT((T), (lf::io::VtkFile::FieldDataArray)(T),
52 (std::string, name),
53 (auto, data)
54);
55
56BOOST_FUSION_ADAPT_TPL_STRUCT((T), (lf::io::VtkFile::ScalarData)(T),
57 (std::string, name)
58 (std::string, lookup_table)
59 (auto, data)
60);
61
62BOOST_FUSION_ADAPT_TPL_STRUCT((T), (lf::io::VtkFile::VectorData)(T),
63 (std::string, name)
64 (auto, data)
65);
66
67// BOOST_FUSION_ADAPT_STRUCT(lf::io::VtkFile::Attributes,
68// (auto, data)
69// );
70
71 // BOOST_FUSION_ADAPT_STRUCT(lf::io::VtkFile::Attributes,
72 // (std::vector<lf::io::VtkFile::ScalarData<int>>, scalar_int_data)
73 // (std::vector<lf::io::VtkFile::ScalarData<unsigned int>>, scalar_unsigned_int_data)
74 // (std::vector<lf::io::VtkFile::ScalarData<long>>, scalar_long_data)
75 // (std::vector<lf::io::VtkFile::ScalarData<unsigned long>>, scalar_unsigned_long_data)
76 // (std::vector<lf::io::VtkFile::ScalarData<float>>, scalar_float_data)
77 // (std::vector<lf::io::VtkFile::ScalarData<double>>, scalar_double_data)
78 // );
79
80BOOST_FUSION_ADAPT_STRUCT(lf::io::VtkFile,
81 (std::string, header)
83 (lf::io::VtkFile::UnstructuredGrid, unstructured_grid)
84 (lf::io::VtkFile::FieldData, field_data)
85 (lf::io::VtkFile::Attributes, point_data)
87);
88
89// clang-format on
90
91namespace /*anonymous*/ {
92void cell_list_size(unsigned int& result, // NOLINT
93 const std::vector<std::vector<unsigned int>>& cells) {
94 result = cells.size();
95 for (const auto& v : cells) {
96 result += v.size();
97 }
98}
99} // namespace
100
101// NOLINTNEXTLINE
102BOOST_PHOENIX_ADAPT_FUNCTION(void, CellListSize, cell_list_size, 2);
103
104namespace lf::io {
105std::ostream& operator<<(std::ostream& stream, const VtkFile::Format& f) {
106 switch (f) {
108 stream << "ASCII";
109 break;
111 stream << "BINARY";
112 break;
113 default:
114 throw base::LfException("Unknown VtkFileFormat.");
115 }
116 return stream;
117}
118
119namespace /*anonymous */ {
120namespace karma = boost::spirit::karma;
121template <class Iterator, bool BINARY>
122struct VtkGrammar : karma::grammar<Iterator, VtkFile()> {
123 // template <class T>
124 // karma::rule<Iterator, std::vector<T>()> vector_size() {
125 // return karma::int_[karma::_1 = boost::phoenix::size(karma::_val)];
126 // }
127
128 // NOLINTNEXTLINE
129 VtkGrammar() : VtkGrammar::base_type(root) {
130 using boost::phoenix::at_c;
131 using boost::phoenix::size;
132 using karma::_1;
133 using karma::_a;
134 using karma::_r1;
135 using karma::_val;
136 using karma::eol;
137 using karma::eps;
138 using karma::lit;
139
140 if (BINARY) {
141 points_vector %= *(karma::big_bin_float << karma::big_bin_float
142 << karma::big_bin_float)
143 << eol;
144 } else {
145 points_vector %= *(karma::float_ << ' ' << karma::float_ << ' '
146 << karma::float_ << karma::eol);
147 }
148
149 points = lit("POINTS ")
150 << karma::uint_[karma::_1 = boost::phoenix::size(karma::_val)]
151 << lit(" float") << karma::eol
152 << points_vector[karma::_1 = karma::_val];
153
154 // cells_vector %= (karma::uint_ % " ") << karma::eol;
155 if (BINARY) {
156 cells_vector %= karma::big_dword(boost::phoenix::size(karma::_val))
157 << *(karma::big_dword);
158 } else {
159 // cells_vector %= karma::duplicate
160 // [karma::uint_[karma::_1 = boost::phoenix::size(karma::_val)]
161 // << lit(' ') << (karma::uint_ % " ") << karma::eol];
162 cells_vector %= lit(size(_val))
163 << lit(' ') << (karma::uint_ % " ") << karma::eol;
164 }
165 cells = lit("CELLS ")
166 << karma::uint_[karma::_1 = boost::phoenix::size(karma::_val)]
167 << lit(' ') << karma::uint_[CellListSize(karma::_1, karma::_val)]
168 << karma::eol << (*cells_vector)[karma::_1 = karma::_val] << eol;
169
170 if (BINARY) {
171 cell_type =
172 karma::big_dword[karma::_1 =
173 boost::phoenix::static_cast_<int>(karma::_val)];
174 } else {
175 cell_type = karma::int_[karma::_1 = boost::phoenix::static_cast_<int>(
176 karma::_val)]
177 << karma::eol;
178 }
179 cell_types = lit("CELL_TYPES ")
180 << karma::uint_[karma::_1 = boost::phoenix::size(karma::_val)]
181 << karma::eol << (*cell_type)[karma::_1 = karma::_val];
182
183 unstructured_grid %= lit("DATASET UNSTRUCTURED_GRID")
184 << karma::eol << points << cells << cell_types;
185
186 // Field data
187 if (BINARY) {
188 field_data_int %= karma::string << lit(" 1 ") << lit(size(at_c<1>(_val)))
189 << lit(" int") << eol << *karma::big_dword
190 << eol;
191 field_data_float %= karma::string
192 << lit(" 1 ") << lit(size(at_c<1>(_val)))
193 << lit(" float") << eol << *(karma::big_bin_float)
194 << eol;
195
196 field_data_double %= karma::string
197 << lit(" 1 ") << lit(size(at_c<1>(_val)))
198 << lit(" double") << eol << *karma::big_bin_double
199 << eol;
200 } else {
201 field_data_int %= karma::string << lit(" 1 ") << lit(size(at_c<1>(_val)))
202 << lit(" int") << eol
203 << (karma::int_ % lit(' ')) << eol;
204 field_data_float %= karma::string
205 << lit(" 1 ") << lit(size(at_c<1>(_val)))
206 << lit(" float") << eol << (karma::float_ % lit(' '))
207 << eol;
208
209 field_data_double %= karma::string << lit(" 1 ")
210 << lit(size(at_c<1>(_val)))
211 << lit(" double") << eol
212 << (karma::double_ % lit(' ')) << eol;
213 }
214 field_data %= lit("FIELD FieldData ")
215 << lit(size(_val)) << eol
216 << *(field_data_int | field_data_float | field_data_double);
217
218 // Point/Cell data
219 if (BINARY) {
220 scalar_data_char %=
221 (lit("SCALARS ") << karma::string << lit(" char 1") << eol
222 << lit("LOOKUP_TABLE ") << karma::string << eol
223 << *karma::byte_)
224 << eol;
225
226 scalar_data_uchar %=
227 (lit("SCALARS ") << karma::string << lit(" char 1") << eol
228 << lit("LOOKUP_TABLE ") << karma::string << eol
229 << *karma::byte_)
230 << eol;
231
232 scalar_data_short %=
233 (lit("SCALARS ") << karma::string << lit(" short 1") << eol
234 << lit("LOOKUP_TABLE ") << karma::string << eol
235 << *karma::big_word)
236 << eol;
237
238 scalar_data_ushort %=
239 (lit("SCALARS ") << karma::string << lit(" unsigned_short 1") << eol
240 << lit("LOOKUP_TABLE ") << karma::string << eol
241 << *karma::big_word)
242 << eol;
243
244 scalar_data_int %=
245 (lit("SCALARS ") << karma::string << lit(" int 1") << eol
246 << lit("LOOKUP_TABLE ") << karma::string << eol
247 << *karma::big_dword)
248 << eol;
249 scalar_data_unsigned_int %=
250 (lit("SCALARS ") << karma::string << lit(" unsigned_int 1") << eol
251 << lit("LOOKUP_TABLE ") << karma::string << eol
252 << *karma::big_dword)
253 << eol;
254
255 // NOLINTNEXTLINE
256 if constexpr (sizeof(long) == 8) {
257 scalar_data_long %=
258 (lit("SCALARS ")
259 << karma::string << lit(" long 1") << eol << lit("LOOKUP_TABLE ")
260 << karma::string << eol << *karma::big_qword)
261 << eol;
262 scalar_data_unsigned_long %=
263 (lit("SCALARS ") << karma::string << lit(" unsigned_long 1") << eol
264 << lit("LOOKUP_TABLE ") << karma::string << eol
265 << *karma::big_qword)
266 << eol;
267 } else {
268 scalar_data_long %=
269 (lit("SCALARS ")
270 << karma::string << lit(" long 1") << eol << lit("LOOKUP_TABLE ")
271 << karma::string << eol << *karma::big_dword)
272 << eol;
273 scalar_data_unsigned_long %=
274 (lit("SCALARS ") << karma::string << lit(" unsigned_long 1") << eol
275 << lit("LOOKUP_TABLE ") << karma::string << eol
276 << *karma::big_dword)
277 << eol;
278 }
279
280 scalar_data_float %=
281 (lit("SCALARS ") << karma::string << lit(" float 1") << eol
282 << lit("LOOKUP_TABLE ") << karma::string << eol
283 << *karma::big_bin_float)
284 << eol;
285
286 scalar_data_double %=
287 (lit("SCALARS ") << karma::string << lit(" double 1") << eol
288 << lit("LOOKUP_TABLE ") << karma::string << eol
289 << *karma::big_bin_double)
290 << eol;
291
292 vector_data_float %= lit("VECTORS ")
293 << karma::string << lit(" float") << eol
294 << *(karma::big_bin_float << karma::big_bin_float
295 << karma::big_bin_float)
296 << eol;
297
298 vector_data_double %= lit("VECTORS ")
299 << karma::string << lit(" double") << eol
300 << *(karma::big_bin_double << karma::big_bin_double
301 << karma::big_bin_double)
302 << eol;
303 } else {
304 scalar_data_char %=
305 (lit("SCALARS ") << karma::string << lit(" char 1") << eol
306 << lit("LOOKUP_TABLE ") << karma::string << eol
307 << (*(karma::short_ % eol)))
308 << eol;
309
310 scalar_data_uchar %=
311 (lit("SCALARS ") << karma::string << lit(" char 1") << eol
312 << lit("LOOKUP_TABLE ") << karma::string << eol
313 << (*(karma::ushort_ % eol)))
314 << eol;
315
316 scalar_data_short %=
317 (lit("SCALARS ") << karma::string << lit(" short 1") << eol
318 << lit("LOOKUP_TABLE ") << karma::string << eol
319 << (*(karma::short_ % eol)))
320 << eol;
321
322 scalar_data_ushort %=
323 (lit("SCALARS ") << karma::string << lit(" unsigned_short 1") << eol
324 << lit("LOOKUP_TABLE ") << karma::string << eol
325 << (*(karma::ushort_ % eol)))
326 << eol;
327
328 scalar_data_int %=
329 (lit("SCALARS ") << karma::string << lit(" int 1") << eol
330 << lit("LOOKUP_TABLE ") << karma::string << eol
331 << (*(karma::int_ % eol)))
332 << eol;
333 scalar_data_unsigned_int %=
334 (lit("SCALARS ") << karma::string << lit(" unsigned_int 1") << eol
335 << lit("LOOKUP_TABLE ") << karma::string << eol
336 << (*(karma::uint_ % eol)))
337 << eol;
338
339 scalar_data_long %=
340 (lit("SCALARS ") << karma::string << lit(" long 1") << eol
341 << lit("LOOKUP_TABLE ") << karma::string << eol
342 << (*(karma::long_ % eol)))
343 << eol;
344
345 scalar_data_unsigned_long %=
346 (lit("SCALARS ") << karma::string << lit(" unsigned_long 1") << eol
347 << lit("LOOKUP_TABLE ") << karma::string << eol
348 << (*(karma::ulong_ % eol)))
349 << eol;
350
351 scalar_data_float %=
352 (lit("SCALARS ") << karma::string << lit(" float 1") << eol
353 << lit("LOOKUP_TABLE ") << karma::string << eol
354 << (*(karma::float_ % eol)))
355 << eol;
356
357 scalar_data_double %=
358 (lit("SCALARS ") << karma::string << lit(" double 1") << eol
359 << lit("LOOKUP_TABLE ") << karma::string << eol
360 << (*(karma::double_ % eol)))
361 << eol;
362
363 vector_data_float %= lit("VECTORS ")
364 << karma::string << lit(" float") << eol
365 << *(karma::float_ << lit(' ') << karma::float_
366 << lit(' ') << karma::float_
367 << eol);
368
369 vector_data_double %= lit("VECTORS ")
370 << karma::string << lit(" double") << eol
371 << *(karma::double_ << lit(' ') << karma::double_
372 << lit(' ') << karma::double_
373 << eol);
374 }
375
376 attributes %=
377 (*(scalar_data_char | scalar_data_uchar | scalar_data_short |
378 scalar_data_ushort | scalar_data_int | scalar_data_unsigned_int |
379 scalar_data_long | scalar_data_unsigned_long | scalar_data_float |
380 scalar_data_double | vector_data_float | vector_data_double));
381
382 root %= lit("# vtk DataFile Version 3.0")
383 << karma::eol << karma::string << karma::eol << karma::stream
384 << karma::eol << unstructured_grid << karma::eol << field_data
385 << lit("POINT_DATA ") << lit(size(at_c<0>(at_c<2>(_val)))) << eol
386 << attributes << lit("CELL_DATA ")
387 << lit(size(at_c<1>(at_c<2>(_val)))) << eol << attributes;
388 }
389
390 karma::rule<Iterator, VtkFile()> root;
391 karma::rule<Iterator, VtkFile()> header;
392 karma::rule<Iterator, VtkFile::UnstructuredGrid()> unstructured_grid;
393 karma::rule<Iterator, std::vector<Eigen::Vector3f>()> points;
394 karma::rule<Iterator, std::vector<Eigen::Vector3f>()> points_vector;
395 karma::rule<Iterator, std::vector<std::vector<unsigned int>>> cells;
396 karma::rule<Iterator, std::vector<unsigned int>> cells_vector;
397 karma::rule<Iterator, std::vector<VtkFile::CellType>> cell_types;
398 karma::rule<Iterator, VtkFile::CellType> cell_type;
399
400 karma::rule<Iterator, VtkFile::FieldData()> field_data;
401 karma::rule<Iterator, VtkFile::FieldDataArray<int>()> field_data_int;
402 karma::rule<Iterator, VtkFile::FieldDataArray<float>()> field_data_float;
403 karma::rule<Iterator, VtkFile::FieldDataArray<double>()> field_data_double;
404
405 karma::rule<Iterator, VtkFile::Attributes()> attributes;
406
407 karma::rule<Iterator, VtkFile::ScalarData<char>()> scalar_data_char;
408 karma::rule<Iterator, VtkFile::ScalarData<unsigned char>()> scalar_data_uchar;
409 karma::rule<Iterator, VtkFile::ScalarData<short>()> scalar_data_short;
410 karma::rule<Iterator, VtkFile::ScalarData<unsigned short>()>
411 scalar_data_ushort;
412 karma::rule<Iterator, VtkFile::ScalarData<int>()> scalar_data_int;
413 karma::rule<Iterator, VtkFile::ScalarData<unsigned int>()>
414 scalar_data_unsigned_int;
415 karma::rule<Iterator, VtkFile::ScalarData<long>()> scalar_data_long;
416 karma::rule<Iterator, VtkFile::ScalarData<unsigned long>()>
417 scalar_data_unsigned_long;
418 karma::rule<Iterator, VtkFile::ScalarData<float>()> scalar_data_float;
419 karma::rule<Iterator, VtkFile::ScalarData<double>()> scalar_data_double;
420 karma::rule<Iterator, VtkFile::VectorData<float>()> vector_data_float;
421 karma::rule<Iterator, VtkFile::VectorData<double>()> vector_data_double;
422};
423
425void ValidateVtkFile(const VtkFile& vtk_file) {
426 // check that header doesn't contain a new line character:
427 if (vtk_file.header.find('\n') != std::string::npos) {
428 throw base::LfException(
429 "Header of vtk file contains a new line character, this is not "
430 "allowed");
431 }
432 if (vtk_file.header.size() > 256) {
433 throw base::LfException("Header of vtk file is longer than 256 characters");
434 }
435 if (vtk_file.unstructured_grid.cell_types.size() !=
436 vtk_file.unstructured_grid.cells.size()) {
437 throw base::LfException(
438 "Mismatch of size of cell_types and cells in VtkFile.");
439 }
440 for (const auto& d : vtk_file.point_data) {
441 boost::apply_visitor(
442 [&](auto e) {
443 if (e.data.size() != vtk_file.unstructured_grid.points.size()) {
444 throw base::LfException("Length of dataset " + e.name +
445 " does not match the number of points.");
446 }
447 },
448 d);
449 }
450 for (const auto& d : vtk_file.cell_data) {
451 boost::apply_visitor(
452 [&](auto e) {
453 if (e.data.size() != vtk_file.unstructured_grid.cells.size()) {
454 throw base::LfException("Length of dataset " + e.name +
455 " does not match the number of cells.");
456 }
457 },
458 d);
459 }
460}
461
463unsigned int NumAuxNodes(base::RefEl ref_el, unsigned char order) {
464 switch (ref_el) {
465 case base::RefEl::kPoint():
466 return 1;
468 return order - 1;
469 case base::RefEl::kTria():
470 return (2 - 3 * order + order * order) / 2;
471 case base::RefEl::kQuad():
472 return (order - 1) * (order - 1);
473 default:
474 LF_VERIFY_MSG(false, "RefElType " << ref_el << " not supported.");
475 }
476}
477
478// Compute aux reference coordinates of lagrange points on segment
479Eigen::MatrixXd AuxNodesSegment(unsigned char order) {
480 LF_ASSERT_MSG(order > 1, "order <= 1");
481 return Eigen::VectorXd::LinSpaced(order - 1, 1. / (order),
482 (order - 1.) / order)
483 .transpose();
484}
485
486// compute aux reference coordinates of lagrange points on Tria:
487// NOLINTNEXTLINE(misc-no-recursion)
488Eigen::MatrixXd AuxNodesTria(unsigned char order) {
489 if (order < 3) {
490 return {};
491 }
492 Eigen::MatrixXd result(2, NumAuxNodes(base::RefEl::kTria(), order));
493 if (order == 3) {
494 result << 1. / 3., 1. / 3.;
495 } else if (order == 4) {
496 result << 0.25, 0.5, 0.25, 0.25, 0.25, 0.5;
497 } else {
498 // assign nodes:
499 result(0, 0) = 1. / order;
500 result(1, 0) = 1. / order;
501 result(0, 1) = (order - 2.) / order;
502 result(1, 1) = 1. / order;
503 result(0, 2) = 1. / order;
504 result(1, 2) = (order - 2.) / order;
505
506 if (order > 4) {
507 auto segment_points = AuxNodesSegment(order - 3).eval();
508 // assign edges:
509 result.block(0, 3, 2, order - 4) =
510 Eigen::Vector2d::UnitX() * segment_points * (order - 3.) / order +
511 Eigen::Vector2d(1. / order, 1. / order) *
512 Eigen::MatrixXd::Ones(1, order - 4);
513 result.block(0, order - 1, 2, order - 4) =
514 Eigen::Vector2d(-1, 1) * (order - 3.) / order * segment_points +
515 Eigen::Vector2d((order - 2.) / order, 1. / order) *
516 Eigen::MatrixXd::Ones(1, order - 4);
517 result.block(0, 2 * order - 5, 2, order - 4) =
518 -Eigen::Vector2d::UnitY() * (order - 3.) / order * segment_points +
519 Eigen::Vector2d(1. / order, (order - 2.) / order) *
520 Eigen::MatrixXd::Ones(1, order - 4);
521 }
522 if (order > 5) {
523 // assign interior points recursively:
524 auto points = AuxNodesTria(order - 3);
525 result.block(0, 3 * order - 9, 2, points.cols()) =
526 AuxNodesTria(order - 3) * (order - 6.) / order +
527 Eigen::Vector2d(2. / order, 2. / order) *
528 Eigen::MatrixXd::Ones(1, points.cols());
529 }
530 }
531 return result;
532}
533
534// compute aux reference coordinates of lagrange points on quad:
535Eigen::MatrixXd AuxNodesQuad(unsigned char order) {
536 Eigen::MatrixXd result(2, NumAuxNodes(base::RefEl::kQuad(), order));
537 result.row(0) =
538 Eigen::VectorXd::LinSpaced(order - 1, 1. / order, (order - 1.) / order)
539 .transpose()
540 .replicate(1, order - 1);
541 result.row(1) =
542 Eigen::kroneckerProduct(Eigen::VectorXd::LinSpaced(order - 1, 1. / order,
543 (order - 1.) / order),
544 Eigen::VectorXd::Constant(order - 1, 1.))
545 .transpose();
546 return result;
547}
548
549} // namespace
550
551void WriteToFile(const VtkFile& vtk_file, const std::string& filename) {
552 std::ofstream file(filename, std::ios_base::out | std::ios_base::binary |
553 std::ios_base::trunc);
554 ValidateVtkFile(vtk_file);
555
556 if (!file.is_open()) {
557 throw base::LfException("Could not open file " + filename +
558 " for writing.");
559 }
560 karma::ostream_iterator<char> outit(file);
561
562 bool result = false;
563 if (vtk_file.format == VtkFile::Format::BINARY) {
564 const VtkGrammar<decltype(outit), true> grammar{};
565 result = karma::generate(outit, grammar, vtk_file);
566 } else {
567 const VtkGrammar<decltype(outit), false> grammar{};
568 result = karma::generate(outit, grammar, vtk_file);
569 }
570
571 file.close();
572 if (!result) {
573 throw base::LfException("Karma error");
574 }
575}
576
577// NOLINTNEXTLINE(readability-function-cognitive-complexity)
578VtkWriter::VtkWriter(std::shared_ptr<const mesh::Mesh> mesh,
579 std::string filename, dim_t codim, unsigned char order)
580 : mesh_(std::move(mesh)),
581 filename_(std::move(filename)),
582 codim_(codim),
583 order_(order),
584 aux_node_offset_{{{mesh_, 0}, {mesh_, 1}}} {
585 auto dim_mesh = mesh_->DimMesh();
586 auto dim_world = mesh_->DimWorld();
587 LF_ASSERT_MSG(dim_world > 0 && dim_world <= 4,
588 "VtkWriter supports only dim_world = 1,2 or 3");
589 LF_ASSERT_MSG(codim >= 0 && codim < dim_mesh, "codim out of bounds.");
590 LF_ASSERT_MSG(order > 0 && order <= 10, "order must lie in [1,10]");
591
592 // initialize reference coordinates of auxilliary nodes (if order_ > 0)
593 if (order_ > 1) {
594 aux_nodes_[base::RefEl::kSegment().Id()] = AuxNodesSegment(order);
595 aux_nodes_[base::RefEl::kTria().Id()] = AuxNodesTria(order);
596 aux_nodes_[base::RefEl::kQuad().Id()] = AuxNodesQuad(order);
597 }
598
599 // calculate total number of nodes (main + auxilliary nodes)
600 unsigned int numNodes = mesh_->NumEntities(dim_mesh);
601 for (auto ref_el :
603 if (ref_el.Dimension() > dim_mesh - codim_) {
604 continue;
605 }
606 numNodes += NumAuxNodes(ref_el, order_) * mesh_->NumEntities(ref_el);
607 }
608
609 // insert main nodes:
610 vtk_file_.unstructured_grid.points.resize(numNodes);
611 const Eigen::Matrix<double, 0, 1> zero;
612 for (const auto* p : mesh_->Entities(dim_mesh)) {
613 auto index = mesh_->Index(*p);
614 Eigen::Vector3f coord;
615 if (dim_world == 1) {
616 coord(0) = static_cast<float>(p->Geometry()->Global(zero)(0));
617 coord(1) = 0.F;
618 coord(2) = 0.F;
619 } else if (dim_world == 2) {
620 coord.topRows<2>() = p->Geometry()->Global(zero).cast<float>();
621 coord(2) = 0.F;
622 } else {
623 coord = p->Geometry()->Global(zero).cast<float>();
624 }
625
626 vtk_file_.unstructured_grid.points[index] = std::move(coord);
627 }
628
629 // insert auxilliary nodes:
630 if (order > 1) {
631 auto index_offset = mesh_->NumEntities(dim_mesh);
632 for (char cd = static_cast<char>(dim_mesh - 1);
633 cd >= static_cast<char>(codim); --cd) {
634 for (const auto* e : mesh_->Entities(cd)) {
635 auto ref_el = e->RefEl();
636 if (ref_el == base::RefEl::kTria() && order < 3) {
637 continue;
638 }
639
640 Eigen::MatrixXf coords(3, NumAuxNodes(ref_el, order));
641
642 if (dim_world == 1) {
643 coords.row(0) =
644 e->Geometry()->Global(aux_nodes_[ref_el.Id()]).cast<float>();
645 coords.row(1).setZero();
646 coords.row(2).setZero();
647 } else if (dim_world == 2) {
648 coords.topRows(2) =
649 e->Geometry()->Global(aux_nodes_[ref_el.Id()]).cast<float>();
650 coords.row(2).setZero();
651 } else {
652 coords = e->Geometry()->Global(aux_nodes_[ref_el.Id()]).cast<float>();
653 }
654 for (Eigen::Index i = 0; i < coords.cols(); ++i) {
655 vtk_file_.unstructured_grid.points[index_offset + i] = coords.col(i);
656 }
657 aux_node_offset_[cd](*e) = index_offset;
658 index_offset += coords.cols();
659 }
660 }
661 }
662
663 // compute number of main/aux nodes for each element:
664 std::array<unsigned int, 5> num_nodes{};
665 num_nodes[base::RefEl::kSegment().Id()] =
666 2 + NumAuxNodes(base::RefEl::kSegment(), order);
667 num_nodes[base::RefEl::kTria().Id()] =
668 3 + 3 * NumAuxNodes(base::RefEl::kSegment(), order) +
669 NumAuxNodes(base::RefEl::kTria(), order);
670 num_nodes[base::RefEl::kQuad().Id()] =
671 4 + 4 * NumAuxNodes(base::RefEl::kSegment(), order) +
672 NumAuxNodes(base::RefEl::kQuad(), order);
673
674 // insert elements:
675 vtk_file_.unstructured_grid.cells.resize(mesh_->NumEntities(codim));
676 vtk_file_.unstructured_grid.cell_types.resize(mesh_->NumEntities(codim));
677 auto points_per_segment = NumAuxNodes(base::RefEl::kSegment(), order);
678 for (const auto* e : mesh_->Entities(codim)) {
679 auto index = mesh_->Index(*e);
680 auto ref_el = e->RefEl();
681 auto& node_indices = vtk_file_.unstructured_grid.cells[index];
682 node_indices.reserve(num_nodes[ref_el.Id()]);
683
684 // node indices that make up this cell:
685 for (const auto* p : e->SubEntities(dim_mesh - codim)) {
686 node_indices.push_back(mesh_->Index(*p));
687 }
688
689 // node indices of segments of this cell:
690 auto addSegmentNodes = [&](const mesh::Entity& s, bool invert) {
691 auto start_index = aux_node_offset_[dim_mesh - 1](s);
692 if (!invert) {
693 for (unsigned int i = 0; i < points_per_segment; ++i) {
694 node_indices.push_back(start_index + i);
695 }
696 } else {
697 for (int i = static_cast<int>(points_per_segment - 1); i >= 0; --i) {
698 node_indices.push_back(start_index + i);
699 }
700 }
701 };
702 switch (ref_el) {
704 addSegmentNodes(*e, false);
705 break;
706 case base::RefEl::kTria(): {
707 auto iterator = e->SubEntities(1).begin();
708 auto o_iterator = e->RelativeOrientations().begin();
709 addSegmentNodes(**iterator,
710 (*o_iterator) == mesh::Orientation::negative);
711 ++iterator;
712 ++o_iterator;
713 addSegmentNodes(**iterator,
714 (*o_iterator) == mesh::Orientation::negative);
715 ++iterator;
716 ++o_iterator;
717 addSegmentNodes(**iterator,
718 (*o_iterator) == mesh::Orientation::negative);
719 break;
720 }
721 case base::RefEl::kQuad(): {
722 auto iterator = e->SubEntities(1).begin();
723 auto o_iterator = e->RelativeOrientations().begin();
724 addSegmentNodes(**iterator,
725 (*o_iterator) == mesh::Orientation::negative);
726 ++iterator;
727 ++o_iterator;
728 addSegmentNodes(**iterator,
729 (*o_iterator) == mesh::Orientation::negative);
730 ++iterator;
731 ++o_iterator;
732 addSegmentNodes(**iterator,
733 (*o_iterator) == mesh::Orientation::positive);
734 ++iterator;
735 ++o_iterator;
736 addSegmentNodes(**iterator,
737 (*o_iterator) == mesh::Orientation::positive);
738 break;
739 }
740 default:
741 LF_VERIFY_MSG(false, "something is wrong");
742 }
743
744 // indices in the interior of the cell:
745 if (dim_mesh - codim > 1) {
746 auto offset = aux_node_offset_[0](*e);
747 for (unsigned i = 0; i < NumAuxNodes(e->RefEl(), order); ++i) {
748 node_indices.push_back(offset + i);
749 }
750 }
751
752 if (order_ == 1) {
753 if (ref_el == base::RefEl::kSegment()) {
754 vtk_file_.unstructured_grid.cell_types[index] =
756 } else if (ref_el == base::RefEl::kTria()) {
757 vtk_file_.unstructured_grid.cell_types[index] =
759 } else if (ref_el == base::RefEl::kQuad()) {
760 vtk_file_.unstructured_grid.cell_types[index] =
762 }
763 } else {
764 switch (ref_el) {
766 vtk_file_.unstructured_grid.cell_types[index] =
767 VtkFile::CellType::VTK_LAGRANGE_CURVE;
768 break;
769 case base::RefEl::kTria():
770 vtk_file_.unstructured_grid.cell_types[index] =
771 VtkFile::CellType::VTK_LAGRANGE_TRIANGLE;
772 break;
773 case base::RefEl::kQuad():
774 vtk_file_.unstructured_grid.cell_types[index] =
775 VtkFile::CellType::VTK_LAGRANGE_QUADRILATERAL;
776 break;
777 default:
778 LF_VERIFY_MSG(false, "Unsupported RefElType " << ref_el);
779 }
780 }
781 }
782}
783
785 const std::string& name, const mesh::utils::MeshDataSet<unsigned char>& mds,
786 unsigned char undefined_value) {
787 WriteScalarPointData(name, mds, undefined_value);
788}
789
790void VtkWriter::WritePointData(const std::string& name,
792 char undefined_value) {
793 WriteScalarPointData(name, mds, undefined_value);
794}
795
797 const std::string& name, const mesh::utils::MeshDataSet<unsigned int>& mds,
798 unsigned int undefined_value) {
799 WriteScalarPointData(name, mds, undefined_value);
800}
801
802void VtkWriter::WritePointData(const std::string& name,
804 int undefined_value) {
805 WriteScalarPointData(name, mds, undefined_value);
806}
807
808void VtkWriter::WritePointData(const std::string& name,
810 float undefined_value) {
811 WriteScalarPointData(name, mds, undefined_value);
812}
813
814void VtkWriter::WritePointData(const std::string& name,
816 double undefined_value) {
817 WriteScalarPointData(name, mds, undefined_value);
818}
819
821 const std::string& name,
823 const Eigen::Vector2d& undefined_value) {
824 WriteVectorPointData<2, double>(name, mds, undefined_value);
825}
826
828 const std::string& name,
830 const Eigen::Vector2f& undefined_value) {
831 WriteVectorPointData<2, float>(name, mds, undefined_value);
832}
833
835 const std::string& name,
837 const Eigen::Vector3d& undefined_value) {
838 WriteVectorPointData<3, double>(name, mds, undefined_value);
839}
840
842 const std::string& name,
844 const Eigen::Vector3f& undefined_value) {
845 WriteVectorPointData<3, float>(name, mds, undefined_value);
846}
847
849 const std::string& name,
851 const Eigen::VectorXd& undefined_value) {
852 WriteVectorPointData<Eigen::Dynamic, double>(name, mds, undefined_value);
853}
854
856 const std::string& name,
858 const Eigen::VectorXf& undefined_value) {
859 WriteVectorPointData<Eigen::Dynamic, float>(name, mds, undefined_value);
860}
861
863 const std::string& name, const mesh::utils::MeshDataSet<unsigned char>& mds,
864 unsigned char undefined_value) {
865 WriteScalarCellData(name, mds, undefined_value);
866}
867
868void VtkWriter::WriteCellData(const std::string& name,
870 char undefined_value) {
871 WriteScalarCellData(name, mds, undefined_value);
872}
873
874void VtkWriter::WriteCellData(const std::string& name,
876 unsigned int undefined_value) {
877 WriteScalarCellData(name, mds, undefined_value);
878}
879
880void VtkWriter::WriteCellData(const std::string& name,
882 int undefined_value) {
883 WriteScalarCellData(name, mds, undefined_value);
884}
885
886void VtkWriter::WriteCellData(const std::string& name,
888 float undefined_value) {
889 WriteScalarCellData(name, mds, undefined_value);
890}
891
892void VtkWriter::WriteCellData(const std::string& name,
894 double undefined_value) {
895 WriteScalarCellData(name, mds, undefined_value);
896}
897
899 const std::string& name,
901 const Eigen::Vector2d& undefined_value) {
902 WriteVectorCellData<2, double>(name, mds, undefined_value);
903}
904
906 const std::string& name,
908 const Eigen::Vector2f& undefined_value) {
909 WriteVectorCellData<2, float>(name, mds, undefined_value);
910}
911
913 const std::string& name,
915 const Eigen::Vector3d& undefined_value) {
916 WriteVectorCellData<3, double>(name, mds, undefined_value);
917}
918
920 const std::string& name,
922 const Eigen::Vector3f& undefined_value) {
923 WriteVectorCellData<3, float>(name, mds, undefined_value);
924}
925
927 const std::string& name,
929 const Eigen::VectorXd& undefined_value) {
930 WriteVectorCellData<Eigen::Dynamic, double>(name, mds, undefined_value);
931}
932
934 const std::string& name,
936 const Eigen::VectorXf& undefined_value) {
937 WriteVectorCellData<Eigen::Dynamic, float>(name, mds, undefined_value);
938}
939
940void VtkWriter::WriteGlobalData(const std::string& name,
941 std::vector<int> data) {
942 WriteFieldData(name, std::move(data));
943}
944
945void VtkWriter::WriteGlobalData(const std::string& name,
946 std::vector<float> data) {
947 WriteFieldData(name, std::move(data));
948}
949
950void VtkWriter::WriteGlobalData(const std::string& name,
951 std::vector<double> data) {
952 WriteFieldData(name, std::move(data));
953}
954
955template <class T>
956void VtkWriter::WriteScalarPointData(const std::string& name,
958 T undefined_value) {
959 LF_ASSERT_MSG(order_ == 1,
960 "WritePointData accepts MeshDataSets only if order = 1. For "
961 "order > 1 you have to provide MeshFunctions.");
964 data.data.resize(mesh_->NumEntities(mesh_->DimMesh()));
965 data.name = name;
966 for (const auto* p : mesh_->Entities(mesh_->DimMesh())) {
967 if (mds.DefinedOn(*p)) {
968 data.data[mesh_->Index(*p)] = mds(*p);
969 } else {
970 data.data[mesh_->Index(*p)] = undefined_value;
971 }
972 }
973 vtk_file_.point_data.push_back(std::move(data));
974}
975
976template <int ROWS, class T>
978 const std::string& name,
979 const mesh::utils::MeshDataSet<Eigen::Matrix<T, ROWS, 1>>& mds,
980 const Eigen::Matrix<T, ROWS, 1>& undefined_value) {
981 LF_ASSERT_MSG(order_ == 1,
982 "WritePointData accepts MeshDataSets only if order = 1. For "
983 "order > 1 you have to provide MeshFunctions.");
986 data.data.resize(mesh_->NumEntities(mesh_->DimMesh()));
987 data.name = name;
988 Eigen::Matrix<T, 3, 1> undefined_value_padded;
989 PadWithZeros<ROWS, T>(undefined_value_padded, undefined_value);
990
991 for (const auto* p : mesh_->Entities(mesh_->DimMesh())) {
992 if (mds.DefinedOn(*p)) {
993 PadWithZeros<ROWS, T>(data.data[mesh_->Index(*p)], mds(*p));
994 } else {
995 data.data[mesh_->Index(*p)] = undefined_value_padded;
996 }
997 }
998 vtk_file_.point_data.push_back(std::move(data));
999}
1000
1001template <class T>
1002void VtkWriter::WriteScalarCellData(const std::string& name,
1003 const mesh::utils::MeshDataSet<T>& mds,
1004 T undefined_value) {
1007 data.data.resize(mesh_->NumEntities(codim_));
1008 data.name = name;
1009 for (const auto* e : mesh_->Entities(codim_)) {
1010 if (mds.DefinedOn(*e)) {
1011 data.data[mesh_->Index(*e)] = mds(*e);
1012 } else {
1013 data.data[mesh_->Index(*e)] = undefined_value;
1014 }
1015 }
1016 vtk_file_.cell_data.push_back(std::move(data));
1017}
1018
1019template <int ROWS, class T>
1021 const std::string& name,
1022 const mesh::utils::MeshDataSet<Eigen::Matrix<T, ROWS, 1>>& mds,
1023 const Eigen::Matrix<T, ROWS, 1>& undefined_value) {
1026 data.data.resize(mesh_->NumEntities(codim_));
1027 data.name = name;
1028 Eigen::Matrix<T, 3, 1> undefined_value_padded;
1029 PadWithZeros<ROWS, T>(undefined_value_padded, undefined_value);
1030
1031 for (const auto* p : mesh_->Entities(codim_)) {
1032 if (mds.DefinedOn(*p)) {
1033 PadWithZeros<ROWS, T>(data.data[mesh_->Index(*p)], mds(*p));
1034 } else {
1035 data.data[mesh_->Index(*p)] = undefined_value_padded;
1036 }
1037 }
1038 vtk_file_.cell_data.push_back(std::move(data));
1039}
1040
1041template <class T>
1042void VtkWriter::WriteFieldData(const std::string& name, std::vector<T> data) {
1044 VtkFile::FieldDataArray<T> vtk_data{};
1045 vtk_data.name = name;
1046 vtk_data.data = std::move(data);
1047 vtk_file_.field_data.push_back(std::move(vtk_data));
1048}
1049
1050} // namespace lf::io
result_type operator()(T &&arg)
Definition vtk_writer.cc:29
VariantVisitor(LAMBDA lambda)
Definition vtk_writer.cc:33
RESULT_TYPE result_type
Definition vtk_writer.cc:26
A simple, general purpose exception class that is thrown by LehrFEM++ if something is wrong.
static constexpr RefEl kSegment()
Returns the (1-dimensional) reference segment.
Definition ref_el.h:153
constexpr unsigned int Id() const
Return a unique id for this reference element.
Definition ref_el.h:494
static constexpr RefEl kPoint()
Returns the (0-dimensional) reference point.
Definition ref_el.h:144
static constexpr RefEl kTria()
Returns the reference triangle.
Definition ref_el.h:161
static constexpr RefEl kQuad()
Returns the reference quadrilateral.
Definition ref_el.h:169
Represents one set of attribute data (can be attached to points or cells)
Definition vtk_writer.h:89
std::vector< T > data
Definition vtk_writer.h:92
std::vector< Eigen::Matrix< T, 3, 1 > > data
Definition vtk_writer.h:108
Representation of a VTK file (only relevant) features (advanced usage)
Definition vtk_writer.h:31
FieldData field_data
Definition vtk_writer.h:140
std::vector< boost::variant< FieldDataArray< int >, FieldDataArray< float >, FieldDataArray< double > > > FieldData
Definition vtk_writer.h:125
std::vector< boost::variant< ScalarData< char >, ScalarData< unsigned char >, ScalarData< short >, ScalarData< unsigned short >, ScalarData< int >, ScalarData< unsigned int >, ScalarData< long >, ScalarData< unsigned long >, ScalarData< float >, ScalarData< double >, VectorData< float >, VectorData< double > > > Attributes
Definition vtk_writer.h:119
Format
Nested types.
Definition vtk_writer.h:37
Attributes cell_data
Definition vtk_writer.h:146
Format format
The format of the file.
Definition vtk_writer.h:135
Attributes point_data
Definition vtk_writer.h:143
void WriteScalarPointData(const std::string &name, const mesh::utils::MeshDataSet< T > &mds, T undefined_value)
VtkWriter(const VtkWriter &)=delete
base::dim_t dim_t
Definition vtk_writer.h:274
void WriteCellData(const std::string &name, const mesh::utils::MeshDataSet< unsigned char > &mds, unsigned char undefined_value=0)
Add a new unsigned char attribute dataset that attaches data to the cells of the mesh (i....
void WriteScalarCellData(const std::string &name, const mesh::utils::MeshDataSet< T > &mds, T undefined_value)
void WriteFieldData(const std::string &name, std::vector< T > data)
void WriteGlobalData(const std::string &name, std::vector< int > data)
Write global data into the vtk file that is not related to the mesh at all.
void WriteVectorPointData(const std::string &name, const mesh::utils::MeshDataSet< Eigen::Matrix< T, ROWS, 1 > > &mds, const Eigen::Matrix< T, ROWS, 1 > &undefined_value)
void WriteVectorCellData(const std::string &name, const mesh::utils::MeshDataSet< Eigen::Matrix< T, ROWS, 1 > > &mds, const Eigen::Matrix< T, ROWS, 1 > &undefined_value)
unsigned char order_
Definition vtk_writer.h:791
void CheckAttributeSetName(const DATA &data, const std::string &name)
Definition vtk_writer.h:974
std::shared_ptr< const mesh::Mesh > mesh_
Definition vtk_writer.h:787
void WritePointData(const std::string &name, const mesh::utils::MeshDataSet< unsigned char > &mds, unsigned char undefined_value=0)
Add a new unsigned char attribute dataset that attaches data to the points/nodes of the mesh.
Interface that specifies how data is stored with an entity.
virtual bool DefinedOn(const Entity &e) const =0
Does the dataset store information with this entity?
Mesh input (from file) and output (in various formats) facilities.
std::ostream & operator<<(std::ostream &stream, GMshFileV2::ElementType et)
Output the element type onto the console:
void WriteToFile(const VtkFile &vtk_file, const std::string &filename)