Line data Source code
1 : /*
2 : * Copyright (c) 2010-2015: G-CSC, Goethe University Frankfurt
3 : * Author: Markus Breit
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 : #ifndef __H__UG__LIB_DISC__IO__VTK_EXPORT_HO__
34 : #define __H__UG__LIB_DISC__IO__VTK_EXPORT_HO__
35 :
36 : #include <cmath> // for ceil, log2
37 : #include <cstddef> // for size_t
38 : #include <limits> // for numeric_limits
39 : #include <ostream> // for string, operator<<, basic_ostream, endl
40 : #include <string> // for char_traits, allocator, operator+, basic_string
41 : #include <vector> // for vector
42 :
43 : #include "common/assert.h" // for UG_ASSERT
44 : #include "common/error.h" // for UG_CATCH_THROW, UG_COND_THROW
45 : #include "common/types.h" // for number
46 : #include "common/math/math_vector_matrix/math_vector.h" // for MathVector
47 : #include "common/util/smart_pointer.h" // for SmartPtr, ConstSmartPtr, make_sp
48 : #include "lib_disc/common/function_group.h" // for FunctionGroup
49 : #include "lib_disc/common/multi_index.h" // for DoFIndex, DoFRef
50 : #include "lib_disc/dof_manager/dof_distribution.h" // for DoFDistribution, DoFDistribution::traits
51 : #include "lib_disc/function_spaces/dof_position_util.h" // for InnerDoFPosition
52 : #include "lib_disc/function_spaces/grid_function_global_user_data.h" // for GlobalGridFunctionNumberData
53 : #include "lib_disc/io/vtkoutput.h" // for VTKOutput
54 : #include "lib_disc/local_finite_element/local_finite_element_id.h" // for LFEID
55 : #include "lib_grid/algorithms/debug_util.h" // for ElementDebugInfo
56 : #include "lib_grid/algorithms/selection_util.h" // for SelectAssociatedGridObjects
57 : #include "lib_grid/attachments/attachment_pipe.h" // for AttachmentAccessor
58 : #include "lib_grid/common_attachments.h" // for AVertex
59 : #include "lib_grid/grid/grid.h" // for Grid::VertexAttachmentAccessor::VertexAttachmentAccessor<TAt...
60 : #include "lib_grid/grid/grid_base_object_traits.h" // for VertexIterator
61 : #include "lib_grid/grid/grid_base_objects.h" // for CustomVertexGroup, Vertex, Edge (ptr only), Face (ptr only)
62 : #ifdef UG_PARALLEL
63 : #include "lib_grid/parallelization/parallel_refinement/parallel_refinement.h" // for ParallelGlobalRefiner_MultiGrid
64 : #else
65 : #include "lib_grid/refinement/global_multi_grid_refiner.h" // for GlobalMultiGridRefiner
66 : #endif
67 : #include "lib_grid/multi_grid.h" // for MultiGrid
68 : #include "lib_grid/tools/selector_grid.h" // for Selector
69 : #include "lib_grid/tools/subset_group.h" // for SubsetGroup
70 :
71 :
72 : namespace ug {
73 :
74 :
75 : template <typename TDomain, class TElem>
76 0 : inline void CopySelectedElements
77 : (
78 : SmartPtr<TDomain> destDom,
79 : SmartPtr<TDomain> srcDom,
80 : Selector& sel,
81 : AVertex& aNewVrt
82 : )
83 : {
84 : typedef typename TDomain::subset_handler_type SH_type;
85 :
86 0 : MultiGrid& srcGrid = *srcDom->grid();
87 0 : MultiGrid& destGrid = *destDom->grid();
88 :
89 0 : ConstSmartPtr<SH_type> srcSH = srcDom->subset_handler();
90 : SmartPtr<SH_type> destSH = destDom->subset_handler();
91 :
92 : Grid::VertexAttachmentAccessor<AVertex> aaNewVrt(srcGrid, aNewVrt);
93 :
94 : CustomVertexGroup vrts;
95 : typedef typename Grid::traits<TElem>::iterator iter_t;
96 :
97 0 : for (iter_t eiter = sel.begin<TElem>(); eiter != sel.end<TElem>(); ++eiter)
98 : {
99 : TElem* e = *eiter;
100 0 : vrts.resize(e->num_vertices());
101 0 : for (size_t iv = 0; iv < e->num_vertices(); ++iv)
102 0 : vrts.set_vertex(iv, aaNewVrt[e->vertex(iv)]);
103 :
104 : TElem* ne;
105 0 : try {ne = *destGrid.create_by_cloning(e, vrts);}
106 0 : UG_CATCH_THROW("New element could not be created.");
107 0 : destSH->assign_subset(ne, srcSH->get_subset_index(e));
108 : }
109 0 : }
110 :
111 : template <typename TDomain>
112 0 : inline void CopySelected
113 : (
114 : SmartPtr<TDomain> destDom,
115 : SmartPtr<TDomain> srcDom,
116 : Selector& sel
117 : )
118 : {
119 : typedef typename TDomain::position_attachment_type APos_type;
120 : typedef typename TDomain::subset_handler_type SH_type;
121 :
122 0 : MultiGrid& srcGrid = *srcDom->grid();
123 0 : MultiGrid& destGrid = *destDom->grid();
124 :
125 0 : ConstSmartPtr<SH_type> srcSH = srcDom->subset_handler();
126 : SmartPtr<SH_type> destSH = destDom->subset_handler();
127 :
128 : APos_type& aPos = srcDom->position_attachment();
129 0 : UG_COND_THROW(!srcGrid.has_vertex_attachment(aPos), "Position attachment required.");
130 : Grid::VertexAttachmentAccessor<APos_type> aaPosSrc(srcGrid, aPos);
131 0 : if (!destGrid.has_vertex_attachment(aPos))
132 0 : destGrid.attach_to_vertices(aPos);
133 : Grid::VertexAttachmentAccessor<APos_type> aaPosDest(destGrid, aPos);
134 :
135 : AVertex aNewVrt;
136 0 : srcGrid.attach_to_vertices(aNewVrt);
137 : Grid::VertexAttachmentAccessor<AVertex> aaNewVrt(srcGrid, aNewVrt);
138 :
139 0 : for (int si = destSH->num_subsets(); si < srcSH->num_subsets(); ++si)
140 0 : destSH->subset_info(si) = srcSH->subset_info(si);
141 :
142 0 : SelectAssociatedGridObjects(sel);
143 :
144 : // create new vertices in destGrid
145 0 : for (VertexIterator viter = sel.begin<Vertex>(); viter != sel.end<Vertex>(); ++viter)
146 : {
147 : Vertex* v = *viter;
148 : Vertex* nv;
149 0 : try {nv = *destGrid.create_by_cloning(v);}
150 0 : UG_CATCH_THROW("New vertex could not be created.");
151 0 : aaNewVrt[v] = nv;
152 : aaPosDest[nv] = aaPosSrc[v];
153 0 : destSH->assign_subset(nv, srcSH->get_subset_index(v));
154 : }
155 :
156 0 : CopySelectedElements<TDomain,Edge>(destDom, srcDom, sel, aNewVrt);
157 0 : CopySelectedElements<TDomain,Face>(destDom, srcDom, sel, aNewVrt);
158 0 : CopySelectedElements<TDomain,Volume>(destDom, srcDom, sel, aNewVrt);
159 :
160 : srcGrid.detach_from_vertices(aNewVrt);
161 0 : }
162 :
163 :
164 : template <typename TElem, typename TGridFunction, typename TGGFND>
165 0 : inline void interpolate_from_original_fct
166 : (
167 : SmartPtr<TGridFunction> u_new,
168 : const TGGFND& u_orig,
169 : size_t fct,
170 : const LFEID& lfeid
171 : )
172 : {
173 : typedef typename TGridFunction::domain_type dom_type;
174 : static const int dim = dom_type::dim;
175 : typedef typename DoFDistribution::traits<TElem>::const_iterator const_iter_type;
176 :
177 : // interpolate subset-wise
178 0 : size_t numSubsets = u_new->num_subsets();
179 0 : for (size_t si = 0; si < numSubsets; ++si)
180 : {
181 0 : if (!u_new->is_def_in_subset(fct, si)) continue;
182 :
183 : const_iter_type elem_iter = u_new->template begin<TElem>(si);
184 : const_iter_type iterEnd = u_new->template end<TElem>(si);
185 :
186 : std::vector<DoFIndex> ind;
187 0 : for (; elem_iter != iterEnd; ++elem_iter)
188 : {
189 : u_new->inner_dof_indices(*elem_iter, fct, ind);
190 :
191 : // get dof positions
192 : std::vector<MathVector<dim> > globPos;
193 0 : InnerDoFPosition<dom_type>(globPos, *elem_iter, *u_new->domain(), lfeid);
194 :
195 : UG_ASSERT(globPos.size() == ind.size(),
196 : "#DoF mismatch: InnerDoFPosition() found " << globPos.size()
197 : << ", but grid function has " << ind.size() << std::endl
198 : << "on " << ElementDebugInfo(*u_new->domain()->grid(), *elem_iter) << ".");
199 :
200 : // write values in new grid function
201 0 : for (size_t dof = 0; dof < ind.size(); ++dof)
202 : {
203 0 : if (!u_orig.evaluate(DoFRef(*u_new, ind[dof]), globPos[dof]))
204 : {
205 0 : DoFRef(*u_new, ind[dof]) = std::numeric_limits<number>::quiet_NaN();
206 : //UG_THROW("Interpolation onto new grid did not succeed.\n"
207 : // "DoF with coords " << globPos[dof] << " is out of range.");
208 : }
209 : }
210 : }
211 : }
212 0 : }
213 :
214 :
215 : template <typename TGridFunction, int trueDim>
216 0 : void vtk_export_ho
217 : (
218 : SmartPtr<TGridFunction> u,
219 : const std::vector<std::string>& vFct,
220 : size_t order,
221 : SmartPtr<VTKOutput<TGridFunction::domain_type::dim> > vtkOutput,
222 : const char* filename,
223 : size_t step,
224 : number time,
225 : const SubsetGroup& ssg
226 : )
227 : {
228 : // for order 1, use the given grid function as is
229 0 : if (order == 1)
230 : {
231 : // print out new grid function on desired subsets
232 0 : bool printAll = ssg.size() == (size_t) u->subset_handler()->num_subsets();
233 0 : if (printAll)
234 0 : vtkOutput->print(filename, *u, step, time);
235 : else
236 : {
237 0 : for (size_t s = 0; s < ssg.size(); ++s)
238 0 : vtkOutput->print_subset(filename, *u, ssg[s], step, time);
239 : }
240 :
241 0 : return;
242 : }
243 :
244 : typedef typename TGridFunction::domain_type dom_type;
245 : typedef typename TGridFunction::algebra_type algebra_type;
246 : typedef typename dom_type::position_attachment_type position_attachment_type;
247 : typedef typename TGridFunction::approximation_space_type approx_space_type;
248 : typedef typename TGridFunction::template dim_traits<trueDim>::grid_base_object elem_type;
249 :
250 : SmartPtr<approx_space_type> srcApproxSpace = u->approx_space();
251 : SmartPtr<dom_type> srcDom = srcApproxSpace->domain();
252 :
253 : // select surface elements in old grid
254 0 : MultiGrid& srcGrid = *srcDom->grid();
255 0 : Selector srcSel(srcGrid);
256 0 : srcSel.select(u->template begin<elem_type>(), u->template end<elem_type>());
257 :
258 : // create new domain
259 : dom_type* dom_ptr;
260 0 : try {dom_ptr = new dom_type();}
261 0 : UG_CATCH_THROW("Temporary domain could not be created.");
262 : SmartPtr<dom_type> destDom = make_sp(dom_ptr);
263 :
264 : // copy grid from old domain to new domain (into a flat grid!)
265 0 : try {CopySelected(destDom, srcDom, srcSel);}
266 0 : UG_CATCH_THROW("Temporary grid could not be created.");
267 :
268 : // refine
269 : #ifdef UG_PARALLEL
270 : ParallelGlobalRefiner_MultiGrid refiner(*destDom->distributed_grid_manager());
271 : #else
272 0 : MultiGrid& destGrid = *destDom->grid();
273 0 : GlobalMultiGridRefiner refiner(destGrid);
274 : #endif
275 0 : size_t numRefs = (size_t) ceil(log2(order));
276 0 : for (size_t iref = 0; iref < numRefs; ++iref)
277 : {
278 0 : try {refiner.refine();}
279 0 : UG_CATCH_THROW("Refinement step " << iref << " could not be carried out.");
280 : }
281 :
282 : // retain function group for functions being exported
283 0 : FunctionGroup fg(srcApproxSpace->dof_distribution_info(), vFct);
284 :
285 : // create approx space and add functions
286 : approx_space_type* approx_ptr;
287 0 : try {approx_ptr = new approx_space_type(destDom, algebra_type::get_type());}
288 0 : UG_CATCH_THROW("Temporary approximation space could not be created.");
289 : SmartPtr<approx_space_type> destApproxSpace = make_sp(approx_ptr);
290 :
291 0 : for (size_t fct = 0; fct < fg.size(); ++fct)
292 : {
293 0 : if (fg.function_pattern()->is_def_everywhere(fg.unique_id(fct)))
294 0 : destApproxSpace->add(fg.name(fct), "Lagrange", 1);
295 : else
296 : {
297 0 : int num_subsets = fg.function_pattern()->num_subsets();
298 : std::string subsets;
299 0 : for (int si = 0; si < num_subsets; ++si)
300 : {
301 0 : if (fg.function_pattern()->is_def_in_subset(fct, si))
302 0 : subsets += std::string(",") + fg.function_pattern()->subset_name(si);
303 : }
304 0 : if (!subsets.empty())
305 0 : subsets.erase(0,1); // delete leading ","
306 :
307 0 : destApproxSpace->add(fg.name(fct), "Lagrange", 1, subsets.c_str());
308 : }
309 : }
310 0 : destApproxSpace->init_top_surface();
311 :
312 : TGridFunction* gridFct_ptr;
313 0 : try {gridFct_ptr = new TGridFunction(destApproxSpace);}
314 0 : UG_CATCH_THROW("Temporary grid function could not be created.");
315 : SmartPtr<TGridFunction> u_new = make_sp(gridFct_ptr);
316 :
317 : // interpolate onto new grid
318 0 : for (size_t fct = 0; fct < fg.size(); ++fct)
319 : {
320 0 : const LFEID lfeid = u_new->dof_distribution()->lfeid(fct);
321 :
322 0 : GlobalGridFunctionNumberData<TGridFunction, trueDim> ggfnd =
323 : GlobalGridFunctionNumberData<TGridFunction, trueDim>(u, fg.name(fct));
324 :
325 : // iterate over DoFs in new function and evaluate
326 : // should be vertices only for Lagrange-1
327 0 : if (u_new->max_dofs(VERTEX))
328 : interpolate_from_original_fct<Vertex, TGridFunction, GlobalGridFunctionNumberData<TGridFunction, trueDim> >
329 0 : (u_new, ggfnd, fct, lfeid);
330 : /*
331 : if (u_new->max_dofs(EDGE))
332 : interpolate_from_original_fct<Edge, TGridFunction, GlobalGridFunctionNumberData<TGridFunction, trueDim> >
333 : (u_new, ggfnd, fct, lfeid);
334 : if (u_new->max_dofs(FACE))
335 : interpolate_from_original_fct<Face, TGridFunction, GlobalGridFunctionNumberData<TGridFunction, trueDim> >
336 : (u_new, ggfnd, fct, lfeid);
337 : if (u_new->max_dofs(VOLUME))
338 : interpolate_from_original_fct<Volume, TGridFunction, GlobalGridFunctionNumberData<TGridFunction, trueDim> >
339 : (u_new, ggfnd, fct, lfeid);
340 : */
341 : }
342 :
343 : #ifdef UG_PARALLEL
344 : // copy storage type and layouts
345 : u_new->set_storage_type(u->get_storage_mask());
346 : u_new->set_layouts(u->layouts());
347 : #endif
348 :
349 : // print out new grid function on desired subsets
350 0 : bool printAll = ssg.size() == (size_t) u->subset_handler()->num_subsets();
351 0 : if (printAll)
352 0 : vtkOutput->print(filename, *u_new, step, time);
353 : else
354 : {
355 0 : for (size_t s = 0; s < ssg.size(); ++s)
356 0 : vtkOutput->print_subset(filename, *u_new, ssg[s], step, time);
357 : }
358 0 : }
359 :
360 :
361 : /// export a solution from a high-order ansatz space to vtk file(s)
362 : /**
363 : * This function will create a temporary domain, copy all elements from the domain
364 : * which the grid function u is defined on to the temporary domain and then refine
365 : * the resulting grid until it has at least as many vertices as the original grid
366 : * functions has unknowns (e.g. a grid for a function of order 2 would be refined
367 : * once, a grid for a function of order 4 would be refined twice, and so on).
368 : * After refinement, a temporary grid function of order 1 (Lagrange) is defined on
369 : * the refined grid and its values interpolated from the original function u.
370 : * The temporary first-order grid function is then exported to vtk using the usual
371 : * mechanisms.
372 : *
373 : * @param u original high-order grid function to be exported
374 : * @param vFct vector of function names (contained in grid function) to be exported
375 : * @param order order to be used
376 : * @param vtkOutput VTKOutput object to use for export of linearized function
377 : * @param filename file name to be used in export
378 : *
379 : * @todo The order parameter might be left out and determined automatically from the
380 : * grid function.
381 : */
382 : template <typename TGridFunction>
383 0 : void vtk_export_ho
384 : (
385 : SmartPtr<TGridFunction> u,
386 : const std::vector<std::string>& vFct,
387 : size_t order,
388 : SmartPtr<VTKOutput<TGridFunction::domain_type::dim> > vtkOutput,
389 : const char* filename,
390 : size_t step,
391 : number time,
392 : const std::vector<std::string>& vSubset
393 : )
394 : {
395 : // construct subset group from given subsets
396 0 : SubsetGroup ssg(u->subset_handler());
397 : try
398 : {
399 0 : if (vSubset.empty())
400 0 : ssg.add_all();
401 : else
402 0 : ssg.add(vSubset);
403 : }
404 0 : UG_CATCH_THROW("Subsets are faulty.");
405 :
406 : // find highest dim that contains any elements
407 0 : MultiGrid& srcGrid = *u->approx_space()->domain()->grid();
408 0 : if (srcGrid.num_volumes())
409 : vtk_export_ho<TGridFunction, TGridFunction::dim >= 3 ? 3 : TGridFunction::dim>
410 0 : (u, vFct, order, vtkOutput, filename, step, time, ssg);
411 0 : else if (srcGrid.num_faces())
412 : vtk_export_ho<TGridFunction, TGridFunction::dim >= 2 ? 2 : TGridFunction::dim>
413 0 : (u, vFct, order, vtkOutput, filename, step, time, ssg);
414 0 : else if (srcGrid.num_edges())
415 : vtk_export_ho<TGridFunction, TGridFunction::dim >= 1 ? 1 : TGridFunction::dim>
416 0 : (u, vFct, order, vtkOutput, filename, step, time, ssg);
417 0 : else if (srcGrid.num_vertices())
418 : vtk_export_ho<TGridFunction, TGridFunction::dim >= 0 ? 0 : TGridFunction::dim>
419 0 : (u, vFct, order, vtkOutput, filename, step, time, ssg);
420 : else
421 0 : vtk_export_ho<TGridFunction, TGridFunction::dim>(u, vFct, order, vtkOutput, filename, step, time, ssg);
422 0 : }
423 :
424 :
425 : template <typename TGridFunction>
426 : void vtk_export_ho
427 : (
428 : SmartPtr<TGridFunction> u,
429 : const std::vector<std::string>& vFct,
430 : size_t order,
431 : SmartPtr<VTKOutput<TGridFunction::domain_type::dim> > vtkOutput,
432 : const char* filename,
433 : const std::vector<std::string>& vSubset
434 : )
435 : {
436 : vtk_export_ho(u, vFct, order, vtkOutput, filename, 0, 0.0, vSubset);
437 : }
438 :
439 :
440 : template <typename TGridFunction>
441 0 : void vtk_export_ho
442 : (
443 : SmartPtr<TGridFunction> u,
444 : const std::vector<std::string>& vFct,
445 : size_t order,
446 : SmartPtr<VTKOutput<TGridFunction::domain_type::dim> > vtkOutput,
447 : const char* filename,
448 : size_t step,
449 : number time
450 : )
451 : {
452 0 : vtk_export_ho(u, vFct, order, vtkOutput, filename, step, time, std::vector<std::string>());
453 0 : }
454 :
455 :
456 : template <typename TGridFunction>
457 0 : void vtk_export_ho
458 : (
459 : SmartPtr<TGridFunction> u,
460 : const std::vector<std::string>& vFct,
461 : size_t order,
462 : SmartPtr<VTKOutput<TGridFunction::domain_type::dim> > vtkOutput,
463 : const char* filename
464 : )
465 : {
466 0 : vtk_export_ho(u, vFct, order, vtkOutput, filename, 0, 0.0, std::vector<std::string>());
467 0 : }
468 :
469 :
470 : } // end namespace ug
471 :
472 :
473 : #endif // __H__UG__LIB_DISC__IO__VTK_EXPORT_HO__
|
