Line data Source code
1 : /*
2 : * Copyright (c) 2010-2015: G-CSC, Goethe University Frankfurt
3 : * Author: Andreas Vogel
4 : *
5 : * This file is part of UG4.
6 : *
7 : * UG4 is free software: you can redistribute it and/or modify it under the
8 : * terms of the GNU Lesser General Public License version 3 (as published by the
9 : * Free Software Foundation) with the following additional attribution
10 : * requirements (according to LGPL/GPL v3 §7):
11 : *
12 : * (1) The following notice must be displayed in the Appropriate Legal Notices
13 : * of covered and combined works: "Based on UG4 (www.ug4.org/license)".
14 : *
15 : * (2) The following notice must be displayed at a prominent place in the
16 : * terminal output of covered works: "Based on UG4 (www.ug4.org/license)".
17 : *
18 : * (3) The following bibliography is recommended for citation and must be
19 : * preserved in all covered files:
20 : * "Reiter, S., Vogel, A., Heppner, I., Rupp, M., and Wittum, G. A massively
21 : * parallel geometric multigrid solver on hierarchically distributed grids.
22 : * Computing and visualization in science 16, 4 (2013), 151-164"
23 : * "Vogel, A., Reiter, S., Rupp, M., Nägel, A., and Wittum, G. UG4 -- a novel
24 : * flexible software system for simulating pde based models on high performance
25 : * computers. Computing and visualization in science 16, 4 (2013), 165-179"
26 : *
27 : * This program is distributed in the hope that it will be useful,
28 : * but WITHOUT ANY WARRANTY; without even the implied warranty of
29 : * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
30 : * GNU Lesser General Public License for more details.
31 : */
32 :
33 : //other libraries
34 : #include <cstdio>
35 : #include <iostream>
36 : #include <cstring>
37 : #include <string>
38 : #include <algorithm>
39 :
40 : // ug4 libraries
41 : #include "common/log.h"
42 : #include "common/util/string_util.h"
43 : #include "common/util/endian_detection.h"
44 : #include "common/profiler/profiler.h"
45 : #include "common/util/provider.h"
46 : #include "lib_disc/common/function_group.h"
47 : #include "lib_disc/common/groups_util.h"
48 : #include "lib_disc/common/multi_index.h"
49 : #include "lib_disc/reference_element/reference_element.h"
50 : #include "lib_disc/local_finite_element/local_finite_element_provider.h"
51 : #include "lib_disc/spatial_disc/disc_util/fv_util.h"
52 : #include "lib_grid/algorithms/debug_util.h"
53 :
54 : #ifdef UG_PARALLEL
55 : #include "pcl/pcl_base.h"
56 : #include "lib_algebra/parallelization/parallel_storage_type.h"
57 : #endif
58 :
59 : // own interface
60 : #include "vtkoutput.h"
61 :
62 : namespace ug{
63 :
64 : /**
65 : * Writes a float in the VTU-compatible ascii format:
66 : */
67 : template <int TDim>
68 0 : VTKFileWriter& VTKOutput<TDim>::
69 : write_asc_float(VTKFileWriter& File, float data)
70 : {
71 0 : if (std::abs (data) < std::numeric_limits<float>::min ()) // a protection against the denormalized floats
72 0 : File << '0';
73 : else
74 0 : File << data;
75 0 : return File;
76 : }
77 :
78 :
79 : /* Functions that write data into the VTU file depending on the type.
80 : * Note that the data may be either scalars, 3-dimensional vectors or 3x3 tensors.
81 : * The binary values should be represented in the Float32 format, the ascii ones must
82 : * be representable in this format.
83 : */
84 :
85 : template <int TDim>
86 0 : void VTKOutput<TDim>::
87 : write_item_to_file(VTKFileWriter& File, float data)
88 : {
89 0 : if(m_bBinary)
90 0 : File << (float) data;
91 : else
92 0 : write_asc_float (File, data) << ' ';
93 0 : }
94 :
95 : // fill position data up with zeros if dim < 3.
96 : template <int TDim>
97 0 : void VTKOutput<TDim>::
98 : write_item_to_file(VTKFileWriter& File, const ug::MathVector<1>& data)
99 : {
100 0 : if(m_bBinary)
101 0 : File << (float) data[0] << (float) 0.f << (float) 0.f;
102 : else
103 0 : write_asc_float (File, data[0]) << " 0 0 ";
104 0 : }
105 :
106 : template <int TDim>
107 0 : void VTKOutput<TDim>::
108 : write_item_to_file(VTKFileWriter& File, const ug::MathVector<2>& data)
109 : {
110 0 : if(m_bBinary)
111 0 : File << (float) data[0] << (float) data[1] << (float) 0.f;
112 : else
113 : {
114 0 : write_asc_float (File, data[0]) << ' ';
115 0 : write_asc_float (File, data[1]) << " 0 ";
116 : }
117 0 : }
118 :
119 : template <int TDim>
120 0 : void VTKOutput<TDim>::
121 : write_item_to_file(VTKFileWriter& File, const ug::MathVector<3>& data)
122 : {
123 0 : if(m_bBinary)
124 0 : File << (float) data[0] << (float) data[1] << (float) data[2];
125 : else
126 : {
127 0 : write_asc_float (File, data[0]) << ' ';
128 0 : write_asc_float (File, data[1]) << ' ';
129 0 : write_asc_float (File, data[2]) << ' ';
130 : }
131 0 : }
132 :
133 : template <int TDim>
134 0 : void VTKOutput<TDim>::
135 : write_item_to_file(VTKFileWriter& File, const ug::MathMatrix<1,1>& data)
136 : {
137 0 : if(m_bBinary)
138 0 : File << (float) data(0,0) << (float) 0.f << (float) 0.f << (float) 0.f << (float) 0.f << (float) 0.f << (float) 0.f << (float) 0.f << (float) 0.f;
139 : else
140 0 : write_asc_float (File, data(0,0)) << " 0 0 0 0 0 0 0 0 ";
141 0 : }
142 :
143 : template <int TDim>
144 0 : void VTKOutput<TDim>::
145 : write_item_to_file(VTKFileWriter& File, const ug::MathMatrix<2,2>& data)
146 : {
147 0 : if(m_bBinary)
148 : {
149 0 : File << (float) data(0,0) << (float) data(0,1) << (float) 0.f;
150 0 : File << (float) data(1,0) << (float) data(1,1) << (float) 0.f << (float) 0.f << (float) 0.f << (float) 0.f;
151 : }
152 : else
153 : {
154 0 : write_asc_float (File, data(0,0)) << ' '; write_asc_float (File, data(0,1)) << " 0 ";
155 0 : write_asc_float (File, data(1,0)) << ' '; write_asc_float (File, data(1,1)) << " 0 0 0 0 ";
156 : }
157 0 : }
158 :
159 : template <int TDim>
160 0 : void VTKOutput<TDim>::
161 : write_item_to_file(VTKFileWriter& File, const ug::MathMatrix<3,3>& data)
162 : {
163 0 : if(m_bBinary)
164 : {
165 0 : File << (float) data(0,0) << (float) data(0,1) << (float) data(0,2);
166 0 : File << (float) data(1,0) << (float) data(1,1) << (float) data(1,2);
167 0 : File << (float) data(2,0) << (float) data(2,1) << (float) data(2,2);
168 : }
169 : else
170 : {
171 0 : write_asc_float (File, data(0,0)) << ' '; write_asc_float (File, data(0,1)) << ' '; write_asc_float (File, data(0,2)) << ' ';
172 0 : write_asc_float (File, data(1,0)) << ' '; write_asc_float (File, data(1,1)) << ' '; write_asc_float (File, data(1,2)) << ' ';
173 0 : write_asc_float (File, data(2,0)) << ' '; write_asc_float (File, data(2,1)) << ' '; write_asc_float (File, data(2,2)) << ' ';
174 : }
175 0 : }
176 :
177 : /**
178 : * VTU output for the data based on the given grid function (u), for all the subsets
179 : * where this function is defined.
180 : */
181 : template <int TDim>
182 : template <typename TFunction>
183 0 : void VTKOutput<TDim>::
184 : print(const char* filename, TFunction& u, int step, number time, bool makeConsistent)
185 : {
186 : PROFILE_FUNC();
187 : #ifdef UG_PARALLEL
188 : if(makeConsistent)
189 : if(!u.change_storage_type(PST_CONSISTENT))
190 : UG_THROW("VTK::print: Cannot change storage type to consistent.");
191 : #endif
192 :
193 : // check functions
194 : bool bEverywhere = true;
195 0 : for(size_t fct = 0; fct < u.num_fct(); ++fct)
196 : {
197 : // check if function is defined everywhere
198 0 : if(!u.is_def_everywhere(fct))
199 : bEverywhere = false;
200 : }
201 :
202 : // in case that all functions are defined everywhere, we write the grid as
203 : // a whole. If not, we must write each subset separately and group the files
204 : // later using a *.pvd file.
205 0 : if(bEverywhere)
206 : {
207 : // write whole grid to a single file
208 : try
209 : {
210 0 : SubsetGroup all_ss (u.domain()->subset_handler ());
211 0 : all_ss.add_all ();
212 0 : print_subsets(filename, u, all_ss, step, time, makeConsistent);
213 : }
214 0 : UG_CATCH_THROW("VTK::print: Failed to write grid.");
215 : }
216 : else
217 : {
218 : // loop subsets
219 0 : for(int si = 0; si < u.num_subsets(); ++si)
220 : {
221 : // write each subset to a single file
222 : try{
223 0 : print_subset(filename, u, si, step, time, makeConsistent);
224 : }
225 0 : UG_CATCH_THROW("VTK::print: Failed to write Subset "<< si << ".");
226 : }
227 :
228 : // write grouping pvd file
229 : try{
230 0 : write_subset_pvd(u.num_subsets(), filename, step, time);
231 : }
232 0 : UG_CATCH_THROW("VTK::print: Failed to write pvd-file.");
233 : }
234 : #ifdef UG_PARALLEL
235 : PCL_DEBUG_BARRIER_ALL();
236 : #endif
237 0 : }
238 :
239 :
240 : /**
241 : * VTU output for the data based on the given grid function (u), for only one given subset
242 : * si. The file name is composed with the number of the subset.
243 : */
244 : template <int TDim>
245 : template <typename TFunction>
246 0 : void VTKOutput<TDim>::
247 : print_subset(const char* filename, TFunction& u, int si, int step, number time, bool makeConsistent)
248 : {
249 : PROFILE_FUNC();
250 :
251 : #ifdef UG_PARALLEL
252 : if(makeConsistent)
253 : if(!u.change_storage_type(PST_CONSISTENT))
254 : UG_THROW("VTK::print_subset: Cannot change storage type to consistent.");
255 : #endif
256 :
257 : // get the grid associated to the solution
258 0 : Grid& grid = *u.domain()->grid();
259 0 : SubsetGroup ssg (u.domain()->subset_handler ());
260 0 : ssg.add (si);
261 :
262 : // attach help indices
263 : typedef ug::Attachment<int> AVrtIndex;
264 : AVrtIndex aVrtIndex;
265 : Grid::VertexAttachmentAccessor<AVrtIndex> aaVrtIndex;
266 : grid.attach_to_vertices(aVrtIndex);
267 : aaVrtIndex.access(grid, aVrtIndex);
268 :
269 : // get rank of process
270 : int rank = 0;
271 : #ifdef UG_PARALLEL
272 : rank = pcl::ProcRank();
273 : #endif
274 :
275 : // get name for *.vtu file
276 : std::string name;
277 : try{
278 0 : vtu_filename(name, filename, rank, si, u.num_subsets()-1, step);
279 : }
280 0 : UG_CATCH_THROW("VTK::print_subset: Failed to write vtu-file.");
281 :
282 :
283 : // open the file
284 : try
285 : {
286 0 : VTKFileWriter File(name.c_str());
287 :
288 : // bool if time point should be written to *.vtu file
289 : // in parallel we must not (!) write it to the *.vtu file, but to the *.pvtu
290 : bool bTimeDep = (step >= 0);
291 : #ifdef UG_PARALLEL
292 : if(pcl::NumProcs() > 1) bTimeDep = false;
293 : #endif
294 :
295 : // header
296 0 : File << VTKFileWriter::normal;
297 0 : File << "<?xml version=\"1.0\"?>\n";
298 :
299 0 : write_comment(File);
300 :
301 0 : File << "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\" byte_order=\"";
302 0 : if(IsLittleEndian()) File << "LittleEndian";
303 : else File << "BigEndian";
304 0 : File << "\">\n";
305 :
306 : // writing time point
307 0 : if(bTimeDep)
308 : {
309 0 : File << " <Time timestep=\""<<time<<"\"/>\n";
310 : }
311 :
312 : // opening the grid
313 0 : File << " <UnstructuredGrid>\n";
314 :
315 : // get dimension of grid-piece
316 0 : int dim = DimensionOfSubset(*u.domain()->subset_handler(), si);
317 :
318 : // write piece of grid
319 0 : if(dim >= 0)
320 : {
321 : try{
322 0 : write_grid_solution_piece(File, aaVrtIndex, grid, u, time, ssg, dim);
323 : }
324 0 : UG_CATCH_THROW("VTK::print_subset: Can not write Subset: "<<si);
325 : }
326 : else
327 : {
328 : // if dim < 0, something is wrong with grid, unless there are no elements in the grid
329 0 : if((si >=0) && u.domain()->subset_handler()->template num<Vertex>(si) != 0)
330 : {
331 0 : UG_THROW("VTK::print_subset: Dimension of grid/subset not"
332 : " detected correctly although grid objects present.");
333 : }
334 :
335 0 : write_empty_grid_piece(File);
336 : }
337 :
338 : // write closing xml tags
339 0 : File << VTKFileWriter::normal;
340 0 : File << " </UnstructuredGrid>\n";
341 0 : File << "</VTKFile>\n";
342 :
343 : // detach help indices
344 : grid.detach_from_vertices(aVrtIndex);
345 :
346 : #ifdef UG_PARALLEL
347 : // write grouping *.pvtu file in parallel case
348 : try{
349 : write_pvtu(u, filename, si, step, time);
350 : }
351 : UG_CATCH_THROW("VTK::print_subset: Failed to write pvtu-file.");
352 : #endif
353 :
354 : // remember time step
355 0 : if(step >= 0)
356 : {
357 : // get vector of time points for the name
358 0 : std::vector<number>& vTimestep = m_mTimestep[filename];
359 :
360 : // resize the vector
361 0 : vTimestep.resize(step+1);
362 :
363 : // write time point
364 0 : vTimestep[step] = time;
365 : }
366 :
367 0 : }
368 0 : UG_CATCH_THROW("VTK::print_subset: Failed to open file: '" << filename << "' for output");
369 : #ifdef UG_PARALLEL
370 : PCL_DEBUG_BARRIER_ALL();
371 : #endif
372 0 : }
373 :
374 :
375 : /**
376 : * VTU output for the data based on the given grid function (u), only for the subsets from
377 : * a given subset group.
378 : */
379 : template <int TDim>
380 : template <typename TFunction>
381 0 : void VTKOutput<TDim>::
382 : print_subsets(const char* filename, TFunction& u, const SubsetGroup& ssGrp, int step, number time, bool makeConsistent)
383 : {
384 : PROFILE_FUNC();
385 :
386 : #ifdef UG_PARALLEL
387 : if(makeConsistent)
388 : if(!u.change_storage_type(PST_CONSISTENT))
389 : UG_THROW("VTK::print_subset: Cannot change storage type to consistent.");
390 : #endif
391 :
392 : // get the grid associated to the solution
393 0 : Grid& grid = *u.domain()->grid();
394 :
395 : // attach help indices
396 : typedef ug::Attachment<int> AVrtIndex;
397 : AVrtIndex aVrtIndex;
398 : Grid::VertexAttachmentAccessor<AVrtIndex> aaVrtIndex;
399 : grid.attach_to_vertices(aVrtIndex);
400 : aaVrtIndex.access(grid, aVrtIndex);
401 :
402 : // get rank of process
403 : int rank = 0;
404 : #ifdef UG_PARALLEL
405 : rank = pcl::ProcRank();
406 : #endif
407 :
408 : // get name for *.vtu file
409 : std::string name;
410 : try{
411 0 : vtu_filename(name, filename, rank, -1, u.num_subsets()-1, step); // "si == -1" because we do not want any subset prefixes!
412 : }
413 0 : UG_CATCH_THROW("VTK::print_subsets: Failed to write vtu-file.");
414 :
415 :
416 : // open the file
417 : try
418 : {
419 0 : VTKFileWriter File(name.c_str());
420 :
421 : // bool if time point should be written to *.vtu file
422 : // in parallel we must not (!) write it to the *.vtu file, but to the *.pvtu
423 : bool bTimeDep = (step >= 0);
424 : #ifdef UG_PARALLEL
425 : if(pcl::NumProcs() > 1) bTimeDep = false;
426 : #endif
427 :
428 : // header
429 0 : File << VTKFileWriter::normal;
430 0 : File << "<?xml version=\"1.0\"?>\n";
431 :
432 0 : write_comment(File);
433 :
434 0 : File << "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\" byte_order=\"";
435 0 : if(IsLittleEndian()) File << "LittleEndian";
436 : else File << "BigEndian";
437 0 : File << "\">\n";
438 :
439 : // writing time point
440 0 : if(bTimeDep)
441 : {
442 0 : File << " <Time timestep=\""<<time<<"\"/>\n";
443 : }
444 :
445 : // opening the grid
446 0 : File << " <UnstructuredGrid>\n";
447 :
448 : // get dimension of grid-piece: the highest dimension of the specified subsets
449 : int dim = -1;
450 0 : for(size_t i = 0; i < ssGrp.size(); i++)
451 : {
452 0 : int ssDim = DimensionOfSubset(*u.domain()->subset_handler(), ssGrp[i]);
453 0 : if(dim < ssDim)
454 : dim = ssDim;
455 : }
456 :
457 : // write piece of grid
458 0 : if(dim >= 0)
459 : {
460 : try{
461 0 : write_grid_solution_piece(File, aaVrtIndex, grid, u, time, ssGrp, dim);
462 : }
463 0 : UG_CATCH_THROW("VTK::print_subsets: Can not write the subsets");
464 : }
465 : else
466 : {
467 0 : UG_THROW("VTK::print_subsets: Dimension of grid/subset not"
468 : " detected correctly although grid objects present.");
469 : }
470 :
471 : // write closing xml tags
472 0 : File << VTKFileWriter::normal;
473 0 : File << " </UnstructuredGrid>\n";
474 0 : File << "</VTKFile>\n";
475 :
476 : // detach help indices
477 : grid.detach_from_vertices(aVrtIndex);
478 :
479 : #ifdef UG_PARALLEL
480 : // write grouping *.pvtu file in parallel case
481 : try{
482 : write_pvtu(u, filename, -1, step, time); // "-1" because we do not want any subset prefixes!
483 : }
484 : UG_CATCH_THROW("VTK::print_subsets: Failed to write pvtu-file.");
485 : #endif
486 :
487 : // remember time step
488 0 : if(step >= 0)
489 : {
490 : // get vector of time points for the name
491 0 : std::vector<number>& vTimestep = m_mTimestep[filename];
492 :
493 : // resize the vector
494 0 : vTimestep.resize(step+1);
495 :
496 : // write time point
497 0 : vTimestep[step] = time;
498 : }
499 :
500 0 : }
501 0 : UG_CATCH_THROW("VTK::print_subset: Failed open file: '" << filename << "' for output");
502 : #ifdef UG_PARALLEL
503 : PCL_DEBUG_BARRIER_ALL();
504 : #endif
505 0 : }
506 :
507 :
508 : /**
509 : * VTU output for the data based on the given grid function (u), only for the subsets from
510 : * a given list of subset names.
511 : */
512 : template <int TDim>
513 : template <typename TFunction>
514 0 : void VTKOutput<TDim>::
515 : print_subsets(const char* filename, TFunction& u, const char* ssNames, int step, number time, bool makeConsistent)
516 : {
517 : std::vector<std::string> v_ss_names;
518 0 : SubsetGroup ssGrp (u.domain()->subset_handler());
519 :
520 0 : TokenizeTrimString (ssNames, v_ss_names);
521 0 : ssGrp.add (v_ss_names);
522 0 : print_subsets(filename, u, ssGrp, step, time, makeConsistent);
523 0 : }
524 :
525 :
526 : ////////////////////////////////////////////////////////////////////////////////
527 : // writing pieces
528 : ////////////////////////////////////////////////////////////////////////////////
529 :
530 :
531 : /**
532 : * This function composes the vtu sections with the coordinates of the grid points and
533 : * the connectivities (cells), a part of the unstructured grid description in the vtu file.
534 : */
535 : template <int TDim>
536 : template <typename T>
537 0 : void VTKOutput<TDim>::
538 : write_points_cells_piece(VTKFileWriter& File,
539 : Grid::VertexAttachmentAccessor<Attachment<int> >& aaVrtIndex,
540 : const Grid::VertexAttachmentAccessor<Attachment<MathVector<TDim> > >& aaPos,
541 : Grid& grid, const T& iterContainer, const SubsetGroup& ssGrp, const int dim,
542 : const int numVert, const int numElem, const int numConn)
543 : {
544 : // write vertices of this piece
545 : try{
546 0 : write_points<T>(File, aaVrtIndex, aaPos, grid, iterContainer, ssGrp, dim, numVert);
547 : }
548 0 : UG_CATCH_THROW("VTK::write_points_cells_piece: Failed to write vertices (Points).");
549 :
550 : // write elements of this piece
551 : try{
552 0 : write_cells(File, aaVrtIndex, grid, iterContainer, ssGrp, dim, numElem, numConn);
553 : }
554 0 : UG_CATCH_THROW("VTK::write_points_cells_piece: Failed to write grid elements (Cells).");
555 0 : }
556 :
557 : /**
558 : * Composes only the vtu-piece of the description of the unstructured grid.
559 : */
560 : template <int TDim>
561 : template <typename T>
562 0 : void VTKOutput<TDim>::
563 : write_grid_piece(VTKFileWriter& File,
564 : Grid::VertexAttachmentAccessor<Attachment<int> >& aaVrtIndex,
565 : const Grid::VertexAttachmentAccessor<Attachment<MathVector<TDim> > >& aaPos,
566 : Grid& grid, const T& iterContainer, SubsetGroup& ssGrp, int dim)
567 : {
568 : // counters
569 0 : int numVert = 0, numElem = 0, numConn = 0;
570 :
571 : // Count needed sizes for vertices, elements and connections
572 : try{
573 0 : count_piece_sizes(grid, iterContainer, ssGrp, dim, numVert, numElem, numConn);
574 : }
575 0 : UG_CATCH_THROW("VTK::write_piece: Failed to count piece sizes.");
576 :
577 : // write the beginning of the piece, indicating the number of vertices
578 : // and the number of elements for this piece of the grid.
579 0 : File << VTKFileWriter::normal;
580 0 : File << " <Piece NumberOfPoints=\""<<numVert<<
581 0 : "\" NumberOfCells=\""<<numElem<<"\">\n";
582 :
583 : // write grid
584 : write_points_cells_piece<T>
585 0 : (File, aaVrtIndex, aaPos, grid, iterContainer, ssGrp, dim, numVert, numElem, numConn);
586 :
587 : // write closing tag
588 0 : File << VTKFileWriter::normal;
589 0 : File << " </Piece>\n";
590 0 : }
591 :
592 : /**
593 : * This function writes a piece of the grid to the vtk file. First the
594 : * geometric data is written (points, cells, connectivity). Then the data
595 : * associated to the points or cells is written.
596 : */
597 : template <int TDim>
598 : template <typename TFunction>
599 0 : void VTKOutput<TDim>::
600 : write_grid_solution_piece(VTKFileWriter& File,
601 : Grid::VertexAttachmentAccessor<Attachment<int> >& aaVrtIndex,
602 : Grid& grid, TFunction& u, number time,
603 : const SubsetGroup& ssGrp, const int dim)
604 : {
605 : // counters
606 0 : int numVert = 0, numElem = 0, numConn = 0;
607 :
608 : // Count needed sizes for vertices, elements and connections
609 : try{
610 0 : count_piece_sizes(grid, u, ssGrp, dim, numVert, numElem, numConn);
611 : }
612 0 : UG_CATCH_THROW("VTK::write_piece: Failed to count piece sizes.");
613 :
614 : // write the beginning of the piece, indicating the number of vertices
615 : // and the number of elements for this piece of the grid.
616 0 : File << VTKFileWriter::normal;
617 0 : File << " <Piece NumberOfPoints=\""<<numVert<<
618 0 : "\" NumberOfCells=\""<<numElem<<"\">\n";
619 :
620 : // write grid
621 : write_points_cells_piece<TFunction>
622 0 : (File, aaVrtIndex, u.domain()->position_accessor(), grid, u, ssGrp, dim, numVert, numElem, numConn);
623 :
624 : // add all components if 'selectAll' chosen
625 0 : if(m_bSelectAll){
626 0 : for(size_t fct = 0; fct < u.num_fct(); ++fct){
627 0 : if(!vtk_name_used(u.name(fct).c_str())){
628 0 : if(LocalFiniteElementProvider::continuous(u.local_finite_element_id(fct))){
629 0 : select_nodal(u.name(fct).c_str(), u.name(fct).c_str());
630 : }else{
631 0 : select_element(u.name(fct).c_str(), u.name(fct).c_str());
632 : }
633 : }
634 : }
635 : }
636 :
637 : // write nodal data
638 0 : write_nodal_values_piece(File, u, time, grid, ssGrp, dim, numVert);
639 :
640 : // write cell data
641 0 : write_cell_values_piece(File, u, time, grid, ssGrp, dim, numElem);
642 :
643 : // write closing tag
644 0 : File << VTKFileWriter::normal;
645 0 : File << " </Piece>\n";
646 0 : }
647 :
648 : ////////////////////////////////////////////////////////////////////////////////
649 : // Sizes
650 : ////////////////////////////////////////////////////////////////////////////////
651 :
652 : template <int TDim>
653 : template <typename TElem, typename T>
654 0 : void VTKOutput<TDim>::
655 : count_sizes(Grid& grid, const T& iterContainer, const int si,
656 : int& numVert, int& numElem, int& numConn)
657 : {
658 : // get reference element
659 : typedef typename reference_element_traits<TElem>::reference_element_type
660 : ref_elem_type;
661 :
662 : // number of corners of element
663 : static const int numCo = ref_elem_type::numCorners;
664 :
665 : // get iterators
666 : typedef typename IteratorProvider<T>::template traits<TElem>::const_iterator const_iterator;
667 0 : const_iterator iterBegin = IteratorProvider<T>::template begin<TElem>(iterContainer, si);
668 0 : const_iterator iterEnd = IteratorProvider<T>::template end<TElem>(iterContainer, si);
669 :
670 : // loop elements
671 0 : for( ; iterBegin != iterEnd; ++iterBegin)
672 : {
673 : // get the element
674 : TElem *elem = *iterBegin;
675 :
676 : // count number of elements and number of connections;
677 : // handle octahedrons separately by splitting into a top and bottom pyramid
678 : if(ref_elem_type::REFERENCE_OBJECT_ID != ROID_OCTAHEDRON)
679 : {
680 0 : ++numElem;
681 0 : numConn += numCo;
682 : }
683 : else
684 : {
685 : // counting top and bottom pyramid
686 0 : numElem += 2;
687 0 : numConn += 10;
688 : }
689 :
690 : // loop vertices of the element
691 0 : for(int i = 0; i < numCo; ++i)
692 : {
693 : // get vertex of the element
694 0 : Vertex* v = GetVertex(elem,i);
695 :
696 : // if this vertex has already been counted, skip it
697 0 : if(grid.is_marked(v)) continue;
698 :
699 : // count vertex and mark it
700 0 : ++numVert;
701 : grid.mark(v);
702 : }
703 : }
704 0 : };
705 :
706 : /**
707 : * Counts VTK grid elements (points, full-dim. elements and connections) in a given
708 : * piece of the ug grid.
709 : */
710 : template <int TDim>
711 : template <typename T>
712 0 : void VTKOutput<TDim>::
713 : count_piece_sizes(Grid& grid, const T& iterContainer, const SubsetGroup& ssGrp, const int dim,
714 : int& numVert, int& numElem, int& numConn)
715 : {
716 0 : numVert = numElem = numConn = 0;
717 :
718 : // reset all marks (we use vertex marking to avoid counting vertices twice)
719 0 : grid.begin_marking();
720 :
721 0 : for(size_t i = 0; i < ssGrp.size(); i++)
722 : {
723 : const int si = ssGrp[i];
724 0 : switch(dim) // switch dimension
725 : {
726 0 : case 0: count_sizes<Vertex, T>(grid, iterContainer, si, numVert, numElem, numConn); break;
727 0 : case 1: count_sizes<RegularEdge, T>(grid, iterContainer, si, numVert, numElem, numConn);
728 0 : count_sizes<ConstrainingEdge, T>(grid, iterContainer, si, numVert, numElem, numConn); break;
729 0 : case 2: count_sizes<Triangle, T>(grid, iterContainer, si, numVert, numElem, numConn);
730 0 : count_sizes<Quadrilateral, T>(grid, iterContainer, si, numVert, numElem, numConn);
731 0 : count_sizes<ConstrainingTriangle, T>(grid, iterContainer, si, numVert, numElem, numConn);
732 0 : count_sizes<ConstrainingQuadrilateral, T>(grid, iterContainer, si, numVert, numElem, numConn); break;
733 0 : case 3: count_sizes<Tetrahedron, T>(grid, iterContainer, si, numVert, numElem, numConn);
734 0 : count_sizes<Pyramid, T>(grid, iterContainer, si, numVert, numElem, numConn);
735 0 : count_sizes<Prism, T>(grid, iterContainer, si, numVert, numElem, numConn);
736 0 : count_sizes<Octahedron, T>(grid, iterContainer, si, numVert, numElem, numConn);
737 0 : count_sizes<Hexahedron, T>(grid, iterContainer, si, numVert, numElem, numConn); break;
738 0 : default: UG_THROW("VTK::count_piece_sizes: Dimension " << dim <<
739 : " is not supported.");
740 : }
741 : }
742 :
743 : // signal end of marking
744 0 : grid.end_marking();
745 0 : }
746 :
747 : ////////////////////////////////////////////////////////////////////////////////
748 : // Points
749 : ////////////////////////////////////////////////////////////////////////////////
750 :
751 : /**
752 : * This method loops all elements of a subset and writes the vertex
753 : * coordinates to a file. All vertices of the element are written reguardless
754 : * if the vertex itself is contained in the subset or another. Written
755 : * vertices are marked and not written again.
756 : */
757 : template <int TDim>
758 : template <typename TElem, typename T>
759 0 : void VTKOutput<TDim>::
760 : write_points_elementwise(VTKFileWriter& File,
761 : Grid::VertexAttachmentAccessor<Attachment<int> >& aaVrtIndex,
762 : const Grid::VertexAttachmentAccessor<Attachment<MathVector<TDim> > >& aaPos,
763 : Grid& grid, const T& iterContainer, const int si, int& n)
764 : {
765 : // get reference element
766 : typedef typename reference_element_traits<TElem>::reference_element_type
767 : ref_elem_type;
768 : // get iterators
769 : typedef typename IteratorProvider<T>::template traits<TElem>::const_iterator const_iterator;
770 0 : const_iterator iterBegin = IteratorProvider<T>::template begin<TElem>(iterContainer, si);
771 0 : const_iterator iterEnd = IteratorProvider<T>::template end<TElem>(iterContainer, si);
772 0 : if(m_bBinary)
773 0 : File << VTKFileWriter::base64_binary;
774 : else
775 0 : File << VTKFileWriter::normal;
776 :
777 : // loop all elements of the subset
778 0 : for( ; iterBegin != iterEnd; ++iterBegin)
779 : {
780 : // get the element
781 : TElem *elem = *iterBegin;
782 :
783 : // loop vertices of the element
784 0 : for(size_t i = 0; i < (size_t) ref_elem_type::numCorners; ++i)
785 : {
786 : // get vertex of element
787 : Vertex* v = GetVertex(elem, i);
788 :
789 : // if vertex has already be handled, skip it
790 0 : if(grid.is_marked(v)) continue;
791 :
792 : // mark the vertex as processed
793 : grid.mark(v);
794 :
795 : // number vertex
796 0 : aaVrtIndex[v] = n++;
797 :
798 : // get position of vertex and write position to stream
799 0 : write_item_to_file(File, aaPos[v]);
800 : }
801 : }
802 0 : }
803 :
804 : /**
805 : * This function writes the vertices of a piece (the specified subsets) of
806 : * the grid to a vtk file.
807 : */
808 : template <int TDim>
809 : template <typename T>
810 0 : void VTKOutput<TDim>::
811 : write_points(VTKFileWriter& File,
812 : Grid::VertexAttachmentAccessor<Attachment<int> >& aaVrtIndex,
813 : const Grid::VertexAttachmentAccessor<Attachment<MathVector<TDim> > >& aaPos,
814 : Grid& grid, const T& iterContainer, const SubsetGroup& ssGrp, const int dim,
815 : const int numVert)
816 : {
817 0 : if(!m_bWriteGrid){
818 0 : return;
819 : }
820 :
821 : // write starting xml tag for points
822 0 : File << VTKFileWriter::normal;
823 0 : File << " <Points>\n";
824 0 : File << " <DataArray type=\"Float32\" NumberOfComponents=\"3\" format="
825 0 : << (m_bBinary ? "\"binary\"" : "\"ascii\"") << ">\n";
826 0 : int n = 3*sizeof(float) * numVert;
827 0 : if(m_bBinary)
828 0 : File << VTKFileWriter::base64_binary << n;
829 :
830 : // reset counter for vertices
831 0 : n = 0;
832 :
833 : // start marking of vertices
834 0 : grid.begin_marking();
835 :
836 : // switch dimension
837 0 : if(numVert > 0)
838 0 : for(size_t i = 0; i < ssGrp.size(); i++){
839 0 : switch(dim){
840 0 : case 0: write_points_elementwise<Vertex,T>(File, aaVrtIndex, aaPos, grid, iterContainer, ssGrp[i], n); break;
841 0 : case 1: write_points_elementwise<RegularEdge,T>(File, aaVrtIndex, aaPos, grid, iterContainer, ssGrp[i], n);
842 0 : write_points_elementwise<ConstrainingEdge,T>(File, aaVrtIndex, aaPos, grid, iterContainer, ssGrp[i], n); break;
843 0 : case 2: write_points_elementwise<Triangle,T>(File, aaVrtIndex, aaPos, grid, iterContainer, ssGrp[i], n);
844 0 : write_points_elementwise<Quadrilateral,T>(File, aaVrtIndex, aaPos, grid, iterContainer, ssGrp[i], n);
845 0 : write_points_elementwise<ConstrainingTriangle,T>(File, aaVrtIndex, aaPos, grid, iterContainer, ssGrp[i], n);
846 0 : write_points_elementwise<ConstrainingQuadrilateral,T>(File, aaVrtIndex, aaPos, grid, iterContainer, ssGrp[i], n); break;
847 0 : case 3: write_points_elementwise<Tetrahedron,T>(File, aaVrtIndex, aaPos, grid, iterContainer, ssGrp[i], n);
848 0 : write_points_elementwise<Pyramid,T>(File, aaVrtIndex, aaPos, grid, iterContainer, ssGrp[i], n);
849 0 : write_points_elementwise<Prism,T>(File, aaVrtIndex, aaPos, grid, iterContainer, ssGrp[i], n);
850 0 : write_points_elementwise<Octahedron,T>(File, aaVrtIndex, aaPos, grid, iterContainer, ssGrp[i], n);
851 0 : write_points_elementwise<Hexahedron,T>(File, aaVrtIndex, aaPos, grid, iterContainer, ssGrp[i], n); break;
852 0 : default: UG_THROW("VTK::write_points: Dimension " << dim <<
853 : " is not supported.");
854 : }
855 : }
856 :
857 : // signal end of marking the grid
858 0 : grid.end_marking();
859 :
860 : // write closing tags
861 0 : File << VTKFileWriter::normal;
862 0 : File << "\n </DataArray>\n";
863 0 : File << " </Points>\n";
864 : }
865 :
866 : ////////////////////////////////////////////////////////////////////////////////
867 : // Cells
868 : ////////////////////////////////////////////////////////////////////////////////
869 :
870 : template <int TDim>
871 0 : void VTKOutput<TDim>::
872 : write_cell_subset_names(VTKFileWriter& File, MGSubsetHandler &sh)
873 : {
874 0 : File << VTKFileWriter::normal;
875 0 : File << " <RegionInfo Name=\"regions\">\n";
876 0 : for(int i = 0; i < sh.num_subsets(); ++i){
877 0 : File << " <Region Name=\"" << sh.get_subset_name(i) << "\"></Region>\n";
878 : }
879 0 : File << " </RegionInfo>\n";
880 0 : }
881 :
882 : /**
883 : * This function writes the elements that are part of the specified subsets.
884 : */
885 : template <int TDim>
886 : template <typename T>
887 0 : void VTKOutput<TDim>::
888 : write_cells(VTKFileWriter& File,
889 : Grid::VertexAttachmentAccessor<Attachment<int> >& aaVrtIndex,
890 : Grid& grid, const T& iterContainer, const SubsetGroup& ssGrp, const int dim,
891 : const int numElem, const int numConn)
892 : {
893 0 : File << VTKFileWriter::normal;
894 :
895 : // write opening tag to indicate that elements will be written
896 0 : File << " <Cells>\n";
897 :
898 : // write connectivities of elements
899 0 : write_cell_connectivity(File, aaVrtIndex, grid, iterContainer, ssGrp, dim, numConn);
900 :
901 : // write offsets for elements (i.e. number of nodes counted up)
902 0 : write_cell_offsets(File, iterContainer, ssGrp, dim, numElem);
903 :
904 : // write a defined type for each cell
905 0 : write_cell_types(File, iterContainer, ssGrp, dim, numElem);
906 :
907 : // write closing tag
908 0 : File << VTKFileWriter::normal;
909 0 : File << " </Cells>\n";
910 0 : }
911 :
912 : ////////////////////////////////////////////////////////////////////////////////
913 : // Connectivity
914 : ////////////////////////////////////////////////////////////////////////////////
915 :
916 :
917 : template <int TDim>
918 : template <class TElem, typename T>
919 0 : void VTKOutput<TDim>::
920 : write_cell_connectivity(VTKFileWriter& File,
921 : Grid::VertexAttachmentAccessor<Attachment<int> >& aaVrtIndex,
922 : Grid& grid, const T& iterContainer, const int si)
923 : {
924 : // get reference element type
925 : typedef typename reference_element_traits<TElem>::reference_element_type
926 : ref_elem_type;
927 :
928 : // get reference element id
929 : static const ReferenceObjectID refID = ref_elem_type::REFERENCE_OBJECT_ID;
930 :
931 : // get iterators
932 : typedef typename IteratorProvider<T>::template traits<TElem>::const_iterator const_iterator;
933 0 : const_iterator iterBegin = IteratorProvider<T>::template begin<TElem>(iterContainer, si);
934 0 : const_iterator iterEnd = IteratorProvider<T>::template end<TElem>(iterContainer, si);
935 :
936 0 : if(m_bBinary)
937 0 : File << VTKFileWriter::base64_binary;
938 : else
939 0 : File << VTKFileWriter::normal;
940 :
941 : // loop all elements
942 0 : for( ; iterBegin != iterEnd; iterBegin++)
943 : {
944 : // get element
945 : TElem* elem = *iterBegin;
946 :
947 : // write ids of the element
948 : if(refID != ROID_PRISM && refID != ROID_OCTAHEDRON)
949 : {
950 0 : for(size_t i=0; i< (size_t) ref_elem_type::numCorners; i++)
951 : {
952 0 : Vertex* vert = elem->vertex(i);
953 0 : int id = aaVrtIndex[vert];
954 0 : File << id;
955 0 : if(!m_bBinary)
956 0 : File << ' ';
957 : }
958 : }
959 : else if(refID == ROID_PRISM)
960 : {
961 0 : int id = aaVrtIndex[elem->vertex(0)]; File << id;
962 0 : if(!m_bBinary) File << ' ';
963 0 : id = aaVrtIndex[elem->vertex(2)]; File << id;
964 0 : if(!m_bBinary) File << ' ';
965 0 : id = aaVrtIndex[elem->vertex(1)]; File << id;
966 0 : if(!m_bBinary) File << ' ';
967 0 : id = aaVrtIndex[elem->vertex(3)]; File << id;
968 0 : if(!m_bBinary) File << ' ';
969 0 : id = aaVrtIndex[elem->vertex(5)]; File << id;
970 0 : if(!m_bBinary) File << ' ';
971 0 : id = aaVrtIndex[elem->vertex(4)]; File << id;
972 0 : if(!m_bBinary) File << ' ';
973 : }
974 : // case == ROID_OCTAHEDRON (handle by a splitting into a top and bottom pyramid)
975 : else
976 : {
977 : // top pyramid
978 0 : int id = aaVrtIndex[elem->vertex(1)]; File << id;
979 0 : if(!m_bBinary) File << ' ';
980 0 : id = aaVrtIndex[elem->vertex(2)]; File << id;
981 0 : if(!m_bBinary) File << ' ';
982 0 : id = aaVrtIndex[elem->vertex(3)]; File << id;
983 0 : if(!m_bBinary) File << ' ';
984 0 : id = aaVrtIndex[elem->vertex(4)]; File << id;
985 0 : if(!m_bBinary) File << ' ';
986 0 : id = aaVrtIndex[elem->vertex(5)]; File << id;
987 0 : if(!m_bBinary) File << ' ';
988 :
989 : // bottom pyramid
990 0 : id = aaVrtIndex[elem->vertex(1)]; File << id;
991 0 : if(!m_bBinary) File << ' ';
992 0 : id = aaVrtIndex[elem->vertex(2)]; File << id;
993 0 : if(!m_bBinary) File << ' ';
994 0 : id = aaVrtIndex[elem->vertex(3)]; File << id;
995 0 : if(!m_bBinary) File << ' ';
996 0 : id = aaVrtIndex[elem->vertex(4)]; File << id;
997 0 : if(!m_bBinary) File << ' ';
998 0 : id = aaVrtIndex[elem->vertex(0)]; File << id;
999 0 : if(!m_bBinary) File << ' ';
1000 : }
1001 : }
1002 0 : }
1003 :
1004 : /**
1005 : * This method writes the 'connectivity' for each element of a subset. The
1006 : * connectivity are the indices of all vertex the element is formed of.
1007 : */
1008 : template <int TDim>
1009 : template <typename T>
1010 0 : void VTKOutput<TDim>::
1011 : write_cell_connectivity(VTKFileWriter& File,
1012 : Grid::VertexAttachmentAccessor<Attachment<int> >& aaVrtIndex,
1013 : Grid& grid, const T& iterContainer, const SubsetGroup& ssGrp,
1014 : const int dim, const int numConn)
1015 : {
1016 0 : File << VTKFileWriter::normal;
1017 : // write opening tag to indicate that connections will be written
1018 0 : File << " <DataArray type=\"Int32\" Name=\"connectivity\" format="
1019 0 : << (m_bBinary ? "\"binary\"": "\"ascii\"") << ">\n";
1020 0 : int n = sizeof(int) * numConn;
1021 :
1022 0 : if(m_bBinary)
1023 0 : File << VTKFileWriter::base64_binary << n;
1024 : // switch dimension
1025 0 : if(numConn > 0)
1026 0 : for(size_t i = 0; i < ssGrp.size(); i++){
1027 0 : switch(dim)
1028 : {
1029 : case 0: break; // no elements -> nothing to do
1030 0 : case 1: write_cell_connectivity<RegularEdge,T>(File, aaVrtIndex, grid, iterContainer, ssGrp[i]);
1031 0 : write_cell_connectivity<ConstrainingEdge,T>(File, aaVrtIndex, grid, iterContainer, ssGrp[i]); break;
1032 0 : case 2: write_cell_connectivity<Triangle,T>(File, aaVrtIndex, grid, iterContainer, ssGrp[i]);
1033 0 : write_cell_connectivity<Quadrilateral,T>(File, aaVrtIndex, grid, iterContainer, ssGrp[i]);
1034 0 : write_cell_connectivity<ConstrainingTriangle,T>(File, aaVrtIndex, grid, iterContainer, ssGrp[i]);
1035 0 : write_cell_connectivity<ConstrainingQuadrilateral,T>(File, aaVrtIndex, grid, iterContainer, ssGrp[i]); break;
1036 0 : case 3: write_cell_connectivity<Tetrahedron,T>(File, aaVrtIndex, grid, iterContainer, ssGrp[i]);
1037 0 : write_cell_connectivity<Pyramid,T>(File, aaVrtIndex, grid, iterContainer, ssGrp[i]);
1038 0 : write_cell_connectivity<Prism,T>(File, aaVrtIndex, grid, iterContainer, ssGrp[i]);
1039 0 : write_cell_connectivity<Octahedron,T>(File, aaVrtIndex, grid, iterContainer, ssGrp[i]);
1040 0 : write_cell_connectivity<Hexahedron,T>(File, aaVrtIndex, grid, iterContainer, ssGrp[i]); break;
1041 0 : default: UG_THROW("VTK::write_cell_connectivity: Dimension " << dim <<
1042 : " is not supported.");
1043 : }
1044 : }
1045 :
1046 : // write closing tag
1047 0 : File << VTKFileWriter::normal;
1048 0 : File << "\n </DataArray>\n";
1049 0 : }
1050 :
1051 : ////////////////////////////////////////////////////////////////////////////////
1052 : // Offset
1053 : ////////////////////////////////////////////////////////////////////////////////
1054 :
1055 :
1056 : template <int TDim>
1057 : template <class TElem, typename T>
1058 0 : void VTKOutput<TDim>::
1059 : write_cell_offsets(VTKFileWriter& File, const T& iterContainer, const int si, int& n)
1060 : {
1061 : // get reference element
1062 : typedef typename reference_element_traits<TElem>::reference_element_type
1063 : ref_elem_type;
1064 :
1065 : // get iterators
1066 : typedef typename IteratorProvider<T>::template traits<TElem>::const_iterator const_iterator;
1067 0 : const_iterator iterBegin = IteratorProvider<T>::template begin<TElem>(iterContainer, si);
1068 0 : const_iterator iterEnd = IteratorProvider<T>::template end<TElem>(iterContainer, si);
1069 :
1070 0 : if(m_bBinary)
1071 0 : File << VTKFileWriter::base64_binary;
1072 : else
1073 0 : File << VTKFileWriter::normal;
1074 :
1075 : // loop all elements
1076 0 : for( ; iterBegin != iterEnd; ++iterBegin)
1077 : {
1078 : // handle octahedron separately by splitting into a top and bottom pyramid
1079 : if(ref_elem_type::REFERENCE_OBJECT_ID != ROID_OCTAHEDRON)
1080 : {
1081 : // increase counter of vertices
1082 0 : n += ref_elem_type::numCorners;
1083 :
1084 : // write offset
1085 0 : File << n;
1086 0 : if(!m_bBinary)
1087 0 : File << ' ';
1088 : }
1089 : else
1090 : {
1091 : // increase counter for top pyramid vertices
1092 0 : n += 5;
1093 :
1094 : // write offset for top pyramid
1095 0 : File << n;
1096 0 : if(!m_bBinary)
1097 0 : File << ' ';
1098 :
1099 : // increase counter for bottom pyramid vertices
1100 0 : n += 5;
1101 :
1102 : // write offset for bottom pyramid
1103 0 : File << n;
1104 0 : if(!m_bBinary)
1105 0 : File << ' ';
1106 : }
1107 : }
1108 0 : }
1109 :
1110 : /**
1111 : * This method writes the 'offset' for each element of a subset.
1112 : */
1113 : template <int TDim>
1114 : template <typename T>
1115 0 : void VTKOutput<TDim>::
1116 : write_cell_offsets(VTKFileWriter& File, const T& iterContainer, const SubsetGroup& ssGrp,
1117 : const int dim, const int numElem)
1118 : {
1119 0 : File << VTKFileWriter::normal;
1120 : // write opening tag indicating that offsets are going to be written
1121 0 : File << " <DataArray type=\"Int32\" Name=\"offsets\" format="
1122 0 : << (m_bBinary ? "\"binary\"": "\"ascii\"") << ">\n";
1123 0 : int n = sizeof(int) * numElem;
1124 0 : if(m_bBinary)
1125 0 : File << VTKFileWriter::base64_binary << n;
1126 :
1127 0 : n = 0;
1128 : // switch dimension
1129 0 : if(numElem > 0)
1130 0 : for(size_t i = 0; i < ssGrp.size(); i++){
1131 0 : switch(dim)
1132 : {
1133 : case 0: break; // no elements -> nothing to do
1134 0 : case 1: write_cell_offsets<RegularEdge,T>(File, iterContainer, ssGrp[i], n);
1135 0 : write_cell_offsets<ConstrainingEdge,T>(File, iterContainer, ssGrp[i], n); break;
1136 0 : case 2: write_cell_offsets<Triangle,T>(File, iterContainer, ssGrp[i], n);
1137 0 : write_cell_offsets<Quadrilateral,T>(File, iterContainer, ssGrp[i], n);
1138 0 : write_cell_offsets<ConstrainingTriangle,T>(File, iterContainer, ssGrp[i], n);
1139 0 : write_cell_offsets<ConstrainingQuadrilateral,T>(File, iterContainer, ssGrp[i], n); break;
1140 0 : case 3: write_cell_offsets<Tetrahedron,T>(File, iterContainer, ssGrp[i], n);
1141 0 : write_cell_offsets<Pyramid,T>(File, iterContainer, ssGrp[i], n);
1142 0 : write_cell_offsets<Prism,T>(File, iterContainer, ssGrp[i], n);
1143 0 : write_cell_offsets<Octahedron,T>(File, iterContainer, ssGrp[i], n);
1144 0 : write_cell_offsets<Hexahedron,T>(File, iterContainer, ssGrp[i], n); break;
1145 0 : default: UG_THROW("VTK::write_cell_offsets: Dimension " << dim <<
1146 : " is not supported.");
1147 : }
1148 : }
1149 :
1150 : // closing tag
1151 0 : File << VTKFileWriter::normal;
1152 0 : File << "\n </DataArray>\n";
1153 0 : }
1154 :
1155 : ////////////////////////////////////////////////////////////////////////////////
1156 : // Types
1157 : ////////////////////////////////////////////////////////////////////////////////
1158 :
1159 :
1160 : template <int TDim>
1161 : template <class TElem, typename T>
1162 0 : void VTKOutput<TDim>::
1163 : write_cell_types(VTKFileWriter& File, const T& iterContainer, const int si)
1164 : {
1165 : // get reference element
1166 : typedef typename reference_element_traits<TElem>::reference_element_type
1167 : ref_elem_type;
1168 : // get object type
1169 : static const ReferenceObjectID refID = ref_elem_type::REFERENCE_OBJECT_ID;
1170 :
1171 : // type
1172 : char type;
1173 :
1174 : // get type, based on reference element type
1175 : switch(refID)
1176 : {
1177 : case ROID_EDGE: type = (char) 3; break;
1178 : case ROID_TRIANGLE: type = (char) 5; break;
1179 : case ROID_QUADRILATERAL: type = (char) 9; break;
1180 : case ROID_TETRAHEDRON: type = (char) 10; break;
1181 : case ROID_PYRAMID: type = (char) 14; break;
1182 : case ROID_PRISM: type = (char) 13; break;
1183 : case ROID_OCTAHEDRON: type = (char) 14; break; // use type id == 14 (PYRAMID) as we handle an octahedron by splitting into a top and bottom pyramid
1184 : case ROID_HEXAHEDRON: type = (char) 12; break;
1185 : default: UG_THROW("Element Type not known.");
1186 : }
1187 :
1188 : // get iterators
1189 : typedef typename IteratorProvider<T>::template traits<TElem>::const_iterator const_iterator;
1190 0 : const_iterator iterBegin = IteratorProvider<T>::template begin<TElem>(iterContainer, si);
1191 0 : const_iterator iterEnd = IteratorProvider<T>::template end<TElem>(iterContainer, si);
1192 :
1193 0 : if(m_bBinary)
1194 0 : File << VTKFileWriter::base64_binary;
1195 : else
1196 0 : File << VTKFileWriter::normal;
1197 : // loop all elements, write type for each element to stream
1198 0 : for( ; iterBegin != iterEnd; ++iterBegin)
1199 : {
1200 : // handle octahedral case separately by splitting into a top and bottom pyramid
1201 : if(ref_elem_type::REFERENCE_OBJECT_ID != ROID_OCTAHEDRON)
1202 : {
1203 0 : if(m_bBinary)
1204 0 : File << type;
1205 : else
1206 0 : File << (int) type << ' ';
1207 : }
1208 : else
1209 : {
1210 0 : if(m_bBinary)
1211 : {
1212 0 : File << type; // top pyramid
1213 0 : File << type; // bottom pyramid
1214 : }
1215 : else
1216 : {
1217 0 : File << (int) type << ' '; // top pyramid
1218 0 : File << (int) type << ' '; // bottom pyramid
1219 : }
1220 : }
1221 : }
1222 0 : }
1223 :
1224 : /**
1225 : * This method writes the 'type' for each element of a subset. The type is
1226 : * a vtk-format specified number for a given element type.
1227 : */
1228 : template <int TDim>
1229 : template <typename T>
1230 0 : void VTKOutput<TDim>::
1231 : write_cell_types(VTKFileWriter& File, const T& iterContainer, const SubsetGroup& ssGrp,
1232 : const int dim, int numElem)
1233 : {
1234 0 : File << VTKFileWriter::normal;
1235 : // write opening tag to indicate that types will be written
1236 0 : File << " <DataArray type=\"Int8\" Name=\"types\" format="
1237 0 : << (m_bBinary ? "\"binary\"": "\"ascii\"") << ">\n";
1238 0 : if(m_bBinary)
1239 0 : File << VTKFileWriter::base64_binary << numElem;
1240 :
1241 : // switch dimension
1242 0 : if(numElem > 0)
1243 0 : for(size_t i = 0; i < ssGrp.size(); i++)
1244 : {
1245 0 : switch(dim)
1246 : {
1247 : case 0: break; // no elements -> nothing to do
1248 0 : case 1: write_cell_types<RegularEdge>(File, iterContainer, ssGrp[i]);
1249 0 : write_cell_types<ConstrainingEdge>(File, iterContainer, ssGrp[i]); break;
1250 0 : case 2: write_cell_types<Triangle>(File, iterContainer, ssGrp[i]);
1251 0 : write_cell_types<Quadrilateral>(File, iterContainer, ssGrp[i]);
1252 0 : write_cell_types<ConstrainingTriangle>(File, iterContainer, ssGrp[i]);
1253 0 : write_cell_types<ConstrainingQuadrilateral>(File, iterContainer, ssGrp[i]); break;
1254 0 : case 3: write_cell_types<Tetrahedron>(File, iterContainer, ssGrp[i]);
1255 0 : write_cell_types<Pyramid>(File, iterContainer, ssGrp[i]);
1256 0 : write_cell_types<Prism>(File, iterContainer, ssGrp[i]);
1257 0 : write_cell_types<Octahedron>(File, iterContainer, ssGrp[i]);
1258 0 : write_cell_types<Hexahedron>(File, iterContainer, ssGrp[i]); break;
1259 0 : default: UG_THROW("VTK::write_cell_types: Dimension " << dim <<
1260 : " is not supported.");
1261 : }
1262 : }
1263 :
1264 : // write closing tag
1265 0 : File << VTKFileWriter::normal;
1266 0 : File << "\n </DataArray>\n";
1267 0 : }
1268 :
1269 : ////////////////////////////////////////////////////////////////////////////////
1270 : // Subset Indices
1271 : ////////////////////////////////////////////////////////////////////////////////
1272 :
1273 :
1274 : template <int TDim>
1275 : template <class TElem, typename T>
1276 : void VTKOutput<TDim>::
1277 : write_cell_subsets(VTKFileWriter& File, const T& iterContainer, const int si, MGSubsetHandler& sh)
1278 : {
1279 : // subset
1280 : int subset;
1281 :
1282 : // get iterators
1283 : typedef typename IteratorProvider<T>::template traits<TElem>::const_iterator const_iterator;
1284 : const_iterator iterBegin = IteratorProvider<T>::template begin<TElem>(iterContainer, si);
1285 : const_iterator iterEnd = IteratorProvider<T>::template end<TElem>(iterContainer, si);
1286 :
1287 : if(m_bBinary)
1288 : File << VTKFileWriter::base64_binary;
1289 : else
1290 : File << VTKFileWriter::normal;
1291 : // loop all elements, write type for each element to stream
1292 : for( ; iterBegin != iterEnd; ++iterBegin)
1293 : {
1294 : subset = (char)sh.get_subset_index(*iterBegin);
1295 : //iterContainer.get_subset_name(2);
1296 :
1297 : if(m_bBinary){
1298 : File << (char) subset << ' ';
1299 : }
1300 : else{
1301 : File << subset << ' ';
1302 : }
1303 : }
1304 : }
1305 :
1306 : /**
1307 : * This method writes the 'region' for each element of a subset. The region is
1308 : * the ug4 subset index for a given element.
1309 : */
1310 : template <int TDim>
1311 : template <typename T>
1312 : void VTKOutput<TDim>::
1313 : write_cell_subsets(VTKFileWriter& File, const T& iterContainer, const SubsetGroup& ssGrp,
1314 : const int dim, const int numElem, MGSubsetHandler& sh)
1315 : {
1316 : File << VTKFileWriter::normal;
1317 : // write opening tag to indicate that types will be written
1318 : File << " <DataArray type=\"Int8\" Name=\"regions\" format="
1319 : << (m_bBinary ? "\"binary\"": "\"ascii\"") << ">\n";
1320 : if(m_bBinary)
1321 : File << VTKFileWriter::base64_binary << numElem;
1322 :
1323 : // switch dimension
1324 : if(numElem > 0)
1325 : for(size_t i = 0; i < ssGrp.size(); i++)
1326 : {
1327 : switch(dim)
1328 : {
1329 : case 0: break; // no elements -> nothing to do
1330 : case 1: write_cell_subsets<RegularEdge>(File, iterContainer, ssGrp[i], sh);
1331 : write_cell_subsets<ConstrainingEdge>(File, iterContainer, ssGrp[i], sh); break;
1332 : case 2: write_cell_subsets<Triangle>(File, iterContainer, ssGrp[i], sh);
1333 : write_cell_subsets<Quadrilateral>(File, iterContainer, ssGrp[i], sh);
1334 : write_cell_subsets<ConstrainingTriangle>(File, iterContainer, ssGrp[i], sh);
1335 : write_cell_subsets<ConstrainingQuadrilateral>(File, iterContainer, ssGrp[i], sh); break;
1336 : case 3: write_cell_subsets<Tetrahedron>(File, iterContainer, ssGrp[i], sh);
1337 : write_cell_subsets<Pyramid>(File, iterContainer, ssGrp[i], sh);
1338 : write_cell_subsets<Prism>(File, iterContainer, ssGrp[i], sh);
1339 : write_cell_subsets<Octahedron>(File, iterContainer, ssGrp[i], sh);
1340 : write_cell_subsets<Hexahedron>(File, iterContainer, ssGrp[i], sh); break;
1341 : default: UG_THROW("VTK::write_cell_subsets: Dimension " << dim <<
1342 : " is not supported.");
1343 : }
1344 : }
1345 :
1346 : // write closing tag
1347 : File << VTKFileWriter::normal;
1348 : File << "\n </DataArray>\n";
1349 : }
1350 :
1351 : ////////////////////////////////////////////////////////////////////////////////
1352 : // Subset Indices
1353 : ////////////////////////////////////////////////////////////////////////////////
1354 :
1355 :
1356 : template <int TDim>
1357 : template <class TElem, typename T>
1358 : void VTKOutput<TDim>::
1359 : write_cell_proc_ranks(VTKFileWriter& File, const T& iterContainer, const int si, MGSubsetHandler& sh)
1360 : {
1361 : // subset
1362 : int rank = 0;
1363 :
1364 : #ifdef UG_PARALLEL
1365 : rank = pcl::ProcRank();
1366 : #endif
1367 :
1368 : // get iterators
1369 : typedef typename IteratorProvider<T>::template traits<TElem>::const_iterator const_iterator;
1370 : const_iterator iterBegin = IteratorProvider<T>::template begin<TElem>(iterContainer, si);
1371 : const_iterator iterEnd = IteratorProvider<T>::template end<TElem>(iterContainer, si);
1372 :
1373 : if(m_bBinary)
1374 : File << VTKFileWriter::base64_binary;
1375 : else
1376 : File << VTKFileWriter::normal;
1377 : // loop all elements, write type for each element to stream
1378 : for( ; iterBegin != iterEnd; ++iterBegin)
1379 : {
1380 :
1381 : if(m_bBinary){
1382 : File << (char) rank << ' ';
1383 : }
1384 : else{
1385 : File << rank << ' ';
1386 : }
1387 : }
1388 : }
1389 :
1390 : /**
1391 : * This method writes the proc ranks for each cell.
1392 : */
1393 : template <int TDim>
1394 : template <typename T>
1395 : void VTKOutput<TDim>::
1396 : write_cell_proc_ranks(VTKFileWriter& File, const T& iterContainer, const SubsetGroup& ssGrp,
1397 : const int dim, const int numElem, MGSubsetHandler& sh)
1398 : {
1399 : File << VTKFileWriter::normal;
1400 : // write opening tag to indicate that types will be written
1401 : File << " <DataArray type=\"Int8\" Name=\"proc_ranks\" format="
1402 : << (m_bBinary ? "\"binary\"": "\"ascii\"") << ">\n";
1403 : if(m_bBinary)
1404 : File << VTKFileWriter::base64_binary << numElem;
1405 :
1406 : // switch dimension
1407 : if(numElem > 0)
1408 : for(size_t i = 0; i < ssGrp.size(); i++)
1409 : {
1410 : switch(dim)
1411 : {
1412 : case 0: break; // no elements -> nothing to do
1413 : case 1: write_cell_proc_ranks<RegularEdge>(File, iterContainer, ssGrp[i], sh);
1414 : write_cell_proc_ranks<ConstrainingEdge>(File, iterContainer, ssGrp[i], sh); break;
1415 : case 2: write_cell_proc_ranks<Triangle>(File, iterContainer, ssGrp[i], sh);
1416 : write_cell_proc_ranks<Quadrilateral>(File, iterContainer, ssGrp[i], sh);
1417 : write_cell_proc_ranks<ConstrainingTriangle>(File, iterContainer, ssGrp[i], sh);
1418 : write_cell_proc_ranks<ConstrainingQuadrilateral>(File, iterContainer, ssGrp[i], sh); break;
1419 : case 3: write_cell_proc_ranks<Tetrahedron>(File, iterContainer, ssGrp[i], sh);
1420 : write_cell_proc_ranks<Pyramid>(File, iterContainer, ssGrp[i], sh);
1421 : write_cell_proc_ranks<Prism>(File, iterContainer, ssGrp[i], sh);
1422 : write_cell_proc_ranks<Octahedron>(File, iterContainer, ssGrp[i], sh);
1423 : write_cell_proc_ranks<Hexahedron>(File, iterContainer, ssGrp[i], sh); break;
1424 : default: UG_THROW("VTK::write_cell_proc_ranks: Dimension " << dim <<
1425 : " is not supported.");
1426 : }
1427 : }
1428 :
1429 : // write closing tag
1430 : File << VTKFileWriter::normal;
1431 : File << "\n </DataArray>\n";
1432 : }
1433 :
1434 : ////////////////////////////////////////////////////////////////////////////////
1435 : // Nodal Value
1436 : ////////////////////////////////////////////////////////////////////////////////
1437 :
1438 : template <int TDim>
1439 : template <typename TElem, typename TFunction, typename TData>
1440 0 : void VTKOutput<TDim>::
1441 : write_nodal_data_elementwise(VTKFileWriter& File, TFunction& u, number time,
1442 : SmartPtr<UserData<TData, TDim> > spData,
1443 : Grid& grid, const int si)
1444 : {
1445 : // get reference element
1446 : typedef typename reference_element_traits<TElem>::reference_element_type
1447 : ref_elem_type;
1448 0 : static const ref_elem_type& refElem = Provider<ref_elem_type>::get();
1449 : static const size_t numCo = ref_elem_type::numCorners;
1450 :
1451 0 : if(m_bBinary)
1452 0 : File << VTKFileWriter::base64_binary;
1453 : else
1454 0 : File << VTKFileWriter::normal;
1455 :
1456 : // get iterators
1457 : typedef typename IteratorProvider<TFunction>::template traits<TElem>::const_iterator const_iterator;
1458 0 : const_iterator iterBegin = IteratorProvider<TFunction>::template begin<TElem>(u, si);
1459 0 : const_iterator iterEnd = IteratorProvider<TFunction>::template end<TElem>(u, si);
1460 :
1461 0 : std::vector<MathVector<TDim> > vCorner(numCo);
1462 :
1463 0 : TData vValue[numCo];
1464 0 : bool bNeedFct = spData->requires_grid_fct();
1465 :
1466 : // loop all elements
1467 0 : for( ; iterBegin != iterEnd; ++iterBegin)
1468 : {
1469 : // get element
1470 : TElem *elem = *iterBegin;
1471 :
1472 : // update the reference mapping for the corners
1473 0 : CollectCornerCoordinates(vCorner, *elem, u.domain()->position_accessor(), true);
1474 :
1475 : // get subset
1476 : int theSI = si;
1477 0 : if(si < 0) theSI = u.domain()->subset_handler()->get_subset_index(elem);
1478 :
1479 : // get local solution if needed
1480 0 : if(bNeedFct)
1481 : {
1482 : // create storage
1483 : LocalIndices ind;
1484 : LocalVector locU;
1485 :
1486 : // get global indices
1487 : u.indices(elem, ind);
1488 :
1489 : // adapt local algebra
1490 0 : locU.resize(ind);
1491 :
1492 : // read local values of u
1493 0 : GetLocalVector(locU, u);
1494 :
1495 : // compute data
1496 : try{
1497 0 : (*spData)(vValue, &vCorner[0], time, theSI, elem,
1498 : &vCorner[0], refElem.corners(), numCo, &locU, NULL);
1499 : }
1500 0 : UG_CATCH_THROW("VTK::write_nodal_data_elementwise: Cannot evaluate data.");
1501 : }
1502 : else
1503 : {
1504 : // compute data
1505 : try{
1506 0 : (*spData)(vValue, &vCorner[0], time, theSI, numCo);
1507 : }
1508 0 : UG_CATCH_THROW("VTK::write_nodal_data_elementwise: Cannot evaluate data.");
1509 : }
1510 :
1511 : // loop vertices of element
1512 0 : for(size_t co = 0; co < numCo; ++co)
1513 : {
1514 : // get vertex of element
1515 : Vertex* v = GetVertex(elem, co);
1516 :
1517 : // if vertex has been handled before, skip
1518 0 : if(grid.is_marked(v)) continue;
1519 :
1520 : // mark as used
1521 : grid.mark(v);
1522 :
1523 : // loop all components
1524 0 : write_item_to_file(File, vValue[co]);
1525 : }
1526 : }
1527 :
1528 0 : }
1529 :
1530 : /**
1531 : * Writes the nodal data to the vtu file.
1532 : */
1533 : template <int TDim>
1534 : template <typename TFunction, typename TData>
1535 0 : void VTKOutput<TDim>::
1536 : write_nodal_data(VTKFileWriter& File, TFunction& u, number time,
1537 : SmartPtr<UserData<TData, TDim> > spData,
1538 : const int numCmp,
1539 : const std::string& name,
1540 : Grid& grid, const SubsetGroup& ssGrp, const int dim, const int numVert)
1541 : {
1542 0 : spData->set_function_pattern(u.function_pattern());
1543 :
1544 : // check that nodal data is possible
1545 0 : if(!spData->continuous())
1546 0 : UG_THROW("VTK: data with name '"<<name<<"' is scheduled for nodal output,"
1547 : " but the data is not continuous. Cannot write it as nodal data.")
1548 :
1549 : // write opening tag
1550 0 : File << VTKFileWriter::normal;
1551 0 : File << " <DataArray type=\"Float32\" Name=\""<<name<<"\" "
1552 0 : "NumberOfComponents=\""<<numCmp<<"\" format="
1553 0 : << (m_bBinary ? "\"binary\"": "\"ascii\"") << ">\n";
1554 :
1555 0 : int n = sizeof(float) * numVert * numCmp;
1556 0 : if(m_bBinary)
1557 0 : File << VTKFileWriter::base64_binary << n;
1558 :
1559 : // start marking of grid
1560 0 : grid.begin_marking();
1561 :
1562 : // switch dimension
1563 0 : for(size_t i = 0; i < ssGrp.size(); i++)
1564 0 : switch(dim)
1565 : {
1566 0 : case 1: write_nodal_data_elementwise<RegularEdge,TFunction,TData>(File, u, time, spData, grid, ssGrp[i]);
1567 0 : write_nodal_data_elementwise<ConstrainingEdge,TFunction,TData>(File, u, time, spData, grid, ssGrp[i]); break;
1568 0 : case 2: write_nodal_data_elementwise<Triangle,TFunction,TData>(File, u, time, spData, grid, ssGrp[i]);
1569 0 : write_nodal_data_elementwise<Quadrilateral,TFunction,TData>(File, u, time, spData, grid, ssGrp[i]);
1570 0 : write_nodal_data_elementwise<ConstrainingTriangle,TFunction,TData>(File, u, time, spData, grid, ssGrp[i]);
1571 0 : write_nodal_data_elementwise<ConstrainingQuadrilateral,TFunction,TData>(File, u, time, spData, grid, ssGrp[i]); break;
1572 0 : case 3: write_nodal_data_elementwise<Tetrahedron,TFunction,TData>(File, u, time, spData, grid, ssGrp[i]);
1573 0 : write_nodal_data_elementwise<Pyramid,TFunction,TData>(File, u, time, spData, grid, ssGrp[i]);
1574 0 : write_nodal_data_elementwise<Prism,TFunction,TData>(File, u, time, spData, grid, ssGrp[i]);
1575 0 : write_nodal_data_elementwise<Octahedron,TFunction,TData>(File, u, time, spData, grid, ssGrp[i]);
1576 0 : write_nodal_data_elementwise<Hexahedron,TFunction,TData>(File, u, time, spData, grid, ssGrp[i]); break;
1577 0 : default: UG_THROW("VTK::write_nodal_data: Dimension " << dim <<
1578 : " is not supported.");
1579 : }
1580 :
1581 : // end marking
1582 0 : grid.end_marking();
1583 :
1584 : // write closing tag
1585 0 : File << VTKFileWriter::normal;
1586 0 : File << "\n </DataArray>\n";
1587 0 : };
1588 :
1589 :
1590 : template <int TDim>
1591 : template <typename TElem, typename TFunction>
1592 0 : void VTKOutput<TDim>::
1593 : write_nodal_values_elementwise(VTKFileWriter& File, TFunction& u,
1594 : const std::vector<size_t>& vFct,
1595 : Grid& grid, const int si)
1596 : {
1597 : // get reference element
1598 : typedef typename reference_element_traits<TElem>::reference_element_type
1599 : ref_elem_type;
1600 0 : if(m_bBinary)
1601 0 : File << VTKFileWriter::base64_binary;
1602 : else
1603 0 : File << VTKFileWriter::normal;
1604 :
1605 : // index vector
1606 : std::vector<DoFIndex> vMultInd;
1607 :
1608 : // get iterators
1609 : typedef typename IteratorProvider<TFunction>::template traits<TElem>::const_iterator const_iterator;
1610 0 : const_iterator iterBegin = IteratorProvider<TFunction>::template begin<TElem>(u, si);
1611 0 : const_iterator iterEnd = IteratorProvider<TFunction>::template end<TElem>(u, si);
1612 :
1613 : // loop all elements
1614 0 : for( ; iterBegin != iterEnd; ++iterBegin)
1615 : {
1616 : // get element
1617 : TElem *elem = *iterBegin;
1618 :
1619 : // loop vertices of element
1620 0 : for(size_t co = 0; co < (size_t) ref_elem_type::numCorners; ++co)
1621 : {
1622 : // get vertex of element
1623 : Vertex* v = GetVertex(elem, co);
1624 :
1625 : // if vertex has been handled before, skip
1626 0 : if(grid.is_marked(v)) continue;
1627 :
1628 : // mark as used
1629 : grid.mark(v);
1630 :
1631 : // loop all components
1632 0 : for(size_t i = 0; i < vFct.size(); ++i)
1633 : {
1634 : // get multi index of vertex for the function
1635 0 : if(u.inner_dof_indices(v, vFct[i], vMultInd) != 1)
1636 0 : UG_THROW("VTK:write_nodal_values_elementwise: "
1637 : "The function component "<<vFct[i]<<" has "<<
1638 : vMultInd.size()<<" DoFs in a vertex of an element of subset "
1639 : <<si<<". To write a component to vtk, exactly "
1640 : "one DoF must be given in any vertex.");
1641 :
1642 : // flush stream
1643 0 : write_item_to_file(File, DoFRef(u, vMultInd[0]));
1644 : }
1645 :
1646 :
1647 :
1648 : // fill with zeros up to 3d if vector type
1649 0 : if(vFct.size() != 1) {
1650 0 : for(size_t i = vFct.size(); i < 3; ++i) {
1651 0 : write_item_to_file(File, 0.f);
1652 : }
1653 : }
1654 : }
1655 : }
1656 :
1657 0 : }
1658 :
1659 : template <int TDim>
1660 : template <typename TFunction>
1661 0 : void VTKOutput<TDim>::
1662 : write_nodal_values(VTKFileWriter& File, TFunction& u,
1663 : const std::vector<size_t>& vFct, const std::string& name,
1664 : Grid& grid, const SubsetGroup& ssGrp, const int dim, const int numVert)
1665 : {
1666 0 : File << VTKFileWriter::normal;
1667 : // write opening tag
1668 0 : File << " <DataArray type=\"Float32\" Name=\""<<name<<"\" "
1669 0 : "NumberOfComponents=\""<<(vFct.size() == 1 ? 1 : 3)<<"\" format="
1670 0 : << (m_bBinary ? "\"binary\"": "\"ascii\"") << ">\n";
1671 :
1672 0 : int n = sizeof(float) * numVert * (vFct.size() == 1 ? 1 : 3);
1673 0 : if(m_bBinary)
1674 0 : File << VTKFileWriter::base64_binary << n;
1675 :
1676 : // start marking of grid
1677 0 : grid.begin_marking();
1678 :
1679 : // switch dimension
1680 0 : for(size_t i = 0; i < ssGrp.size(); i++)
1681 0 : switch(dim)
1682 : {
1683 0 : case 0: write_nodal_values_elementwise<Vertex>(File, u, vFct, grid, ssGrp[i]); break;
1684 0 : case 1: write_nodal_values_elementwise<RegularEdge>(File, u, vFct, grid, ssGrp[i]);
1685 0 : write_nodal_values_elementwise<ConstrainingEdge>(File, u, vFct, grid, ssGrp[i]); break;
1686 0 : case 2: write_nodal_values_elementwise<Triangle>(File, u, vFct, grid, ssGrp[i]);
1687 0 : write_nodal_values_elementwise<Quadrilateral>(File, u, vFct, grid, ssGrp[i]);
1688 0 : write_nodal_values_elementwise<ConstrainingTriangle>(File, u, vFct, grid, ssGrp[i]);
1689 0 : write_nodal_values_elementwise<ConstrainingQuadrilateral>(File, u, vFct, grid, ssGrp[i]); break;
1690 0 : case 3: write_nodal_values_elementwise<Tetrahedron>(File, u, vFct, grid, ssGrp[i]);
1691 0 : write_nodal_values_elementwise<Pyramid>(File, u, vFct, grid, ssGrp[i]);
1692 0 : write_nodal_values_elementwise<Prism>(File, u, vFct, grid, ssGrp[i]);
1693 0 : write_nodal_values_elementwise<Octahedron>(File, u, vFct, grid, ssGrp[i]);
1694 0 : write_nodal_values_elementwise<Hexahedron>(File, u, vFct, grid, ssGrp[i]); break;
1695 0 : default: UG_THROW("VTK::write_nodal_values: Dimension " << dim <<
1696 : " is not supported.");
1697 : }
1698 :
1699 : // end marking
1700 0 : grid.end_marking();
1701 :
1702 : // write closing tag
1703 0 : File << VTKFileWriter::normal;
1704 0 : File << "\n </DataArray>\n";
1705 0 : };
1706 :
1707 : template <int TDim>
1708 : template <typename TFunction>
1709 0 : void VTKOutput<TDim>::
1710 : write_nodal_values_piece(VTKFileWriter& File, TFunction& u, number time, Grid& grid,
1711 : const SubsetGroup& ssGrp, const int dim, const int numVert)
1712 : {
1713 0 : if(!m_vSymbFct.empty()){
1714 0 : for(std::map<std::string, std::vector<std::string> >::const_iterator iter =
1715 0 : m_vSymbFct.begin(); iter != m_vSymbFct.end(); ++iter){
1716 0 : const std::vector<std::string>& symbNames = (*iter).second;
1717 0 : const std::string& vtkName = (*iter).first;
1718 :
1719 : bool bContinuous = true;
1720 0 : for(size_t i = 0; i < symbNames.size(); ++i){
1721 : size_t fct = u.fct_id_by_name(symbNames[i].c_str());
1722 0 : if(!LocalFiniteElementProvider::continuous(u.local_finite_element_id(fct))){
1723 : bContinuous = false; break;
1724 : }
1725 : }
1726 :
1727 0 : if(bContinuous){
1728 0 : m_vSymbFctNodal[vtkName] = symbNames;
1729 : }
1730 : }
1731 : }
1732 :
1733 : // check if something to do
1734 0 : if(m_vSymbFctNodal.empty() && m_vScalarNodalData.empty() && m_vVectorNodalData.empty())
1735 : return;
1736 :
1737 : // write opening tag to indicate point data
1738 0 : File << VTKFileWriter::normal;
1739 0 : File << " <PointData>\n";
1740 :
1741 : // loop all selected symbolic names
1742 0 : for(std::map<std::string, std::vector<std::string> >::const_iterator iter =
1743 0 : m_vSymbFctNodal.begin(); iter != m_vSymbFctNodal.end(); ++iter)
1744 : {
1745 : // get symb function
1746 : const std::vector<std::string>& symbNames = (*iter).second;
1747 0 : const std::string& vtkName = (*iter).first;
1748 :
1749 : // create function group
1750 0 : std::vector<size_t> fctGrp(symbNames.size());
1751 0 : for(size_t i = 0; i < symbNames.size(); ++i)
1752 0 : fctGrp[i] = u.fct_id_by_name(symbNames[i].c_str());
1753 :
1754 : // check that all functions are contained in subset
1755 : bool bContained = true;
1756 0 : for(size_t i = 0; i < fctGrp.size(); ++i)
1757 0 : for (size_t j = 0; j < ssGrp.size(); ++j)
1758 0 : if(!u.is_def_in_subset(fctGrp[i], ssGrp[j]))
1759 : bContained = false;
1760 :
1761 0 : if(!bContained) continue;
1762 :
1763 : // write scalar value of this function
1764 : try{
1765 0 : write_nodal_values(File, u, fctGrp, vtkName, grid, ssGrp, dim, numVert);
1766 : }
1767 0 : UG_CATCH_THROW("VTK::write_piece: Failed write nodal Values.");
1768 : }
1769 :
1770 : // loop all scalar data
1771 0 : for(ScalarDataIterator iter = m_vScalarNodalData.begin();
1772 0 : iter != m_vScalarNodalData.end(); ++iter)
1773 : {
1774 : // get symb function
1775 : SmartPtr<UserData<number,TDim> > spData = (*iter).second;
1776 0 : const std::string& vtkName = (*iter).first;
1777 :
1778 : // write scalar value of this data
1779 : try{
1780 : write_nodal_data<TFunction,number>
1781 0 : (File, u, time, spData, 1, vtkName, grid, ssGrp, dim, numVert);
1782 : }
1783 0 : UG_CATCH_THROW("VTK::write_piece: Failed write nodal scalar Data.");
1784 : }
1785 :
1786 : // loop all vector data
1787 0 : for(VectorDataIterator iter = m_vVectorNodalData.begin();
1788 0 : iter != m_vVectorNodalData.end(); ++iter)
1789 : {
1790 : // get symb function
1791 : SmartPtr<UserData<MathVector<TDim>,TDim> > spData = (*iter).second;
1792 0 : const std::string& vtkName = (*iter).first;
1793 :
1794 : // write scalar value of this data
1795 : try{
1796 : write_nodal_data<TFunction,MathVector<TDim> >
1797 0 : (File, u, time, spData, 3, vtkName, grid, ssGrp, dim, numVert);
1798 : }
1799 0 : UG_CATCH_THROW("VTK::write_piece: Failed write nodal vector Data.");
1800 : }
1801 :
1802 : // loop all vector data
1803 0 : for(MatrixDataIterator iter = m_vMatrixNodalData.begin();
1804 0 : iter != m_vMatrixNodalData.end(); ++iter)
1805 : {
1806 : // get symb function
1807 : SmartPtr<UserData<MathMatrix<TDim,TDim>,TDim> > spData = (*iter).second;
1808 0 : const std::string& vtkName = (*iter).first;
1809 :
1810 : // write scalar value of this data
1811 : try{
1812 : write_nodal_data<TFunction,MathMatrix<TDim,TDim> >
1813 0 : (File, u, time, spData, 9, vtkName, grid, ssGrp, dim, numVert);
1814 : }
1815 0 : UG_CATCH_THROW("VTK::write_piece: Failed write nodal matrix Data.");
1816 : }
1817 :
1818 : // write closing tag
1819 0 : File << VTKFileWriter::normal;
1820 0 : File << " </PointData>\n";
1821 : }
1822 :
1823 :
1824 : ////////////////////////////////////////////////////////////////////////////////
1825 : // Cell Value
1826 : ////////////////////////////////////////////////////////////////////////////////
1827 :
1828 : template <int TDim>
1829 : template <typename TElem, typename TFunction, typename TData>
1830 0 : void VTKOutput<TDim>::
1831 : write_cell_data_elementwise(VTKFileWriter& File, TFunction& u, number time,
1832 : SmartPtr<UserData<TData, TDim> > spData,
1833 : Grid& grid, const int si)
1834 : {
1835 : // get reference element
1836 : typedef typename reference_element_traits<TElem>::reference_element_type
1837 : ref_elem_type;
1838 : static const int refDim = reference_element_traits<TElem>::dim;
1839 0 : static const ref_elem_type& refElem = Provider<ref_elem_type>::get();
1840 : static const size_t numCo = ref_elem_type::numCorners;
1841 :
1842 0 : if(m_bBinary)
1843 0 : File << VTKFileWriter::base64_binary;
1844 : else
1845 0 : File << VTKFileWriter::normal;
1846 :
1847 : // get iterators
1848 : typedef typename IteratorProvider<TFunction>::template traits<TElem>::const_iterator const_iterator;
1849 0 : const_iterator iterBegin = IteratorProvider<TFunction>::template begin<TElem>(u, si);
1850 0 : const_iterator iterEnd = IteratorProvider<TFunction>::template end<TElem>(u, si);
1851 :
1852 : MathVector<refDim> localIP;
1853 : MathVector<TDim> globIP;
1854 0 : AveragePositions(localIP, refElem.corners(), numCo);
1855 :
1856 0 : std::vector<MathVector<TDim> > vCorner(numCo);
1857 :
1858 : TData value;
1859 0 : bool bNeedFct = spData->requires_grid_fct();
1860 :
1861 : // loop all elements
1862 0 : for( ; iterBegin != iterEnd; ++iterBegin)
1863 : {
1864 : // get element
1865 : TElem *elem = *iterBegin;
1866 :
1867 : // update the reference mapping for the corners
1868 0 : CollectCornerCoordinates(vCorner, *elem, u.domain()->position_accessor(), true);
1869 :
1870 : // compute global integration points
1871 0 : AveragePositions(globIP, &vCorner[0], numCo);
1872 :
1873 : // get subset
1874 : int theSI = si;
1875 0 : if(si < 0) theSI = u.domain()->subset_handler()->get_subset_index(elem);
1876 :
1877 : // get local solution if needed
1878 0 : if(bNeedFct)
1879 : {
1880 : // create storage
1881 : LocalIndices ind;
1882 : LocalVector locU;
1883 :
1884 : // get global indices
1885 : u.indices(elem, ind);
1886 :
1887 : // adapt local algebra
1888 0 : locU.resize(ind);
1889 :
1890 : // read local values of u
1891 0 : GetLocalVector(locU, u);
1892 :
1893 : // compute data
1894 : try{
1895 : (*spData)(value, globIP, time, theSI, elem,
1896 : &vCorner[0], localIP, &locU);
1897 : }
1898 0 : UG_CATCH_THROW("VTK::write_cell_data_elementwise: Cannot evaluate data.");
1899 : }
1900 : else
1901 : {
1902 : // compute data
1903 : try{
1904 0 : (*spData)(value, globIP, time, theSI);
1905 : }
1906 0 : UG_CATCH_THROW("VTK::write_cell_data_elementwise: Cannot evaluate data.");
1907 : }
1908 :
1909 : // handle octahedral case separately by splitting into top and bottom pyramid with the same values
1910 : if(ref_elem_type::REFERENCE_OBJECT_ID != ROID_OCTAHEDRON)
1911 : {
1912 0 : write_item_to_file(File, value);
1913 : }
1914 : else
1915 : {
1916 0 : write_item_to_file(File, value); // top pyramid
1917 0 : write_item_to_file(File, value); // bottom pyramid
1918 : }
1919 : }
1920 :
1921 0 : }
1922 :
1923 :
1924 : template <int TDim>
1925 : template <typename TFunction, typename TData>
1926 0 : void VTKOutput<TDim>::
1927 : write_cell_data(VTKFileWriter& File, TFunction& u, number time,
1928 : SmartPtr<UserData<TData, TDim> > spData,
1929 : const int numCmp,
1930 : const std::string& name,
1931 : Grid& grid, const SubsetGroup& ssGrp, const int dim, const int numElem)
1932 : {
1933 0 : spData->set_function_pattern(u.function_pattern());
1934 :
1935 : // write opening tag
1936 0 : File << VTKFileWriter::normal;
1937 0 : File << " <DataArray type=\"Float32\" Name=\""<<name<<"\" "
1938 0 : "NumberOfComponents=\""<<numCmp<<"\" format="
1939 0 : << (m_bBinary ? "\"binary\"": "\"ascii\"") << ">\n";
1940 :
1941 0 : int n = sizeof(float) * numElem * numCmp;
1942 0 : if(m_bBinary)
1943 0 : File << VTKFileWriter::base64_binary << n;
1944 :
1945 : // switch dimension
1946 0 : for(size_t i = 0; i < ssGrp.size(); i++)
1947 0 : switch(dim)
1948 : {
1949 0 : case 1: write_cell_data_elementwise<RegularEdge,TFunction,TData>(File, u, time, spData, grid, ssGrp[i]);
1950 0 : write_cell_data_elementwise<ConstrainingEdge,TFunction,TData>(File, u, time, spData, grid, ssGrp[i]); break;
1951 0 : case 2: write_cell_data_elementwise<Triangle,TFunction,TData>(File, u, time, spData, grid, ssGrp[i]);
1952 0 : write_cell_data_elementwise<Quadrilateral,TFunction,TData>(File, u, time, spData, grid, ssGrp[i]);
1953 0 : write_cell_data_elementwise<ConstrainingTriangle,TFunction,TData>(File, u, time, spData, grid, ssGrp[i]);
1954 0 : write_cell_data_elementwise<ConstrainingQuadrilateral,TFunction,TData>(File, u, time, spData, grid, ssGrp[i]); break;
1955 0 : case 3: write_cell_data_elementwise<Tetrahedron,TFunction,TData>(File, u, time, spData, grid, ssGrp[i]);
1956 0 : write_cell_data_elementwise<Pyramid,TFunction,TData>(File, u, time, spData, grid, ssGrp[i]);
1957 0 : write_cell_data_elementwise<Prism,TFunction,TData>(File, u, time, spData, grid, ssGrp[i]);
1958 0 : write_cell_data_elementwise<Octahedron,TFunction,TData>(File, u, time, spData, grid, ssGrp[i]);
1959 0 : write_cell_data_elementwise<Hexahedron,TFunction,TData>(File, u, time, spData, grid, ssGrp[i]); break;
1960 0 : default: UG_THROW("VTK::write_cell_data: Dimension " << dim <<
1961 : " is not supported.");
1962 : }
1963 :
1964 : // write closing tag
1965 0 : File << VTKFileWriter::normal;
1966 0 : File << "\n </DataArray>\n";
1967 0 : };
1968 :
1969 :
1970 : template <int TDim>
1971 : template <typename TElem, typename TFunction>
1972 0 : void VTKOutput<TDim>::
1973 : write_cell_values_elementwise(VTKFileWriter& File, TFunction& u,
1974 : const std::vector<size_t>& vFct,
1975 : Grid& grid, const int si)
1976 : {
1977 : // get reference element
1978 : typedef typename reference_element_traits<TElem>::reference_element_type
1979 : ref_elem_type;
1980 :
1981 0 : static const ref_elem_type& refElem = Provider<ref_elem_type>::get();
1982 : static const ReferenceObjectID roid = ref_elem_type::REFERENCE_OBJECT_ID;
1983 : static const int dim = ref_elem_type::dim;
1984 : static const size_t numCo = ref_elem_type::numCorners;
1985 :
1986 : // index vector
1987 : std::vector<DoFIndex> vMultInd;
1988 :
1989 : // get iterators
1990 : typedef typename IteratorProvider<TFunction>::template traits<TElem>::const_iterator const_iterator;
1991 0 : const_iterator iterBegin = IteratorProvider<TFunction>::template begin<TElem>(u, si);
1992 0 : const_iterator iterEnd = IteratorProvider<TFunction>::template end<TElem>(u, si);
1993 :
1994 0 : if(m_bBinary)
1995 0 : File << VTKFileWriter::base64_binary;
1996 : else
1997 0 : File << VTKFileWriter::normal;
1998 :
1999 : MathVector<dim> localIP;
2000 0 : AveragePositions(localIP, refElem.corners(), numCo);
2001 :
2002 : // request for trial space
2003 : try{
2004 :
2005 : // avoid passing this code, since LocalShapeFunctionSet might not exist for RefElem-Type
2006 0 : if(iterBegin == iterEnd) return;
2007 :
2008 0 : std::vector<std::vector<number> > vvShape(vFct.size());
2009 0 : std::vector<size_t> vNsh(vFct.size());
2010 :
2011 0 : for(size_t f = 0; f < vFct.size(); ++f)
2012 : {
2013 0 : const LFEID lfeID = u.local_finite_element_id(vFct[f]);
2014 : const LocalShapeFunctionSet<dim>& lsfs
2015 : = LocalFiniteElementProvider::get<dim>(roid, lfeID);
2016 :
2017 0 : vNsh[f] = lsfs.num_sh();
2018 0 : lsfs.shapes(vvShape[f], localIP);
2019 : }
2020 :
2021 : // loop all elements
2022 0 : for( ; iterBegin != iterEnd; ++iterBegin)
2023 : {
2024 : // get element
2025 : TElem *elem = *iterBegin;
2026 :
2027 : // loop all components
2028 0 : for(size_t f = 0; f < vFct.size(); ++f)
2029 : {
2030 : // get multi index of vertex for the function
2031 0 : if(u.dof_indices(elem, vFct[f], vMultInd) != vNsh[f])
2032 0 : UG_THROW("VTK:write_cell_values_elementwise: "
2033 : "Number of shape functions for component "<<vFct[f]<<
2034 : " does not match number of DoFs");
2035 :
2036 : number ipVal = 0.0;
2037 0 : for(size_t sh = 0; sh < vNsh[f]; ++sh)
2038 : {
2039 0 : ipVal += DoFRef(u, vMultInd[sh]) * vvShape[f][sh];
2040 : }
2041 :
2042 : // flush stream
2043 : write_item_to_file(File, ipVal);
2044 :
2045 : // flush stream a second time in case of an octahedron as we handle it by splitting into a top and bottom pyramid
2046 : if(roid == ROID_OCTAHEDRON)
2047 : write_item_to_file(File, ipVal); // bottom pyramid
2048 : }
2049 :
2050 : // fill with zeros up to 3d if vector type
2051 0 : if(vFct.size() != 1){
2052 0 : for(size_t i = vFct.size(); i < 3; ++i) {
2053 0 : write_item_to_file(File, 0.f);
2054 :
2055 : // flush stream a second time in case of an octahedron as we handle it by splitting into a top and bottom pyramid
2056 : if(roid == ROID_OCTAHEDRON)
2057 0 : write_item_to_file(File, 0.f); // bottom pyramid
2058 : }
2059 : }
2060 : }
2061 0 : }
2062 0 : UG_CATCH_THROW("VTK: Could not find Shape function Set.");
2063 :
2064 : }
2065 :
2066 : template <int TDim>
2067 : template <typename TFunction>
2068 0 : void VTKOutput<TDim>::
2069 : write_cell_values(VTKFileWriter& File, TFunction& u,
2070 : const std::vector<size_t>& vFct, const std::string& name,
2071 : Grid& grid, const SubsetGroup& ssGrp, const int dim, const int numElem)
2072 : {
2073 : // write opening tag
2074 0 : File << VTKFileWriter::normal;
2075 0 : File << " <DataArray type=\"Float32\" Name=\""<<name<<"\" "
2076 0 : "NumberOfComponents=\""<<(vFct.size() == 1 ? 1 : 3)<<"\" format="
2077 0 : << (m_bBinary ? "\"binary\"": "\"ascii\"") << ">\n";
2078 :
2079 0 : int n = sizeof(float) * numElem * (vFct.size() == 1 ? 1 : 3);
2080 0 : if(m_bBinary)
2081 0 : File << VTKFileWriter::base64_binary << n;
2082 :
2083 : // switch dimension
2084 0 : for(size_t i = 0; i < ssGrp.size(); i++)
2085 0 : switch(dim)
2086 : {
2087 0 : case 1: write_cell_values_elementwise<RegularEdge>(File, u, vFct, grid, ssGrp[i]);
2088 0 : write_cell_values_elementwise<ConstrainingEdge>(File, u, vFct, grid, ssGrp[i]); break;
2089 0 : case 2: write_cell_values_elementwise<Triangle>(File, u, vFct, grid, ssGrp[i]);
2090 0 : write_cell_values_elementwise<Quadrilateral>(File, u, vFct, grid, ssGrp[i]);
2091 0 : write_cell_values_elementwise<ConstrainingTriangle>(File, u, vFct, grid, ssGrp[i]);
2092 0 : write_cell_values_elementwise<ConstrainingQuadrilateral>(File, u, vFct, grid, ssGrp[i]); break;
2093 0 : case 3: write_cell_values_elementwise<Tetrahedron>(File, u, vFct, grid, ssGrp[i]);
2094 0 : write_cell_values_elementwise<Pyramid>(File, u, vFct, grid, ssGrp[i]);
2095 0 : write_cell_values_elementwise<Prism>(File, u, vFct, grid, ssGrp[i]);
2096 0 : write_cell_values_elementwise<Octahedron>(File, u, vFct, grid, ssGrp[i]);
2097 0 : write_cell_values_elementwise<Hexahedron>(File, u, vFct, grid, ssGrp[i]); break;
2098 0 : default: UG_THROW("VTK::write_cell_values: Dimension " << dim <<
2099 : " is not supported.");
2100 : }
2101 :
2102 : // write closing tag
2103 0 : File << VTKFileWriter::normal;
2104 0 : File << "\n </DataArray>\n";
2105 0 : };
2106 :
2107 : template <int TDim>
2108 : template <typename TFunction>
2109 0 : void VTKOutput<TDim>::
2110 : write_cell_values_piece(VTKFileWriter& File, TFunction& u, number time, Grid& grid,
2111 : const SubsetGroup& ssGrp, const int dim, const int numElem)
2112 : {
2113 0 : if(!m_vSymbFct.empty()){
2114 0 : for(std::map<std::string, std::vector<std::string> >::const_iterator iter =
2115 0 : m_vSymbFct.begin(); iter != m_vSymbFct.end(); ++iter){
2116 0 : const std::vector<std::string>& symbNames = (*iter).second;
2117 0 : const std::string& vtkName = (*iter).first;
2118 :
2119 : bool bContinuous = true;
2120 0 : for(size_t i = 0; i < symbNames.size(); ++i){
2121 : size_t fct = u.fct_id_by_name(symbNames[i].c_str());
2122 0 : if(!LocalFiniteElementProvider::continuous(u.local_finite_element_id(fct))){
2123 : bContinuous = false; break;
2124 : }
2125 : }
2126 :
2127 0 : if(!bContinuous){
2128 0 : m_vSymbFctElem[vtkName] = symbNames;
2129 : }
2130 : }
2131 : }
2132 :
2133 : // check if something to do
2134 0 : if(m_vSymbFctElem.empty() && m_vScalarElemData.empty() && m_vVectorElemData.empty())
2135 : return;
2136 :
2137 : // write opening tag to indicate point data
2138 0 : File << VTKFileWriter::normal;
2139 0 : File << " <CellData>\n";
2140 :
2141 : // loop all selected symbolic names
2142 0 : for(ComponentsIterator iter = m_vSymbFctElem.begin();
2143 0 : iter != m_vSymbFctElem.end(); ++iter)
2144 : {
2145 : // get symb function
2146 : const std::vector<std::string>& symbNames = (*iter).second;
2147 0 : const std::string& vtkName = (*iter).first;
2148 :
2149 : // create function group
2150 0 : std::vector<size_t> fctGrp(symbNames.size());
2151 0 : for(size_t i = 0; i < symbNames.size(); ++i)
2152 0 : fctGrp[i] = u.fct_id_by_name(symbNames[i].c_str());
2153 :
2154 : // check that all functions are contained in subset
2155 : bool bContained = true;
2156 0 : for(size_t i = 0; i < fctGrp.size(); ++i)
2157 0 : for(size_t j = 0; j < ssGrp.size(); ++j)
2158 0 : if(!u.is_def_in_subset(fctGrp[i], ssGrp[j]))
2159 : bContained = false;
2160 :
2161 0 : if(!bContained) continue;
2162 :
2163 : // write scalar value of this function
2164 : try{
2165 0 : write_cell_values(File, u, fctGrp, vtkName, grid, ssGrp, dim, numElem);
2166 : }
2167 0 : UG_CATCH_THROW("VTK::write_piece: Failed write cell Values.");
2168 : }
2169 :
2170 : // loop all scalar data
2171 0 : for(ScalarDataIterator iter = m_vScalarElemData.begin();
2172 0 : iter != m_vScalarElemData.end(); ++iter)
2173 : {
2174 : // get symb function
2175 : SmartPtr<UserData<number,TDim> > spData = (*iter).second;
2176 0 : const std::string& vtkName = (*iter).first;
2177 :
2178 : // write scalar value of this data
2179 : try{
2180 : write_cell_data<TFunction,number>
2181 0 : (File, u, time, spData, 1, vtkName, grid, ssGrp, dim, numElem);
2182 : }
2183 0 : UG_CATCH_THROW("VTK::write_piece: Failed to write cell scalar Data.");
2184 : }
2185 :
2186 : // loop all vector data
2187 0 : for(VectorDataIterator iter = m_vVectorElemData.begin();
2188 0 : iter != m_vVectorElemData.end(); ++iter)
2189 : {
2190 : // get symb function
2191 : SmartPtr<UserData<MathVector<TDim>,TDim> > spData = (*iter).second;
2192 0 : const std::string& vtkName = (*iter).first;
2193 :
2194 : // write scalar value of this data
2195 : try{
2196 : write_cell_data<TFunction,MathVector<TDim> >
2197 0 : (File, u, time, spData, 3, vtkName, grid, ssGrp, dim, numElem);
2198 : }
2199 0 : UG_CATCH_THROW("VTK::write_piece: Failed to write cell vector Data.");
2200 : }
2201 :
2202 : // loop all vector data
2203 0 : for(MatrixDataIterator iter = m_vMatrixElemData.begin();
2204 0 : iter != m_vMatrixElemData.end(); ++iter)
2205 : {
2206 : // get symb function
2207 : SmartPtr<UserData<MathMatrix<TDim,TDim>,TDim> > spData = (*iter).second;
2208 0 : const std::string& vtkName = (*iter).first;
2209 :
2210 : // write scalar value of this data
2211 : try{
2212 : write_cell_data<TFunction,MathMatrix<TDim,TDim> >
2213 0 : (File, u, time, spData, 9, vtkName, grid, ssGrp, dim, numElem);
2214 : }
2215 0 : UG_CATCH_THROW("VTK::write_piece: Failed to write cell matrix Data.");
2216 : }
2217 :
2218 : // write closing tag
2219 0 : File << VTKFileWriter::normal;
2220 0 : File << " </CellData>\n";
2221 : }
2222 :
2223 : ////////////////////////////////////////////////////////////////////////////////
2224 : // Grouping Files
2225 : ////////////////////////////////////////////////////////////////////////////////
2226 :
2227 : template <int TDim>
2228 : template <typename TFunction>
2229 : void VTKOutput<TDim>::
2230 : write_pvtu(TFunction& u, const std::string& filename,
2231 : int si, int step, number time)
2232 : {
2233 : #ifdef UG_PARALLEL
2234 : // File pointer
2235 : FILE* file;
2236 :
2237 : // file name
2238 : std::string name;
2239 :
2240 : // get and check number of procs (only for numProcs > 1 we write the pvtu)
2241 : int numProcs = pcl::NumProcs();
2242 : if(numProcs == 1) return;
2243 :
2244 : // check if this proc is output proc
2245 : bool isOutputProc = GetLogAssistant().is_output_process();
2246 :
2247 : // max subset
2248 : int maxSi = u.num_subsets() - 1;
2249 :
2250 : // only the master process writes this file
2251 : if (isOutputProc)
2252 : {
2253 : // get name for *.pvtu file
2254 : pvtu_filename(name, filename, si, maxSi, step);
2255 :
2256 : // open file
2257 : file = fopen(name.c_str(), "w");
2258 : if (file == NULL)
2259 : UG_THROW("VTKOutput: Cannot print to file.");
2260 :
2261 : // Write to file
2262 : fprintf(file, "<?xml version=\"1.0\"?>\n");
2263 :
2264 : // Write comment
2265 : write_comment_printf(file);
2266 :
2267 : fprintf(file, "<VTKFile type=\"PUnstructuredGrid\" version=\"0.1\">\n");
2268 : fprintf(file, " <Time timestep=\"%.17g\"/>\n", time);
2269 : fprintf(file, " <PUnstructuredGrid GhostLevel=\"0\">\n");
2270 : fprintf(file, " <PPoints>\n");
2271 : fprintf(file, " <PDataArray type=\"Float32\" NumberOfComponents=\"3\"/>\n");
2272 : fprintf(file, " </PPoints>\n");
2273 :
2274 : // Node Data
2275 : if(!m_vSymbFctNodal.empty() || !m_vScalarNodalData.empty() || !m_vVectorNodalData.empty())
2276 : {
2277 : fprintf(file, " <PPointData>\n");
2278 : // for(size_t sym = 0; sym < m_vSymbFctNodal.size(); ++sym)
2279 : for(ComponentsIterator iter = m_vSymbFctNodal.begin();
2280 : iter != m_vSymbFctNodal.end(); ++iter)
2281 : {
2282 : // get symb function
2283 : const std::vector<std::string>& symbNames = (*iter).second;
2284 : const std::string& vtkName = (*iter).first;
2285 :
2286 : // create function group
2287 : std::vector<size_t> fctGrp(symbNames.size());
2288 : for(size_t i = 0; i < symbNames.size(); ++i)
2289 : fctGrp[i] = u.fct_id_by_name(symbNames[i].c_str());
2290 :
2291 : // check that all functions are contained in subset
2292 : bool bContained = true;
2293 : for(size_t i = 0; i < fctGrp.size(); ++i)
2294 : if(!u.is_def_in_subset(fctGrp[i], si))
2295 : bContained = false;
2296 :
2297 : if(!bContained) continue;
2298 :
2299 : fprintf(file, " <PDataArray type=\"Float32\" Name=\"%s\" "
2300 : "NumberOfComponents=\"%d\"/>\n",
2301 : vtkName.c_str(), (fctGrp.size() == 1 ? 1 : 3));
2302 : }
2303 :
2304 : // loop all scalar data
2305 : for(ScalarDataIterator iter = m_vScalarNodalData.begin();
2306 : iter != m_vScalarNodalData.end(); ++iter)
2307 : {
2308 : // get symb function
2309 : const std::string& vtkName = (*iter).first;
2310 :
2311 : fprintf(file, " <PDataArray type=\"Float32\" Name=\"%s\" "
2312 : "NumberOfComponents=\"%d\"/>\n",
2313 : vtkName.c_str(), 1);
2314 : }
2315 :
2316 : // loop all vector data
2317 : for(VectorDataIterator iter = m_vVectorNodalData.begin();
2318 : iter != m_vVectorNodalData.end(); ++iter)
2319 : {
2320 : // get symb function
2321 : const std::string& vtkName = (*iter).first;
2322 :
2323 : fprintf(file, " <PDataArray type=\"Float32\" Name=\"%s\" "
2324 : "NumberOfComponents=\"%d\"/>\n",
2325 : vtkName.c_str(), 3);
2326 : }
2327 :
2328 : // loop all matrix data
2329 : for(MatrixDataIterator iter = m_vMatrixNodalData.begin();
2330 : iter != m_vMatrixNodalData.end(); ++iter)
2331 : {
2332 : // get symb function
2333 : const std::string& vtkName = (*iter).first;
2334 :
2335 : fprintf(file, " <PDataArray type=\"Float32\" Name=\"%s\" "
2336 : "NumberOfComponents=\"%d\"/>\n",
2337 : vtkName.c_str(), 9);
2338 : }
2339 : fprintf(file, " </PPointData>\n");
2340 : }
2341 :
2342 : // Elem Data
2343 : if(!m_vSymbFctElem.empty() || !m_vScalarElemData.empty() || !m_vVectorElemData.empty() || m_bWriteSubsetIndices || m_bWriteProcRanks) //TODO: cleanup
2344 : {
2345 : fprintf(file, " <PCellData>\n");
2346 : // for(size_t sym = 0; sym < m_vSymbFctElem.size(); ++sym)
2347 : for(ComponentsIterator iter = m_vSymbFctElem.begin();
2348 : iter != m_vSymbFctElem.end(); ++iter)
2349 : {
2350 : // get symb function
2351 : const std::vector<std::string>& symbNames = (*iter).second;//m_vSymbFctElem[sym].first;
2352 : const std::string& vtkName = (*iter).first;//m_vSymbFctElem[sym].second;
2353 :
2354 : // create function group
2355 : std::vector<size_t> fctGrp(symbNames.size());
2356 : for(size_t i = 0; i < symbNames.size(); ++i)
2357 : fctGrp[i] = u.fct_id_by_name(symbNames[i].c_str());
2358 :
2359 : // check that all functions are contained in subset
2360 : bool bContained = true;
2361 : for(size_t i = 0; i < fctGrp.size(); ++i)
2362 : if(!u.is_def_in_subset(fctGrp[i], si))
2363 : bContained = false;
2364 :
2365 : if(!bContained) continue;
2366 :
2367 : fprintf(file, " <PDataArray type=\"Float32\" Name=\"%s\" "
2368 : "NumberOfComponents=\"%d\"/>\n",
2369 : vtkName.c_str(), (fctGrp.size() == 1 ? 1 : 3));
2370 : }
2371 :
2372 : // loop all scalar data
2373 : for(ScalarDataIterator iter = m_vScalarElemData.begin();
2374 : iter != m_vScalarElemData.end(); ++iter)
2375 : {
2376 : // get symb function
2377 : const std::string& vtkName = (*iter).first;
2378 :
2379 : fprintf(file, " <PDataArray type=\"Float32\" Name=\"%s\" "
2380 : "NumberOfComponents=\"%d\"/>\n",
2381 : vtkName.c_str(), 1);
2382 : }
2383 :
2384 : //TODO: cleanup!!
2385 : if(m_bWriteSubsetIndices){
2386 : fprintf(file, " <DataArray type=\"Int8\" Name=\"regions\" "
2387 : "NumberOfComponents=\"1\"/>\n");
2388 : }
2389 : if(m_bWriteProcRanks){
2390 : fprintf(file, " <DataArray type=\"Int8\" Name=\"proc_ranks\" "
2391 : "NumberOfComponents=\"1\"/>\n");
2392 : }
2393 :
2394 : // loop all vector data
2395 : for(VectorDataIterator iter = m_vVectorElemData.begin();
2396 : iter != m_vVectorElemData.end(); ++iter)
2397 : {
2398 : // get symb function
2399 : const std::string& vtkName = (*iter).first;
2400 :
2401 : fprintf(file, " <PDataArray type=\"Float32\" Name=\"%s\" "
2402 : "NumberOfComponents=\"%d\"/>\n",
2403 : vtkName.c_str(), 3);
2404 : }
2405 :
2406 : // loop all vector data
2407 : for(MatrixDataIterator iter = m_vMatrixElemData.begin();
2408 : iter != m_vMatrixElemData.end(); ++iter)
2409 : {
2410 : // get symb function
2411 : const std::string& vtkName = (*iter).first;
2412 :
2413 : fprintf(file, " <PDataArray type=\"Float32\" Name=\"%s\" "
2414 : "NumberOfComponents=\"%d\"/>\n",
2415 : vtkName.c_str(), 9);
2416 : }
2417 : fprintf(file, " </PCellData>\n");
2418 : }
2419 :
2420 : // include files from all procs
2421 : for (int i = 0; i < numProcs; i++) {
2422 : vtu_filename(name, filename, i, si, maxSi, step);
2423 : name = FilenameWithoutPath(name);
2424 : fprintf(file, " <Piece Source=\"%s\"/>\n", name.c_str());
2425 : }
2426 :
2427 : // write closing tags
2428 : fprintf(file, " </PUnstructuredGrid>\n");
2429 : fprintf(file, "</VTKFile>\n");
2430 : fclose(file);
2431 : }
2432 :
2433 : PCL_DEBUG_BARRIER_ALL();
2434 : #endif
2435 : }
2436 :
2437 :
2438 : template <int TDim>
2439 : template <typename TFunction>
2440 0 : void VTKOutput<TDim>::
2441 : write_time_pvd(const char* filename, TFunction& u)
2442 : {
2443 : // File
2444 : FILE* file;
2445 :
2446 : // filename
2447 : std::string name;
2448 :
2449 : // get some numbers
2450 0 : bool isOutputProc = GetLogAssistant().is_output_process();
2451 : int numProcs = 1;
2452 : #ifdef UG_PARALLEL
2453 : numProcs = pcl::NumProcs();
2454 : #endif
2455 :
2456 : // get time steps
2457 0 : std::vector<number>& vTimestep = m_mTimestep[filename];
2458 :
2459 : // change locale to ensure decimal . is really a .
2460 0 : char* oldLocale = setlocale (LC_ALL, NULL);
2461 0 : setlocale(LC_NUMERIC, "C");
2462 :
2463 0 : if (isOutputProc)
2464 : {
2465 : // get file name
2466 0 : pvd_filename(name, filename);
2467 :
2468 : // open file
2469 0 : file = fopen(name.c_str(), "w");
2470 0 : if (file == NULL)
2471 0 : UG_THROW("Cannot print to file.");
2472 :
2473 : // Write to file
2474 : fprintf(file, "<?xml version=\"1.0\"?>\n");
2475 :
2476 : // Write comment
2477 0 : write_comment_printf(file);
2478 :
2479 : fprintf(file, "<VTKFile type=\"Collection\" version=\"0.1\">\n");
2480 : fprintf(file, " <Collection>\n");
2481 :
2482 : // check functions
2483 : bool bEverywhere = true;
2484 0 : for(size_t fct = 0; fct < u.num_fct(); ++fct)
2485 : {
2486 : // check if function is defined everywhere
2487 0 : if(!u.is_def_everywhere(fct))
2488 : bEverywhere = false;
2489 : }
2490 :
2491 : // include files from all procs
2492 0 : if(bEverywhere)
2493 : {
2494 0 : for(int step = 0; step < (int)vTimestep.size(); ++step)
2495 : {
2496 0 : vtu_filename(name, filename, 0, -1, 0, step);
2497 : if(numProcs > 1) pvtu_filename(name, filename, -1, 0, step);
2498 :
2499 0 : name = FilenameWithoutPath(name);
2500 0 : fprintf(file, " <DataSet timestep=\"%.17g\" part=\"%d\" file=\"%s\"/>\n",
2501 : vTimestep[step], 0, name.c_str());
2502 : }
2503 : }
2504 : else
2505 : {
2506 0 : for(int step = 0; step < (int)vTimestep.size(); ++step)
2507 0 : for(int si = 0; si < u.num_subsets(); ++si)
2508 : {
2509 0 : vtu_filename(name, filename, 0, si, u.num_subsets()-1, step);
2510 : if(numProcs > 1) pvtu_filename(name, filename, si, u.num_subsets()-1, step);
2511 :
2512 0 : name = FilenameWithoutPath(name);
2513 0 : fprintf(file, " <DataSet timestep=\"%.17g\" part=\"%d\" file=\"%s\"/>\n",
2514 : vTimestep[step], si, name.c_str());
2515 : }
2516 : }
2517 :
2518 : // close file
2519 : fprintf(file, " </Collection>\n");
2520 : fprintf(file, "</VTKFile>\n");
2521 0 : fclose(file);
2522 : }
2523 :
2524 : // restore old locale
2525 0 : setlocale(LC_NUMERIC, oldLocale);
2526 :
2527 : #ifdef UG_PARALLEL
2528 : PCL_DEBUG_BARRIER_ALL();
2529 : #endif
2530 0 : }
2531 :
2532 :
2533 : template <int TDim>
2534 : template <typename TFunction>
2535 0 : void VTKOutput<TDim>::
2536 : write_time_processwise_pvd(const char* filename, TFunction& u)
2537 : {
2538 : // File
2539 : FILE* file;
2540 :
2541 : // filename
2542 : std::string name;
2543 :
2544 : // get some numbers
2545 0 : bool isOutputProc = GetLogAssistant().is_output_process();
2546 : int numProcs = 1;
2547 : #ifdef UG_PARALLEL
2548 : numProcs = pcl::NumProcs();
2549 : #endif
2550 :
2551 : // get time steps
2552 0 : std::vector<number>& vTimestep = m_mTimestep[filename];
2553 :
2554 : // change locale to ensure decimal . is really a .
2555 0 : char* oldLocale = setlocale (LC_ALL, NULL);
2556 0 : setlocale(LC_NUMERIC, "C");
2557 :
2558 : if (isOutputProc && numProcs > 1)
2559 : {
2560 : // adjust filename
2561 : std::string procName = filename;
2562 : procName.append("_processwise");
2563 : pvd_filename(name, procName);
2564 :
2565 : // open file
2566 : file = fopen(name.c_str(), "w");
2567 : if (file == NULL)
2568 : UG_THROW("Cannot print to file.");
2569 :
2570 : // Write to file
2571 : fprintf(file, "<?xml version=\"1.0\"?>\n");
2572 :
2573 : // Write comment
2574 : write_comment_printf(file);
2575 :
2576 : fprintf(file, "<VTKFile type=\"Collection\" version=\"0.1\">\n");
2577 : fprintf(file, " <Collection>\n");
2578 :
2579 : // include files from all procs
2580 : for(int step = 0; step < (int)vTimestep.size(); ++step)
2581 : for (int rank = 0; rank < numProcs; rank++)
2582 : for(int si = 0; si < u.num_subsets(); ++si)
2583 : {
2584 : vtu_filename(name, filename, rank, si, u.num_subsets()-1, step);
2585 : if(numProcs > 1) pvtu_filename(name, filename, si, u.num_subsets()-1, step);
2586 :
2587 : name = FilenameWithoutPath(name);
2588 : fprintf(file, " <DataSet timestep=\"%.17g\" part=\"%d\" file=\"%s\"/>\n",
2589 : vTimestep[step], rank, name.c_str());
2590 : }
2591 :
2592 : // end file
2593 : fprintf(file, " </Collection>\n");
2594 : fprintf(file, "</VTKFile>\n");
2595 : fclose(file);
2596 : }
2597 :
2598 : // restore old locale
2599 0 : setlocale(LC_NUMERIC, oldLocale);
2600 :
2601 : #ifdef UG_PARALLEL
2602 : PCL_DEBUG_BARRIER_ALL();
2603 : #endif
2604 0 : }
2605 :
2606 :
2607 :
2608 : template <int TDim>
2609 : template <typename TFunction>
2610 0 : void VTKOutput<TDim>::
2611 : write_time_pvd_subset(const char* filename, TFunction& u, int si)
2612 : {
2613 : // File
2614 : FILE* file;
2615 :
2616 : // filename
2617 : std::string name;
2618 :
2619 : // get some numbers
2620 0 : bool isOutputProc = GetLogAssistant().is_output_process();
2621 : int numProcs = 1;
2622 : #ifdef UG_PARALLEL
2623 : numProcs = pcl::NumProcs();
2624 : #endif
2625 :
2626 : // get time steps
2627 0 : std::vector<number>& vTimestep = m_mTimestep[filename];
2628 :
2629 : // change locale to ensure decimal . is really a .
2630 0 : char* oldLocale = setlocale (LC_ALL, NULL);
2631 0 : setlocale(LC_NUMERIC, "C");
2632 :
2633 0 : if (isOutputProc)
2634 : {
2635 : // get file name
2636 0 : pvd_filename(name, filename);
2637 :
2638 : // open file
2639 0 : file = fopen(name.c_str(), "w");
2640 0 : if (file == NULL)
2641 0 : UG_THROW("Cannot print to file.");
2642 :
2643 : // Write to file
2644 : fprintf(file, "<?xml version=\"1.0\"?>\n");
2645 :
2646 : // Write comment
2647 0 : write_comment_printf(file);
2648 :
2649 : fprintf(file, "<VTKFile type=\"Collection\" version=\"0.1\">\n");
2650 : fprintf(file, " <Collection>\n");
2651 :
2652 : // include files from all procs
2653 0 : for(int step = 0; step < (int)vTimestep.size(); ++step)
2654 : {
2655 0 : vtu_filename(name, filename, 0, si, u.num_subsets()-1, step);
2656 : if(numProcs > 1) pvtu_filename(name, filename, si, u.num_subsets()-1, step);
2657 :
2658 0 : name = FilenameWithoutPath(name);
2659 0 : fprintf(file, " <DataSet timestep=\"%g\" part=\"%d\" file=\"%s\"/>\n",
2660 : vTimestep[step], si, name.c_str());
2661 : }
2662 :
2663 : // close file
2664 : fprintf(file, " </Collection>\n");
2665 : fprintf(file, "</VTKFile>\n");
2666 0 : fclose(file);
2667 : }
2668 :
2669 : // restore old locale
2670 0 : setlocale(LC_NUMERIC, oldLocale);
2671 :
2672 : #ifdef UG_PARALLEL
2673 : PCL_DEBUG_BARRIER_ALL();
2674 : #endif
2675 0 : }
2676 :
2677 : }
2678 :
2679 :
|
