FEAT 3
Finite Element Analysis Toolbox
Loading...
Searching...
No Matches
ext_vtk_writer.hpp
1// FEAT3: Finite Element Analysis Toolbox, Version 3
2// Copyright (C) 2010 by Stefan Turek & the FEAT group
3// FEAT3 is released under the GNU General Public License version 3,
4// see the file 'copyright.txt' in the top level directory for details.
5
6#pragma once
7
8// includes, FEAT
9#include <kernel/geometry/export_vtk.hpp>
10#include <kernel/geometry/reference_cell_factory.hpp>
11#include <kernel/space/eval_data.hpp>
12#include <kernel/util/math.hpp>
13
14// includes, system
15#include <iostream>
16#include <fstream>
17
18namespace FEAT
19{
20 namespace Space
21 {
23 namespace Intern
24 {
25 template<bool _enable>
26 struct ExtVtkValueWriter
27 {
28 template<typename ExtVtkWriter_>
29 static void write(const ExtVtkWriter_&, std::ostream&) {}
30 };
31
32 template<bool _enable>
33 struct ExtVtkGradWriter
34 {
35 template<typename ExtVtkWriter_>
36 static void write(const ExtVtkWriter_&, std::ostream&) {}
37 };
38
39 template<bool _enable>
40 struct ExtVtkHessWriter
41 {
42 template<typename ExtVtkWriter_>
43 static void write(const ExtVtkWriter_&, std::ostream&) {}
44 };
45
46 template<>
47 struct ExtVtkValueWriter<true>
48 {
49 template<typename ExtVtkWriter_>
50 static void write(const ExtVtkWriter_& writer, std::ostream& os)
51 {
52 writer.write_values(os);
53 }
54 };
55
56 template<>
57 struct ExtVtkGradWriter<true>
58 {
59 template<typename ExtVtkWriter_>
60 static void write(const ExtVtkWriter_& writer, std::ostream& os)
61 {
62 writer.write_gradients(os);
63 }
64 };
65
66 template<>
67 struct ExtVtkHessWriter<true>
68 {
69 template<typename ExtVtkWriter_>
70 static void write(const ExtVtkWriter_& writer, std::ostream& os)
71 {
72 writer.write_hessians(os);
73 }
74 };
75 } // namespace Intern
77
86 template<typename Trafo_>
87 class ExtVtkWriter
88 {
89 friend struct Intern::ExtVtkValueWriter<true>;
90 friend struct Intern::ExtVtkGradWriter<true>;
91 friend struct Intern::ExtVtkHessWriter<true>;
92
93 public:
94 typedef Trafo_ TrafoType;
95 typedef typename TrafoType::MeshType MeshType;
96 typedef typename MeshType::ShapeType ShapeType;
97 typedef typename Geometry::ConformalMesh<ShapeType> RefMeshType;
98
99 protected:
100 typedef typename TrafoType::template Evaluator<>::Type TrafoEval;
101
102 protected:
103 const TrafoType& _trafo;
104 std::ofstream _ofs;
105 RefMeshType* _ref_mesh;
106
107 public:
108 explicit ExtVtkWriter(const TrafoType& trafo, Index num_refines = 1) :
109 _trafo(trafo),
110 _ref_mesh(nullptr)
111 {
112 // create refined reference cell mesh
113 typedef Geometry::ReferenceCellFactory<ShapeType> RefCellFactory;
114 typedef Geometry::StandardRefinery<RefMeshType> RefCellRefinery;
115 RefCellFactory ref_factory;
116 RefMeshType* ref_mesh = new RefMeshType(ref_factory);
117 for(Index i(0); i < num_refines; ++i)
118 {
119 RefMeshType* old_mesh = ref_mesh;
120 {
121 RefCellRefinery refinery(*old_mesh);
122 ref_mesh = new RefMeshType(refinery);
123 }
124 delete old_mesh;
125 }
126 _ref_mesh = ref_mesh;
127 }
128
129 virtual ~ExtVtkWriter()
130 {
131 if(_ref_mesh != nullptr)
132 {
133 delete _ref_mesh;
134 }
135 }
136
140 bool open(String filename)
141 {
142 _ofs.open(filename.c_str());
143 if(!_ofs.is_open() || !_ofs.good())
144 return false;
145
146 // write VTK header
147 _ofs << "# vtk DataFile Version 2.0\n";
148 _ofs << "Generated by FEAT v" << version_major << "." << version_minor << "." << version_patch << "\n";
149 _ofs << "ASCII\n";
150
151 // write mesh type
152 _ofs << "DATASET UNSTRUCTURED_GRID\n";
153
154 // write vertices
155 write_vertices();
156
157 // write indices
158 write_indices();
159
160 // write point data header
161 _ofs << "POINT_DATA " << (_trafo.get_mesh().get_num_entities(ShapeType::dimension)
162 * _ref_mesh->get_num_entities(0)) << "\n";
163
164 // okay
165 return true;
166 }
167
171 void close()
172 {
173 _ofs.close();
174 }
175
188 template<typename Space_, typename VectorType_>
189 void write_values_blocked(String name, const Space_& space, const VectorType_& v)
190 {
191 typedef Space_ SpaceType;
192 typedef typename SpaceType::template Evaluator<TrafoEval>::Type SpaceEval;
193 //static_assert(SpaceEval::can_value != 0, "space cannot evaluate basis function values!");
194 typedef typename SpaceType::DofMappingType DofMappingType;
195 typedef typename SpaceEval::template ConfigTraits<SpaceTags::value> SpaceConfigTraits;
196 typedef typename SpaceConfigTraits::EvalDataType SpaceData;
197 static constexpr TrafoTags trafo_value_config = SpaceConfigTraits::trafo_config;
198 typedef typename TrafoEval::template ConfigTraits<trafo_value_config>::EvalDataType TrafoData;
199 typedef typename VectorType_::ValueType ValueType;
200 Tiny::Vector<ValueType, SpaceEval::max_local_dofs> loc_vec;
201
202 DofMappingType dof_map(space);
203 TrafoEval trafo_eval(_trafo);
204 SpaceEval space_eval(space);
205 TrafoData trafo_data;
206 SpaceData space_data;
207
208 // write basic info
209 _ofs << "VECTORS " << name << " double\n";
210
211 // loop over all cells
212 for(Index cell(trafo_eval.begin()); cell != trafo_eval.end(); ++cell)
213 {
214 // gather local vector
215 dof_map.prepare(cell);
216 for(int i(0); i < dof_map.get_num_local_dofs(); ++i)
217 {
218 loc_vec[i] = v(dof_map.get_index(i));
219 }
220 dof_map.finish();
221
222 // prepare evaluators
223 trafo_eval.prepare(cell);
224 space_eval.prepare(trafo_eval);
225
226 // loop over all points
227 for(Index pt(0); pt < _ref_mesh->get_num_entities(0); ++pt)
228 {
229 // get domain point
230 typename TrafoEval::DomainPointType dom_point;
231 get_dom_point(dom_point, pt);
232
233 // evaluate trafo and space
234 trafo_eval(trafo_data, dom_point);
235 space_eval(space_data, trafo_data);
236
237 // compute value
238 ValueType value(0);
239 for(int i(0); i < space_eval.get_num_local_dofs(); ++i)
240 value += space_data.phi[i].value * loc_vec[i];
241
242 // write value
243 for(int d(0); d < ValueType::n; ++d)
244 _ofs << " " << value(d);
245
246 for(int d(ValueType::n); d < 3; ++d)
247 _ofs << " 0.0";
248 _ofs << "\n";
249 }
250
251 // finish evaluators
252 space_eval.finish();
253 trafo_eval.finish();
254 }
255 }
256
269 template<typename Space_, typename T_>
270 void write_values(String name, const Space_& space, const T_* data)
271 {
272 XASSERT(data != nullptr);
273 typedef Space_ SpaceType;
274 typedef typename SpaceType::template Evaluator<TrafoEval>::Type SpaceEval;
275 static_assert(*(SpaceEval::eval_caps & SpaceTags::value), "space cannot evaluate basis function values!");
276 typedef typename SpaceType::DofMappingType DofMappingType;
277 typedef typename SpaceEval::template ConfigTraits<SpaceTags::value> SpaceConfigTraits;
278 typedef typename SpaceConfigTraits::EvalDataType SpaceData;
279 static constexpr TrafoTags trafo_value_config = SpaceConfigTraits::trafo_config;
280 typedef typename TrafoEval::template ConfigTraits<trafo_value_config>::EvalDataType TrafoData;
281 typedef typename SpaceEval::DataType DataType;
282 Tiny::Vector<DataType, SpaceEval::max_local_dofs> loc_vec;
283
284 DofMappingType dof_map(space);
285 TrafoEval trafo_eval(_trafo);
286 SpaceEval space_eval(space);
287 TrafoData trafo_data;
288 SpaceData space_data;
289
290 // write basic info
291 _ofs << "SCALARS " << name << " double 1\n";
292 _ofs << "LOOKUP_TABLE default\n";
293
294 // loop over all cells
295 for(Index cell(trafo_eval.begin()); cell != trafo_eval.end(); ++cell)
296 {
297 // gather local vector
298 dof_map.prepare(cell);
299 for(int i(0); i < dof_map.get_num_local_dofs(); ++i)
300 {
301 loc_vec[i] = DataType(data[dof_map.get_index(i)]);
302 }
303 dof_map.finish();
304
305 // prepare evaluators
306 trafo_eval.prepare(cell);
307 space_eval.prepare(trafo_eval);
308
309 // loop over all points
310 for(Index pt(0); pt < _ref_mesh->get_num_entities(0); ++pt)
311 {
312 // get domain point
313 typename TrafoEval::DomainPointType dom_point;
314 get_dom_point(dom_point, pt);
315
316 // evaluate trafo and space
317 trafo_eval(trafo_data, dom_point);
318 space_eval(space_data, trafo_data);
319
320 // compute value
321 DataType value(0);
322 for(int i(0); i < space_eval.get_num_local_dofs(); ++i)
323 value += loc_vec[i] * space_data.phi[i].value;
324
325 // write value
326 _ofs << value << "\n";
327 }
328
329 // finish evaluators
330 space_eval.finish();
331 trafo_eval.finish();
332 }
333 }
334
347 template<typename Space_, typename T_>
348 void write_gradients(String name, const Space_& space, const T_* data)
349 {
350 XASSERT(data != nullptr);
351 typedef Space_ SpaceType;
352 typedef typename SpaceType::template Evaluator<TrafoEval>::Type SpaceEval;
353 static_assert(*(SpaceEval::eval_caps & SpaceTags::grad), "space cannot evaluate basis function gradients!");
354 typedef typename SpaceType::DofMappingType DofMappingType;
355 typedef typename SpaceEval::template ConfigTraits<SpaceTags::grad> SpaceConfigTraits;
356 typedef typename SpaceConfigTraits::EvalDataType SpaceData;
357 static constexpr TrafoTags trafo_grad_config = SpaceConfigTraits::trafo_config;
358 typedef typename TrafoEval::template ConfigTraits<trafo_grad_config>::EvalDataType TrafoData;
359 typedef typename SpaceEval::DataType DataType;
360 typedef typename SpaceEval::SpaceEvalTraits SpaceEvalTraits;
361 typedef typename SpaceEvalTraits::BasisGradientType BasisGradientType;
362 Tiny::Vector<DataType, SpaceEval::max_local_dofs> loc_vec;
363
364 DofMappingType dof_map(space);
365 TrafoEval trafo_eval(_trafo);
366 SpaceEval space_eval(space);
367 TrafoData trafo_data;
368 SpaceData space_data;
369
370 // write basic info
371 _ofs << "VECTORS " << name << " double\n";
372
373 // loop over all cells
374 for(Index cell(trafo_eval.begin()); cell != trafo_eval.end(); ++cell)
375 {
376 // gather local vector
377 dof_map.prepare(cell);
378 for(int i(0); i < dof_map.get_num_local_dofs(); ++i)
379 {
380 loc_vec[i] = DataType(data[dof_map.get_index(i)]);
381 }
382 dof_map.finish();
383
384 // prepare evaluators
385 trafo_eval.prepare(cell);
386 space_eval.prepare(trafo_eval);
387
388 // loop over all points
389 for(Index pt(0); pt < _ref_mesh->get_num_entities(0); ++pt)
390 {
391 // get domain point
392 typename TrafoEval::DomainPointType dom_point;
393 get_dom_point(dom_point, pt);
394
395 // evaluate trafo and space
396 trafo_eval(trafo_data, dom_point);
397 space_eval(space_data, trafo_data);
398
399 // compute value
400 BasisGradientType grad;
401 grad.format();
402 for(int i(0); i < space_eval.get_num_local_dofs(); ++i)
403 grad += loc_vec[i] * space_data.phi[i].grad;
404
405 // write value
406 _ofs << grad[0];
407 for(int i(1); i < ShapeType::dimension; ++i)
408 _ofs << " " << grad[i];
409 for(int i(ShapeType::dimension); i < 3; ++i)
410 _ofs << " 0.0";
411 _ofs << "\n";
412 }
413
414 // finish evaluators
415 space_eval.finish();
416 trafo_eval.finish();
417 }
418
419 }
420
433 template<typename Space_, typename T_>
434 void write_hessians(String name, const Space_& space, const T_* data)
435 {
436 XASSERT(data != nullptr);
437 typedef Space_ SpaceType;
438 typedef typename SpaceType::template Evaluator<TrafoEval>::Type SpaceEval;
439 static_assert(*(SpaceEval::eval_caps & SpaceTags::hess), "space cannot evaluate basis function hessians!");
440 typedef typename SpaceType::DofMappingType DofMappingType;
441 typedef typename SpaceEval::template ConfigTraits<SpaceTags::hess> SpaceConfigTraits;
442 typedef typename SpaceConfigTraits::EvalDataType SpaceData;
443 static constexpr TrafoTags trafo_hess_config = SpaceConfigTraits::trafo_config;
444 typedef typename TrafoEval::template ConfigTraits<trafo_hess_config>::EvalDataType TrafoData;
445 typedef typename SpaceEval::DataType DataType;
446 typedef typename SpaceEval::SpaceEvalTraits SpaceEvalTraits;
447 typedef typename SpaceEvalTraits::BasisHessianType BasisHessianType;
448 Tiny::Vector<DataType, SpaceEval::max_local_dofs> loc_vec;
449
450 DofMappingType dof_map(space);
451 TrafoEval trafo_eval(_trafo);
452 SpaceEval space_eval(space);
453 TrafoData trafo_data;
454 SpaceData space_data;
455
456 // write basic info
457 _ofs << "TENSORS " << name << " double\n";
458
459 // loop over all cells
460 for(Index cell(trafo_eval.begin()); cell != trafo_eval.end(); ++cell)
461 {
462 // gather local vector
463 dof_map.prepare(cell);
464 for(int i(0); i < dof_map.get_num_local_dofs(); ++i)
465 {
466 loc_vec[i] = DataType(data[dof_map.get_index(i)]);
467 }
468 dof_map.finish();
469
470 // prepare evaluators
471 trafo_eval.prepare(cell);
472 space_eval.prepare(trafo_eval);
473
474 // loop over all points
475 for(Index pt(0); pt < _ref_mesh->get_num_entities(0); ++pt)
476 {
477 // get domain point
478 typename TrafoEval::DomainPointType dom_point;
479 get_dom_point(dom_point, pt);
480
481 // evaluate trafo and space
482 trafo_eval(trafo_data, dom_point);
483 space_eval(space_data, trafo_data);
484
485 // compute value
486 BasisHessianType hess;
487 hess.format();
488 for(int i(0); i < space_eval.get_num_local_dofs(); ++i)
489 hess += loc_vec[i] * space_data.phi[i].hess;
490
491 // write value
492 for(int i(0); i < ShapeType::dimension; ++i)
493 {
494 _ofs << hess(i,0);
495 for(int j(1); j < ShapeType::dimension; ++j)
496 _ofs << " " << hess(i,j);
497 for(int j(ShapeType::dimension); j < 3; ++j)
498 _ofs << " 0.0";
499 _ofs << "\n";
500 }
501 for(int i(ShapeType::dimension); i < 3; ++i)
502 {
503 _ofs << "0.0";
504 for(int j(1); j < 3; ++j)
505 _ofs << " 0.0";
506 _ofs << "\n";
507 }
508 }
509
510 // finish evaluators
511 space_eval.finish();
512 trafo_eval.finish();
513 }
514 }
515
516 protected:
517 void write_vertices()
518 {
519 // fetch index set and vertex set
520 const TrafoType& trafo = _trafo;
521 const MeshType& mesh = trafo.get_mesh();
522
523 // write vertex set
524 const Index num_cells = mesh.get_num_entities(ShapeType::dimension);
525 const Index num_ref_verts = num_cells * _ref_mesh->get_num_entities(0);
526 _ofs << "POINTS " << num_ref_verts << " double\n";
527
528 // create trafo evaluator and data
529 typedef typename TrafoEval::DomainPointType DomainPointType;
530 typedef typename TrafoEval::template ConfigTraits<TrafoTags::img_point>::EvalDataType TrafoData;
531 DomainPointType dom_point;
532 TrafoEval trafo_eval(trafo);
533 TrafoData trafo_data;
534
535 // loop over all cells
536 for(Index cell(trafo_eval.begin()); cell != trafo_eval.end(); ++cell)
537 {
538 // prepare trafo
539 trafo_eval.prepare(cell);
540
541 // loop over all reference vertices
542 for(Index i(0); i < _ref_mesh->get_num_entities(0); ++i)
543 {
544 // get reference vertex
545 get_dom_point(dom_point, i);
546
547 // apply trafo
548 trafo_eval(trafo_data, dom_point);
549
550 // write image coords
551 _ofs << trafo_data.img_point[0];
552 for(int j(1); j < ShapeType::dimension; ++j)
553 _ofs << " " << trafo_data.img_point[j];
554 for(int j(ShapeType::dimension); j < 3; ++j)
555 _ofs << " 0.0";
556 _ofs << "\n";
557 }
558
559 // finish trafo
560 trafo_eval.finish();
561 }
562 }
563
564 void write_indices()
565 {
566 // fetch index set and vertex set
567 const TrafoType& trafo = _trafo;
568 const MeshType& mesh = trafo.get_mesh();
569
570 // write index set
571 const Index num_cells = mesh.get_num_entities(ShapeType::dimension);
572 const Index num_ref_verts = _ref_mesh->get_num_entities(0);
573 const Index num_ref_cells = _ref_mesh->get_num_entities(ShapeType::dimension);
574 const typename RefMeshType::template IndexSet<ShapeType::dimension,0>::Type& idx =
575 _ref_mesh->template get_index_set<ShapeType::dimension, 0>();
576 int num_idx = idx.get_num_indices();
577
578 typedef Geometry::Intern::VTKShape<ShapeType> VTKHelperType;
579
580 _ofs << "CELLS " << (num_cells*num_ref_cells) << " " << (Index(num_idx+1)*num_cells*num_ref_cells) << "\n";
581
582 // loop over all mesh cells
583 for(Index cell(0); cell < num_cells; ++cell)
584 {
585 // compute vertex offset
586 Index voff = num_ref_verts * cell;
587 for(Index i(0); i < num_ref_cells; ++i)
588 {
589 _ofs << num_idx;
590 for(int j(0); j < num_idx; ++j)
591 {
592 _ofs << " " << voff + idx(i,VTKHelperType::map(j));
593 }
594 _ofs << "\n";
595 }
596 }
597
598 // write cell types
599 _ofs << "CELL_TYPES " << (num_cells*num_ref_cells) << "\n";
600
601 // loop over all mesh cells
602 for(Index cell(0); cell < num_cells*num_ref_cells; ++cell)
603 {
604 _ofs << VTKHelperType::type << "\n";
605 }
606 }
607
608 void get_dom_point(typename TrafoEval::DomainPointType& dom_point, Index i) const
609 {
610 // get reference mesh vertex set
611 const typename RefMeshType::VertexSetType& ref_vtx = _ref_mesh->get_vertex_set();
612
613 // set domain point
614 for(int j(0); j < ShapeType::dimension; ++j)
615 dom_point[j] = ref_vtx[i][j];
616 }
617 };
618 } // namespace Space
619} // namespace FEAT
XASSERT
#define XASSERT(expr)
Assertion macro definition.
Definition: assertion.hpp:262
FEAT::Geometry::ConformalMesh
Conformal mesh class template.
Definition: conformal_mesh.hpp:43
FEAT::Geometry::ConformalMesh::get_num_entities
Index get_num_entities(int dim) const
Returns the number of entities.
Definition: conformal_mesh.hpp:259
FEAT::Geometry::ConformalMesh::VertexSetType
VertexSet< num_coords_, Coord_ > VertexSetType
Vertex set type.
Definition: conformal_mesh.hpp:54
FEAT::Geometry::ConformalMesh::get_vertex_set
VertexSetType & get_vertex_set()
Returns a reference to the vertex set of the mesh.
Definition: conformal_mesh.hpp:507
FEAT::Geometry::StandardRefinery
Standard Refinery class template.
Definition: factory.hpp:58
FEAT::Space::ExtVtkWriter
Extended VTK writer class template.
Definition: ext_vtk_writer.hpp:88
FEAT::Space::ExtVtkWriter::close
void close()
Closes the VTK file.
Definition: ext_vtk_writer.hpp:171
FEAT::Space::ExtVtkWriter::write_gradients
void write_gradients(String name, const Space_ &space, const T_ *data)
Writes a finite element function gradient field to the VTK file.
Definition: ext_vtk_writer.hpp:348
FEAT::Space::ExtVtkWriter::write_values_blocked
void write_values_blocked(String name, const Space_ &space, const VectorType_ &v)
Writes a finite element function to the VTK file.
Definition: ext_vtk_writer.hpp:189
FEAT::Space::ExtVtkWriter::open
bool open(String filename)
Opens a VTK files and writes basic information.
Definition: ext_vtk_writer.hpp:140
FEAT::Space::ExtVtkWriter::write_hessians
void write_hessians(String name, const Space_ &space, const T_ *data)
Writes a finite element function Hessian field to the VTK file.
Definition: ext_vtk_writer.hpp:434
FEAT::Space::ExtVtkWriter::write_values
void write_values(String name, const Space_ &space, const T_ *data)
Writes a finite element function to the VTK file.
Definition: ext_vtk_writer.hpp:270
FEAT::String
String class implementation.
Definition: string.hpp:46
FEAT::Tiny::Vector
Tiny Vector class template.
Definition: tiny_algebra.hpp:52
FEAT
FEAT namespace.
Definition: adjactor.hpp:12
FEAT::version_major
static constexpr int version_major
FEAT major version number.
Definition: base_header.hpp:110
FEAT::version_patch
static constexpr int version_patch
FEAT patch version number.
Definition: base_header.hpp:114
FEAT::version_minor
static constexpr int version_minor
FEAT minor version number.
Definition: base_header.hpp:112
FEAT::SpaceTags::value
@ value
specifies whether the space should supply basis function values
FEAT::SpaceTags::hess
@ hess
specifies whether the space should supply basis function hessians
FEAT::SpaceTags::grad
@ grad
specifies whether the space should supply basis function gradients
FEAT::Index
std::uint64_t Index
Index data type.
Definition: base_header.hpp:122
FEAT::TrafoTags
TrafoTags
Trafo configuration tags enum.
Definition: eval_tags.hpp:22
FEAT::TrafoTags::dom_point
@ dom_point
specifies whether the trafo should supply domain point coordinates