LehrFEM++ 1.0.0
A simple Finite Element Library for teaching
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:
487Eigen::MatrixXd AuxNodesTria(unsigned char order) {
488 if (order < 3) {
489 return Eigen::MatrixXd();
490 }
491 Eigen::MatrixXd result(2, NumAuxNodes(base::RefEl::kTria(), order));
492 if (order == 3) {
493 result << 1. / 3., 1. / 3.;
494 } else if (order == 4) {
495 result << 0.25, 0.5, 0.25, 0.25, 0.25, 0.5;
496 } else {
497 // assign nodes:
498 result(0, 0) = 1. / order;
499 result(1, 0) = 1. / order;
500 result(0, 1) = (order - 2.) / order;
501 result(1, 1) = 1. / order;
502 result(0, 2) = 1. / order;
503 result(1, 2) = (order - 2.) / order;
504
505 if (order > 4) {
506 auto segment_points = AuxNodesSegment(order - 3).eval();
507 // assign edges:
508 result.block(0, 3, 2, order - 4) =
509 Eigen::Vector2d::UnitX() * segment_points * (order - 3.) / order +
510 Eigen::Vector2d(1. / order, 1. / order) *
511 Eigen::MatrixXd::Ones(1, order - 4);
512 result.block(0, order - 1, 2, order - 4) =
513 Eigen::Vector2d(-1, 1) * (order - 3.) / order * segment_points +
514 Eigen::Vector2d((order - 2.) / order, 1. / order) *
515 Eigen::MatrixXd::Ones(1, order - 4);
516 result.block(0, 2 * order - 5, 2, order - 4) =
517 -Eigen::Vector2d::UnitY() * (order - 3.) / order * segment_points +
518 Eigen::Vector2d(1. / order, (order - 2.) / order) *
519 Eigen::MatrixXd::Ones(1, order - 4);
520 }
521 if (order > 5) {
522 // assign interior points recursively:
523 auto points = AuxNodesTria(order - 3);
524 result.block(0, 3 * order - 9, 2, points.cols()) =
525 AuxNodesTria(order - 3) * (order - 6.) / order +
526 Eigen::Vector2d(2. / order, 2. / order) *
527 Eigen::MatrixXd::Ones(1, points.cols());
528 }
529 }
530 return result;
531}
532
533// compute aux reference coordinates of lagrange points on quad:
534Eigen::MatrixXd AuxNodesQuad(unsigned char order) {
535 Eigen::MatrixXd result(2, NumAuxNodes(base::RefEl::kQuad(), order));
536 result.row(0) =
537 Eigen::VectorXd::LinSpaced(order - 1, 1. / order, (order - 1.) / order)
538 .transpose()
539 .replicate(1, order - 1);
540 result.row(1) =
541 Eigen::kroneckerProduct(Eigen::VectorXd::LinSpaced(order - 1, 1. / order,
542 (order - 1.) / order),
543 Eigen::VectorXd::Constant(order - 1, 1.))
544 .transpose();
545 return result;
546}
547
548} // namespace
549
550void WriteToFile(const VtkFile& vtk_file, const std::string& filename) {
551 std::ofstream file(filename, std::ios_base::out | std::ios_base::binary |
552 std::ios_base::trunc);
553 ValidateVtkFile(vtk_file);
554
555 if (!file.is_open()) {
556 throw base::LfException("Could not open file " + filename +
557 " for writing.");
558 }
559 karma::ostream_iterator<char> outit(file);
560
561 bool result = false;
562 if (vtk_file.format == VtkFile::Format::BINARY) {
563 VtkGrammar<decltype(outit), true> grammar{};
564 result = karma::generate(outit, grammar, vtk_file);
565 } else {
566 VtkGrammar<decltype(outit), false> grammar{};
567 result = karma::generate(outit, grammar, vtk_file);
568 }
569
570 file.close();
571 if (!result) {
572 throw base::LfException("Karma error");
573 }
574}
575
576VtkWriter::VtkWriter(std::shared_ptr<const mesh::Mesh> mesh,
577 std::string filename, dim_t codim, unsigned char order)
578 : mesh_(std::move(mesh)),
579 filename_(std::move(filename)),
580 codim_(codim),
581 order_(order),
582 aux_node_offset_{{{mesh_, 0}, {mesh_, 1}}} {
583 auto dim_mesh = mesh_->DimMesh();
584 auto dim_world = mesh_->DimWorld();
585 LF_ASSERT_MSG(dim_world > 0 && dim_world <= 4,
586 "VtkWriter supports only dim_world = 1,2 or 3");
587 LF_ASSERT_MSG(codim >= 0 && codim < dim_mesh, "codim out of bounds.");
588 LF_ASSERT_MSG(order > 0 && order <= 10, "order must lie in [1,10]");
589
590 // initialize reference coordinates of auxilliary nodes (if order_ > 0)
591 if (order_ > 1) {
592 aux_nodes_[base::RefEl::kSegment().Id()] = AuxNodesSegment(order);
593 aux_nodes_[base::RefEl::kTria().Id()] = AuxNodesTria(order);
594 aux_nodes_[base::RefEl::kQuad().Id()] = AuxNodesQuad(order);
595 }
596
597 // calculate total number of nodes (main + auxilliary nodes)
598 unsigned int numNodes = mesh_->NumEntities(dim_mesh);
599 for (auto ref_el :
601 if (ref_el.Dimension() > dim_mesh - codim_) {
602 continue;
603 }
604 numNodes += NumAuxNodes(ref_el, order_) * mesh_->NumEntities(ref_el);
605 }
606
607 // insert main nodes:
608 vtk_file_.unstructured_grid.points.resize(numNodes);
609 Eigen::Matrix<double, 0, 1> zero;
610 for (const auto* p : mesh_->Entities(dim_mesh)) {
611 auto index = mesh_->Index(*p);
612 Eigen::Vector3f coord;
613 if (dim_world == 1) {
614 coord(0) = p->Geometry()->Global(zero)(0);
615 coord(1) = 0.F;
616 coord(2) = 0.F;
617 } else if (dim_world == 2) {
618 coord.topRows<2>() = p->Geometry()->Global(zero).cast<float>();
619 coord(2) = 0.F;
620 } else {
621 coord = p->Geometry()->Global(zero).cast<float>();
622 }
623
624 vtk_file_.unstructured_grid.points[index] = std::move(coord);
625 }
626
627 // insert auxilliary nodes:
628 if (order > 1) {
629 auto index_offset = mesh_->NumEntities(dim_mesh);
630 for (char cd = static_cast<char>(dim_mesh - 1);
631 cd >= static_cast<char>(codim); --cd) {
632 for (const auto* e : mesh_->Entities(cd)) {
633 auto ref_el = e->RefEl();
634 if (ref_el == base::RefEl::kTria() && order < 3) {
635 continue;
636 }
637
638 Eigen::MatrixXf coords(3, NumAuxNodes(ref_el, order));
639
640 if (dim_world == 1) {
641 coords.row(0) =
642 e->Geometry()->Global(aux_nodes_[ref_el.Id()]).cast<float>();
643 coords.row(1).setZero();
644 coords.row(2).setZero();
645 } else if (dim_world == 2) {
646 coords.topRows(2) =
647 e->Geometry()->Global(aux_nodes_[ref_el.Id()]).cast<float>();
648 coords.row(2).setZero();
649 } else {
650 coords = e->Geometry()->Global(aux_nodes_[ref_el.Id()]).cast<float>();
651 }
652 for (Eigen::Index i = 0; i < coords.cols(); ++i) {
653 vtk_file_.unstructured_grid.points[index_offset + i] = coords.col(i);
654 }
655 aux_node_offset_[cd](*e) = index_offset;
656 index_offset += coords.cols();
657 }
658 }
659 }
660
661 // compute number of main/aux nodes for each element:
662 std::array<unsigned int, 5> num_nodes{};
663 num_nodes[base::RefEl::kSegment().Id()] =
664 2 + NumAuxNodes(base::RefEl::kSegment(), order);
665 num_nodes[base::RefEl::kTria().Id()] =
666 3 + 3 * NumAuxNodes(base::RefEl::kSegment(), order) +
667 NumAuxNodes(base::RefEl::kTria(), order);
668 num_nodes[base::RefEl::kQuad().Id()] =
669 4 + 4 * NumAuxNodes(base::RefEl::kSegment(), order) +
670 NumAuxNodes(base::RefEl::kQuad(), order);
671
672 // insert elements:
673 vtk_file_.unstructured_grid.cells.resize(mesh_->NumEntities(codim));
674 vtk_file_.unstructured_grid.cell_types.resize(mesh_->NumEntities(codim));
675 auto points_per_segment = NumAuxNodes(base::RefEl::kSegment(), order);
676 for (const auto* e : mesh_->Entities(codim)) {
677 auto index = mesh_->Index(*e);
678 auto ref_el = e->RefEl();
679 auto& node_indices = vtk_file_.unstructured_grid.cells[index];
680 node_indices.reserve(num_nodes[ref_el.Id()]);
681
682 // node indices that make up this cell:
683 for (const auto* p : e->SubEntities(dim_mesh - codim)) {
684 node_indices.push_back(mesh_->Index(*p));
685 }
686
687 // node indices of segments of this cell:
688 auto addSegmentNodes = [&](const mesh::Entity& s, bool invert) {
689 auto start_index = aux_node_offset_[dim_mesh - 1](s);
690 if (!invert) {
691 for (unsigned int i = 0; i < points_per_segment; ++i) {
692 node_indices.push_back(start_index + i);
693 }
694 } else {
695 for (int i = static_cast<int>(points_per_segment - 1); i >= 0; --i) {
696 node_indices.push_back(start_index + i);
697 }
698 }
699 };
700 switch (ref_el) {
702 addSegmentNodes(*e, false);
703 break;
704 case base::RefEl::kTria(): {
705 const auto* iterator = e->SubEntities(1).begin();
706 const auto* o_iterator = e->RelativeOrientations().begin();
707 addSegmentNodes(**iterator,
708 (*o_iterator) == mesh::Orientation::negative);
709 ++iterator;
710 ++o_iterator;
711 addSegmentNodes(**iterator,
712 (*o_iterator) == mesh::Orientation::negative);
713 ++iterator;
714 ++o_iterator;
715 addSegmentNodes(**iterator,
716 (*o_iterator) == mesh::Orientation::negative);
717 break;
718 }
719 case base::RefEl::kQuad(): {
720 const auto* iterator = e->SubEntities(1).begin();
721 const auto* o_iterator = e->RelativeOrientations().begin();
722 addSegmentNodes(**iterator,
723 (*o_iterator) == mesh::Orientation::negative);
724 ++iterator;
725 ++o_iterator;
726 addSegmentNodes(**iterator,
727 (*o_iterator) == mesh::Orientation::negative);
728 ++iterator;
729 ++o_iterator;
730 addSegmentNodes(**iterator,
731 (*o_iterator) == mesh::Orientation::positive);
732 ++iterator;
733 ++o_iterator;
734 addSegmentNodes(**iterator,
735 (*o_iterator) == mesh::Orientation::positive);
736 break;
737 }
738 default:
739 LF_VERIFY_MSG(false, "something is wrong");
740 }
741
742 // indices in the interior of the cell:
743 if (dim_mesh - codim > 1) {
744 auto offset = aux_node_offset_[0](*e);
745 for (unsigned i = 0; i < NumAuxNodes(e->RefEl(), order); ++i) {
746 node_indices.push_back(offset + i);
747 }
748 }
749
750 if (order_ == 1) {
751 if (ref_el == base::RefEl::kSegment()) {
752 vtk_file_.unstructured_grid.cell_types[index] =
754 } else if (ref_el == base::RefEl::kTria()) {
755 vtk_file_.unstructured_grid.cell_types[index] =
757 } else if (ref_el == base::RefEl::kQuad()) {
758 vtk_file_.unstructured_grid.cell_types[index] =
760 }
761 } else {
762 switch (ref_el) {
764 vtk_file_.unstructured_grid.cell_types[index] =
765 VtkFile::CellType::VTK_LAGRANGE_CURVE;
766 break;
767 case base::RefEl::kTria():
768 vtk_file_.unstructured_grid.cell_types[index] =
769 VtkFile::CellType::VTK_LAGRANGE_TRIANGLE;
770 break;
771 case base::RefEl::kQuad():
772 vtk_file_.unstructured_grid.cell_types[index] =
773 VtkFile::CellType::VTK_LAGRANGE_QUADRILATERAL;
774 break;
775 default:
776 LF_VERIFY_MSG(false, "Unsupported RefElType " << ref_el);
777 }
778 }
779 }
780}
781
783 const std::string& name, const mesh::utils::MeshDataSet<unsigned char>& mds,
784 unsigned char undefined_value) {
785 WriteScalarPointData(name, mds, undefined_value);
786}
787
788void VtkWriter::WritePointData(const std::string& name,
790 char undefined_value) {
791 WriteScalarPointData(name, mds, undefined_value);
792}
793
795 const std::string& name, const mesh::utils::MeshDataSet<unsigned int>& mds,
796 unsigned int undefined_value) {
797 WriteScalarPointData(name, mds, undefined_value);
798}
799
800void VtkWriter::WritePointData(const std::string& name,
802 int undefined_value) {
803 WriteScalarPointData(name, mds, undefined_value);
804}
805
806void VtkWriter::WritePointData(const std::string& name,
808 float undefined_value) {
809 WriteScalarPointData(name, mds, undefined_value);
810}
811
812void VtkWriter::WritePointData(const std::string& name,
814 double undefined_value) {
815 WriteScalarPointData(name, mds, undefined_value);
816}
817
819 const std::string& name,
821 const Eigen::Vector2d& undefined_value) {
822 WriteVectorPointData<2, double>(name, mds, undefined_value);
823}
824
826 const std::string& name,
828 const Eigen::Vector2f& undefined_value) {
829 WriteVectorPointData<2, float>(name, mds, undefined_value);
830}
831
833 const std::string& name,
835 const Eigen::Vector3d& undefined_value) {
836 WriteVectorPointData<3, double>(name, mds, undefined_value);
837}
838
840 const std::string& name,
842 const Eigen::Vector3f& undefined_value) {
843 WriteVectorPointData<3, float>(name, mds, undefined_value);
844}
845
847 const std::string& name,
849 const Eigen::VectorXd& undefined_value) {
850 WriteVectorPointData<Eigen::Dynamic, double>(name, mds, undefined_value);
851}
852
854 const std::string& name,
856 const Eigen::VectorXf& undefined_value) {
857 WriteVectorPointData<Eigen::Dynamic, float>(name, mds, undefined_value);
858}
859
861 const std::string& name, const mesh::utils::MeshDataSet<unsigned char>& mds,
862 unsigned char undefined_value) {
863 WriteScalarCellData(name, mds, undefined_value);
864}
865
866void VtkWriter::WriteCellData(const std::string& name,
868 char undefined_value) {
869 WriteScalarCellData(name, mds, undefined_value);
870}
871
872void VtkWriter::WriteCellData(const std::string& name,
874 unsigned int undefined_value) {
875 WriteScalarCellData(name, mds, undefined_value);
876}
877
878void VtkWriter::WriteCellData(const std::string& name,
880 int undefined_value) {
881 WriteScalarCellData(name, mds, undefined_value);
882}
883
884void VtkWriter::WriteCellData(const std::string& name,
886 float undefined_value) {
887 WriteScalarCellData(name, mds, undefined_value);
888}
889
890void VtkWriter::WriteCellData(const std::string& name,
892 double undefined_value) {
893 WriteScalarCellData(name, mds, undefined_value);
894}
895
897 const std::string& name,
899 const Eigen::Vector2d& undefined_value) {
900 WriteVectorCellData<2, double>(name, mds, undefined_value);
901}
902
904 const std::string& name,
906 const Eigen::Vector2f& undefined_value) {
907 WriteVectorCellData<2, float>(name, mds, undefined_value);
908}
909
911 const std::string& name,
913 const Eigen::Vector3d& undefined_value) {
914 WriteVectorCellData<3, double>(name, mds, undefined_value);
915}
916
918 const std::string& name,
920 const Eigen::Vector3f& undefined_value) {
921 WriteVectorCellData<3, float>(name, mds, undefined_value);
922}
923
925 const std::string& name,
927 const Eigen::VectorXd& undefined_value) {
928 WriteVectorCellData<Eigen::Dynamic, double>(name, mds, undefined_value);
929}
930
932 const std::string& name,
934 const Eigen::VectorXf& undefined_value) {
935 WriteVectorCellData<Eigen::Dynamic, float>(name, mds, undefined_value);
936}
937
938void VtkWriter::WriteGlobalData(const std::string& name,
939 std::vector<int> data) {
940 WriteFieldData(name, std::move(data));
941}
942
943void VtkWriter::WriteGlobalData(const std::string& name,
944 std::vector<float> data) {
945 WriteFieldData(name, std::move(data));
946}
947
948void VtkWriter::WriteGlobalData(const std::string& name,
949 std::vector<double> data) {
950 WriteFieldData(name, std::move(data));
951}
952
953template <class T>
954void VtkWriter::WriteScalarPointData(const std::string& name,
956 T undefined_value) {
957 LF_ASSERT_MSG(order_ == 1,
958 "WritePointData accepts MeshDataSets only if order = 1. For "
959 "order > 1 you have to provide MeshFunctions.");
962 data.data.resize(mesh_->NumEntities(mesh_->DimMesh()));
963 data.name = name;
964 for (const auto* p : mesh_->Entities(mesh_->DimMesh())) {
965 if (mds.DefinedOn(*p)) {
966 data.data[mesh_->Index(*p)] = mds(*p);
967 } else {
968 data.data[mesh_->Index(*p)] = undefined_value;
969 }
970 }
971 vtk_file_.point_data.push_back(std::move(data));
972}
973
974template <int ROWS, class T>
976 const std::string& name,
977 const mesh::utils::MeshDataSet<Eigen::Matrix<T, ROWS, 1>>& mds,
978 const Eigen::Matrix<T, ROWS, 1>& undefined_value) {
979 LF_ASSERT_MSG(order_ == 1,
980 "WritePointData accepts MeshDataSets only if order = 1. For "
981 "order > 1 you have to provide MeshFunctions.");
984 data.data.resize(mesh_->NumEntities(mesh_->DimMesh()));
985 data.name = name;
986 Eigen::Matrix<T, 3, 1> undefined_value_padded;
987 PadWithZeros<ROWS, T>(undefined_value_padded, undefined_value);
988
989 for (const auto* p : mesh_->Entities(mesh_->DimMesh())) {
990 if (mds.DefinedOn(*p)) {
991 PadWithZeros<ROWS, T>(data.data[mesh_->Index(*p)], mds(*p));
992 } else {
993 data.data[mesh_->Index(*p)] = undefined_value_padded;
994 }
995 }
996 vtk_file_.point_data.push_back(std::move(data));
997}
998
999template <class T>
1000void VtkWriter::WriteScalarCellData(const std::string& name,
1001 const mesh::utils::MeshDataSet<T>& mds,
1002 T undefined_value) {
1005 data.data.resize(mesh_->NumEntities(codim_));
1006 data.name = name;
1007 for (const auto* e : mesh_->Entities(codim_)) {
1008 if (mds.DefinedOn(*e)) {
1009 data.data[mesh_->Index(*e)] = mds(*e);
1010 } else {
1011 data.data[mesh_->Index(*e)] = undefined_value;
1012 }
1013 }
1014 vtk_file_.cell_data.push_back(std::move(data));
1015}
1016
1017template <int ROWS, class T>
1019 const std::string& name,
1020 const mesh::utils::MeshDataSet<Eigen::Matrix<T, ROWS, 1>>& mds,
1021 const Eigen::Matrix<T, ROWS, 1>& undefined_value) {
1024 data.data.resize(mesh_->NumEntities(codim_));
1025 data.name = name;
1026 Eigen::Matrix<T, 3, 1> undefined_value_padded;
1027 PadWithZeros<ROWS, T>(undefined_value_padded, undefined_value);
1028
1029 for (const auto* p : mesh_->Entities(codim_)) {
1030 if (mds.DefinedOn(*p)) {
1031 PadWithZeros<ROWS, T>(data.data[mesh_->Index(*p)], mds(*p));
1032 } else {
1033 data.data[mesh_->Index(*p)] = undefined_value_padded;
1034 }
1035 }
1036 vtk_file_.cell_data.push_back(std::move(data));
1037}
1038
1039template <class T>
1040void VtkWriter::WriteFieldData(const std::string& name, std::vector<T> data) {
1042 VtkFile::FieldDataArray<T> vtk_data{};
1043 vtk_data.name = name;
1044 vtk_data.data = std::move(data);
1045 vtk_file_.field_data.push_back(std::move(vtk_data));
1046}
1047
1048} // namespace lf::io
result_type operator()(T &&arg)
Definition: vtk_writer.cc:29
VariantVisitor(LAMBDA lambda)
Definition: vtk_writer.cc:33
LAMBDA lambda_
Definition: vtk_writer.cc:36
RESULT_TYPE result_type
Definition: vtk_writer.cc:26
A simple, general purpose exception class that is thrown by LehrFEM++ if something is wrong.
Definition: lf_exception.h:21
static constexpr RefEl kSegment()
Returns the (1-dimensional) reference segment.
Definition: ref_el.h:150
constexpr unsigned int Id() const
Return a unique id for this reference element.
Definition: ref_el.h:491
static constexpr RefEl kPoint()
Returns the (0-dimensional) reference point.
Definition: ref_el.h:141
static constexpr RefEl kTria()
Returns the reference triangle.
Definition: ref_el.h:158
static constexpr RefEl kQuad()
Returns the reference quadrilateral.
Definition: ref_el.h:166
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:127
Format
Nested types.
Definition: vtk_writer.h:37
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:123
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)
Definition: vtk_writer.cc:954
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....
Definition: vtk_writer.cc:860
void WriteScalarCellData(const std::string &name, const mesh::utils::MeshDataSet< T > &mds, T undefined_value)
Definition: vtk_writer.cc:1000
void WriteFieldData(const std::string &name, std::vector< T > data)
Definition: vtk_writer.cc:1040
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.
Definition: vtk_writer.cc:938
void WriteVectorPointData(const std::string &name, const mesh::utils::MeshDataSet< Eigen::Matrix< T, ROWS, 1 > > &mds, const Eigen::Matrix< T, ROWS, 1 > &undefined_value)
Definition: vtk_writer.cc:975
void WriteVectorCellData(const std::string &name, const mesh::utils::MeshDataSet< Eigen::Matrix< T, ROWS, 1 > > &mds, const Eigen::Matrix< T, ROWS, 1 > &undefined_value)
Definition: vtk_writer.cc:1018
unsigned char order_
Definition: vtk_writer.h:795
void CheckAttributeSetName(const DATA &data, const std::string &name)
Definition: vtk_writer.h:977
std::shared_ptr< const mesh::Mesh > mesh_
Definition: vtk_writer.h:791
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.
Definition: vtk_writer.cc:782
Interface that specifies how data is stored with an entity.
Definition: mesh_data_set.h:36
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.
Definition: gmsh_file_v2.cc:35
std::ostream & operator<<(std::ostream &stream, GMshFileV2::ElementType et)
Output the element type onto the console:
Definition: gmsh_file_v2.cc:38
void WriteToFile(const VtkFile &vtk_file, const std::string &filename)
Definition: vtk_writer.cc:550
span_constexpr std::size_t size(span< T, Extent > const &spn)
Definition: span.h:1463