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 : #ifndef __H__UG__LIB_DISC__FUNCTION_SPACE__GRID_FUNCTION_UTIL__
34 : #define __H__UG__LIB_DISC__FUNCTION_SPACE__GRID_FUNCTION_UTIL__
35 :
36 : #include <vector>
37 : #include <string>
38 : #include <cmath> // for isinf, isnan
39 :
40 :
41 : #include "common/util/file_util.h"
42 :
43 : #include "lib_algebra/cpu_algebra/sparsematrix_print.h"
44 : #include "lib_algebra/operator/interface/matrix_operator.h"
45 : #include "lib_algebra/operator/debug_writer.h"
46 : #include "lib_algebra/operator/vector_writer.h"
47 : #include "lib_algebra/common/matrixio/matrix_io_mtx.h"
48 : #include "lib_algebra/common/connection_viewer_output.h"
49 : #include "lib_algebra/common/csv_gnuplot_output.h"
50 : #include "lib_disc/io/vtkoutput.h"
51 : #include "lib_disc/spatial_disc/constraints/constraint_interface.h"
52 : #include "lib_disc/dof_manager/dof_distribution.h"
53 : #include "lib_disc/spatial_disc/user_data/user_data.h"
54 : #include "lib_disc/common/groups_util.h"
55 : #include "lib_disc/common/geometry_util.h"
56 : #include "lib_disc/function_spaces/integrate.h"
57 : #include "lib_grid/algorithms/debug_util.h" // for ElementDebugInfo
58 : #include "lib_grid/tools/periodic_boundary_manager.h"
59 :
60 : #include "grid_function.h"
61 : #include "dof_position_util.h"
62 :
63 : #ifdef UG_PARALLEL
64 : #include "pcl/pcl.h"
65 : #endif
66 :
67 : namespace ug {
68 :
69 : #ifndef isnan
70 : using boost::math::isnan;
71 : #endif
72 :
73 : #ifndef isinf
74 : using boost::math::isinf;
75 : #endif
76 :
77 :
78 : template <typename TBaseElem, typename TGridFunction>
79 0 : static void ScaleGFOnElems
80 : (
81 : ConstSmartPtr<DoFDistribution> dd,
82 : SmartPtr<TGridFunction> vecOut,
83 : ConstSmartPtr<TGridFunction> vecIn,
84 : const std::vector<number>& vScale
85 : )
86 : {
87 : typename DoFDistribution::traits<TBaseElem>::const_iterator iter, iterEnd;
88 : std::vector<DoFIndex> vInd;
89 :
90 : try
91 : {
92 : // iterate all elements (including SHADOW_RIM_COPY!)
93 0 : iter = dd->template begin<TBaseElem>(SurfaceView::ALL);
94 0 : iterEnd = dd->template end<TBaseElem>(SurfaceView::ALL);
95 0 : for (; iter != iterEnd; ++iter)
96 : {
97 0 : for (size_t fi = 0; fi < dd->num_fct(); ++fi)
98 : {
99 0 : size_t nInd = dd->inner_dof_indices(*iter, fi, vInd);
100 :
101 : // remember multi indices
102 0 : for (size_t dof = 0; dof < nInd; ++dof)
103 0 : DoFRef(*vecOut, vInd[dof]) = vScale[fi] * DoFRef(*vecIn, vInd[dof]);
104 : }
105 : }
106 : }
107 0 : UG_CATCH_THROW("Error while scaling vector.")
108 0 : }
109 :
110 :
111 : /**
112 : * \brief Scales all functions contained in a grid function.
113 : *
114 : * Each function has a separate scaling factor.
115 : *
116 : * \param scaledVecOut the scaled grid function (output)
117 : * \param vecIn the original grid function (input)
118 : * \param scalingFactors vector of scales for each of the composite functions
119 : */
120 : template <typename TGridFunction>
121 0 : void ScaleGF
122 : (
123 : SmartPtr<TGridFunction> scaledVecOut,
124 : ConstSmartPtr<TGridFunction> vecIn,
125 : const std::vector<number>& scalingFactors
126 : )
127 : {
128 : // check that the correct numbers of scaling factors are given
129 : size_t n = scalingFactors.size();
130 0 : UG_COND_THROW(n != vecIn->num_fct(), "Number of scaling factors (" << n << ") "
131 : "does not match number of functions given in dimless vector (" << vecIn->num_fct() << ").");
132 :
133 : // check that input and output vectors have the same number of components and dofs
134 0 : UG_COND_THROW(n != scaledVecOut->num_fct(), "Input and output vectors do not have "
135 : "the same number of functions (" << n << " vs. " << scaledVecOut->num_fct() << ").");
136 0 : for (size_t fct = 0; fct < n; ++fct)
137 : {
138 0 : UG_COND_THROW(vecIn->num_dofs(fct) != scaledVecOut->num_dofs(fct),
139 : "Input and output vectors do not have the same number of DoFs for function " << fct
140 : << " (" << vecIn->num_dofs(fct) << " vs. " << scaledVecOut->num_dofs(fct) << ").");
141 : }
142 :
143 : ConstSmartPtr<DoFDistribution> dd = vecIn->dof_distribution();
144 :
145 0 : if (dd->max_dofs(VERTEX))
146 0 : ScaleGFOnElems<Vertex, TGridFunction>(dd, scaledVecOut, vecIn, scalingFactors);
147 0 : if (dd->max_dofs(EDGE))
148 0 : ScaleGFOnElems<Edge, TGridFunction>(dd, scaledVecOut, vecIn, scalingFactors);
149 0 : if (dd->max_dofs(FACE))
150 0 : ScaleGFOnElems<Face, TGridFunction>(dd, scaledVecOut, vecIn, scalingFactors);
151 0 : if (dd->max_dofs(VOLUME))
152 0 : ScaleGFOnElems<Volume, TGridFunction>(dd, scaledVecOut, vecIn, scalingFactors);
153 0 : }
154 :
155 :
156 :
157 :
158 : ////////////////////////////////////////////////////////////////////////////////
159 : // AverageComponent
160 : ////////////////////////////////////////////////////////////////////////////////
161 :
162 : /**
163 : * Subtracts a given value from all DoFs of a given grid function associated with
164 : * a given element type.
165 : *
166 : * \tparam TGridFunction grid function type
167 : * \tparam TBaseElem type of the elements the DoFs are associated with
168 : */
169 : template<typename TGridFunction, typename TBaseElem>
170 0 : void SubtractValueFromComponent(SmartPtr<TGridFunction> spGF, size_t fct, number sub)
171 : {
172 : typedef TGridFunction GF;
173 : typedef typename GF::template traits<TBaseElem>::const_iterator iter_type;
174 :
175 : iter_type iter = spGF->template begin<TBaseElem>();
176 : iter_type iterEnd = spGF->template end<TBaseElem>();
177 :
178 0 : PeriodicBoundaryManager* pbm = spGF->domain()->grid()->periodic_boundary_manager();
179 :
180 : // loop elems
181 : std::vector<DoFIndex> vMultInd;
182 0 : for(; iter != iterEnd; ++iter)
183 : {
184 : // get element
185 : TBaseElem* elem = *iter;
186 :
187 : // skip periodic ghosts
188 0 : if (pbm && pbm->is_slave(elem)) continue;
189 :
190 : // get global indices
191 : spGF->inner_dof_indices(elem, fct, vMultInd);
192 :
193 : // sum up value
194 0 : for(size_t i = 0; i < vMultInd.size(); ++i)
195 : {
196 0 : DoFRef(*spGF, vMultInd[i]) -= sub;
197 : }
198 : }
199 0 : }
200 :
201 : /**
202 : * Subtracts the average of a given grid function (computed by taking into account
203 : * the shape functions), shifted by the value of 'mean', from all the DoFs of this
204 : * grid function. Typically, it should result in a grid function with the average
205 : * equal to 'mean'.
206 : *
207 : * \remark Note that mean is not the desired average but the one multiplied by the area of the domain.
208 : *
209 : * \tparam TGridFunction type of the grid function
210 : */
211 : template<typename TGridFunction>
212 0 : void AdjustMeanValue
213 : (
214 : SmartPtr<TGridFunction> spGF, ///< the grid function
215 : const std::vector<std::string>& vCmp, ///< the component to adjust the average for
216 : number mean ///< the desired mean value
217 : )
218 : {
219 : typedef TGridFunction GF;
220 : PROFILE_FUNC_GROUP("gmg");
221 :
222 0 : if(vCmp.empty())
223 0 : return;
224 :
225 0 : if(spGF.invalid())
226 0 : UG_THROW("AverageComponent: expects a valid GridFunction.");
227 :
228 : ConstSmartPtr<DoFDistributionInfo> ddinfo =
229 0 : spGF->approx_space()->dof_distribution_info();
230 :
231 : // compute integral of components
232 0 : const number area = Integral(1.0, spGF);
233 0 : std::vector<number> vIntegral(vCmp.size(), 0.0);
234 0 : for(size_t f = 0; f < vCmp.size(); f++){
235 : const size_t fct = spGF->fct_id_by_name(vCmp[f].c_str());
236 0 : vIntegral[f] = Integral(spGF, vCmp[f].c_str(), NULL, ddinfo->lfeid(fct).order());
237 : }
238 :
239 : // subtract value
240 0 : for(size_t f = 0; f < vCmp.size(); f++)
241 : {
242 0 : const number sub = (vIntegral[f] - mean) / area;
243 : const size_t fct = spGF->fct_id_by_name(vCmp[f].c_str());
244 :
245 0 : if(ddinfo->max_fct_dofs(fct, VERTEX)) SubtractValueFromComponent<GF, Vertex>(spGF, fct, sub);
246 0 : if(ddinfo->max_fct_dofs(fct, EDGE)) SubtractValueFromComponent<GF, Edge>(spGF, fct, sub);
247 0 : if(ddinfo->max_fct_dofs(fct, FACE)) SubtractValueFromComponent<GF, Face>(spGF, fct, sub);
248 0 : if(ddinfo->max_fct_dofs(fct, VOLUME)) SubtractValueFromComponent<GF, Volume>(spGF, fct, sub);
249 : }
250 0 : }
251 :
252 : /**
253 : * Subtracts the average of a given grid function (computed by taking into account
254 : * the shape functions) from all the DoFs of this grid function. Typically, it
255 : * should result in a grid function with the zero average.
256 : *
257 : * \tparam TGridFunction type of the grid function
258 : */
259 : template<typename TGridFunction>
260 0 : void AdjustMeanValue(SmartPtr<TGridFunction> spGF, const std::vector<std::string>& vCmp)
261 : {
262 0 : AdjustMeanValue(spGF, vCmp, 0.0);
263 0 : }
264 :
265 : /**
266 : * Subtracts the average of a given grid function (computed by taking into account
267 : * the shape functions) from all the DoFs of this grid function. Typically, it
268 : * should result in a grid function with the zero average.
269 : *
270 : * \tparam TGridFunction type of the grid function
271 : */
272 : template<typename TGridFunction>
273 0 : void AdjustMeanValue(SmartPtr<TGridFunction> spGF, const std::string& fcts)
274 : {
275 0 : AdjustMeanValue(spGF, TokenizeTrimString(fcts), 0.0);
276 0 : }
277 :
278 : /**
279 : * Subtracts the average of a given grid function (computed by taking into account
280 : * the shape functions), shifted by the value of 'mean', from all the DoFs of this
281 : * grid function. Typically, it should result in a grid function with the average
282 : * equal to 'mean'.
283 : *
284 : * \remark Note that mean is not the desired average but the one multiplied by the area of the domain.
285 : *
286 : * \tparam TGridFunction type of the grid function
287 : */
288 : template<typename TGridFunction>
289 0 : void AdjustMeanValue(SmartPtr<TGridFunction> spGF, const std::string& fcts, number mean)
290 : {
291 0 : AdjustMeanValue(spGF, TokenizeTrimString(fcts), mean);
292 0 : }
293 :
294 : ////////////////////////////////////////////////////////////////////////////////
295 : // Summing up components
296 : ////////////////////////////////////////////////////////////////////////////////
297 :
298 : /**
299 : * Sums values of a grid function at all given elements (like vertices)
300 : *
301 : * \tparam TGridFunction the grid function type
302 : * \tparam TBaseElem the given element type (like 'Vertex')
303 : */
304 : template <typename TGridFunction, typename TBaseElem>
305 0 : number SumGFValuesAt
306 : (
307 : TGridFunction * u, ///< the grid function
308 : size_t fct ///< index of the function
309 : )
310 : {
311 : typedef typename TGridFunction::template traits<TBaseElem>::const_iterator t_elem_iter;
312 :
313 : std::vector<DoFIndex> ind;
314 : number sum (1);
315 :
316 : sum = 0;
317 : // Loop the elements in the subset
318 0 : for (t_elem_iter vi = u->template begin<TBaseElem> ();
319 0 : vi != u->template end<TBaseElem> (); ++vi)
320 : {
321 : TBaseElem * vert = *vi;
322 :
323 : // indices at this element
324 : u->inner_dof_indices (vert, fct, ind);
325 :
326 0 : if (ind.size () != 1)
327 0 : UG_THROW ("SumGFValuesAt: The function must be scalar!");
328 :
329 : // add the contribution
330 0 : sum += DoFRef (*u, ind [0]);
331 : }
332 :
333 : #ifdef UG_PARALLEL
334 : pcl::ProcessCommunicator procComm;
335 : sum = procComm.allreduce (sum, PCL_RO_SUM);
336 : #endif
337 :
338 0 : return sum;
339 : }
340 :
341 : /**
342 : * Sums values of a grid function at all given elements (like vertices).
343 : *
344 : * This version takes symbolic data.
345 : *
346 : * \tparam TGridFunction the grid function type
347 : * \tparam TBaseElem the given element type (like 'Vertex')
348 : */
349 : template <typename TGridFunction, typename TBaseElem>
350 0 : number SumGFValuesAt
351 : (
352 : TGridFunction * u, ///< the grid function
353 : const char * fct_names ///< index of the function
354 : )
355 : {
356 : // Get the function index
357 : std::vector<std::string> vfctNames;
358 0 : TokenizeString (fct_names, vfctNames);
359 0 : if (vfctNames.size () != 1)
360 0 : UG_THROW ("SumGFValuesAt: Exactly one function name must be specified.");
361 0 : FunctionGroup fctGroup (u->function_pattern ());
362 0 : fctGroup.add (vfctNames [0]);
363 :
364 : // Compute the sum
365 0 : return SumGFValuesAt<TGridFunction, TBaseElem> (u, fctGroup[0]);
366 0 : }
367 :
368 : /**
369 : * Sums values of a grid function at given elements (like vertices) in given subsets
370 : *
371 : * \tparam TGridFunction the grid function type
372 : * \tparam TBaseElem the given element type (like 'Vertex')
373 : */
374 : template <typename TGridFunction, typename TBaseElem>
375 0 : number SumGFValuesAt
376 : (
377 : TGridFunction * u, ///< the grid function
378 : size_t fct, ///< index of the function
379 : SubsetGroup & ssGroup ///< the subsets
380 : )
381 : {
382 : typedef typename TGridFunction::template traits<TBaseElem>::const_iterator t_elem_iter;
383 :
384 : std::vector<DoFIndex> ind;
385 : number sum (1);
386 :
387 : // Loop the selected subsets
388 : sum = 0;
389 0 : for (size_t i = 0; i < ssGroup.size (); i++)
390 : {
391 : int ssi = ssGroup [i];
392 :
393 : // Loop the elements in the subset
394 0 : for (t_elem_iter vi = u->template begin<TBaseElem> (ssi);
395 0 : vi != u->template end<TBaseElem> (ssi); ++vi)
396 : {
397 : TBaseElem * vert = *vi;
398 :
399 : // indices at this element
400 : u->inner_dof_indices (vert, fct, ind);
401 :
402 0 : if (ind.size () != 1)
403 0 : UG_THROW ("SumGFValuesAt: The function must be scalar!");
404 :
405 : // add the contribution
406 0 : sum += DoFRef (*u, ind [0]);
407 : }
408 : }
409 :
410 : #ifdef UG_PARALLEL
411 : pcl::ProcessCommunicator procComm;
412 : sum = procComm.allreduce (sum, PCL_RO_SUM);
413 : #endif
414 :
415 0 : return sum;
416 : }
417 :
418 : /**
419 : * Sums values of a grid function at given elements (like vertices) in given subsets.
420 : *
421 : * This version takes symbolic data.
422 : *
423 : * \tparam TGridFunction the grid function type
424 : * \tparam TBaseElem the given element type (like 'Vertex')
425 : */
426 : template <typename TGridFunction, typename TBaseElem>
427 0 : number SumGFValuesAt
428 : (
429 : TGridFunction * u, ///< the grid function
430 : const char * fct_names, ///< index of the function
431 : const char * subset_names ///< the subsets
432 : )
433 : {
434 : // Get the function index
435 : std::vector<std::string> vfctNames;
436 0 : TokenizeString (fct_names, vfctNames);
437 0 : if (vfctNames.size () != 1)
438 0 : UG_THROW ("SumGFValuesAt: Exactly one function name must be specified.");
439 0 : FunctionGroup fctGroup (u->function_pattern ());
440 0 : fctGroup.add (vfctNames [0]);
441 :
442 : // Get the subset indices
443 : std::vector<std::string> vssNames;
444 0 : TokenizeString (subset_names, vssNames);
445 0 : SubsetGroup ssGroup (u->domain()->subset_handler ());
446 0 : ssGroup.add (vssNames);
447 :
448 : // Compute the sum
449 0 : return SumGFValuesAt<TGridFunction, TBaseElem> (u, fctGroup[0], ssGroup);
450 0 : }
451 :
452 : ////////////////////////////////////////////////////////////////////////////////
453 : // Checking components for nan and inf
454 : ////////////////////////////////////////////////////////////////////////////////
455 :
456 : /**
457 : * Checks values of a grid function at given elements (like vertices) for nan and inf
458 : *
459 : * \tparam TGridFunction the grid function type
460 : * \tparam TBaseElem the given element type (like 'Vertex')
461 : */
462 : template <typename TGridFunction, typename TBaseElem>
463 0 : bool CheckGFforNaN
464 : (
465 : const TGridFunction * u, ///< the grid function
466 : size_t fct ///< index of the function
467 : )
468 : {
469 : typedef typename TGridFunction::template traits<TBaseElem>::const_iterator t_elem_iter;
470 :
471 : std::vector<DoFIndex> ind;
472 :
473 : // Loop the elements in the subset
474 0 : for (t_elem_iter vi = u->template begin<TBaseElem> ();
475 0 : vi != u->template end<TBaseElem> (); ++vi)
476 : {
477 : TBaseElem * vert = *vi;
478 :
479 : // indices at this element
480 : u->inner_dof_indices (vert, fct, ind);
481 :
482 0 : if (ind.size () != 1)
483 0 : UG_THROW ("CheckGFforNaN: The function must be scalar!");
484 :
485 : // check the value
486 0 : number value = DoFRef (*u, ind [0]);
487 0 : if (isnan (value))
488 : {
489 0 : int si = u->domain()->subset_handler()->get_subset_index (vert);
490 0 : UG_LOG ("nan at index " << ind [0] << ", grid data idx " << vert->grid_data_index ());
491 0 : if (si >= 0)
492 0 : UG_LOG (", subset " << u->domain()->subset_handler()->get_subset_name (si))
493 0 : UG_LOG ('\n');
494 0 : return true;
495 : }
496 0 : if (isinf (value))
497 : {
498 0 : int si = u->domain()->subset_handler()->get_subset_index (vert);
499 0 : UG_LOG ("inf at index " << ind [0] << ", grid data idx " << vert->grid_data_index ());
500 0 : if (si >= 0)
501 0 : UG_LOG (", subset " << u->domain()->subset_handler()->get_subset_name (si))
502 0 : UG_LOG ('\n');
503 : return true;
504 : }
505 : }
506 :
507 0 : return false;
508 : }
509 :
510 : /**
511 : * Checks values of a grid function at given elements (like vertices) for nan and inf
512 : *
513 : * This version takes symbolic data.
514 : *
515 : * \tparam TGridFunction the grid function type
516 : * \tparam TBaseElem the given element type (like 'Vertex')
517 : */
518 : template <typename TGridFunction, typename TBaseElem>
519 0 : bool CheckGFforNaN
520 : (
521 : const TGridFunction * u, ///< the grid function
522 : const char * fct_names ///< names of the functions
523 : )
524 : {
525 : // Get the function index
526 : std::vector<std::string> vfctNames;
527 0 : TokenizeString (fct_names, vfctNames);
528 0 : FunctionGroup fctGroup (u->function_pattern ());
529 0 : for (size_t i = 0; i < vfctNames.size (); i++)
530 0 : fctGroup.add (vfctNames [i]);
531 :
532 : // Check the functions
533 : bool result = false;
534 0 : for (size_t i = 0; i < vfctNames.size (); i++)
535 : {
536 : UG_LOG ("Checking " << vfctNames[i] << " ... ");
537 0 : if (CheckGFforNaN<TGridFunction, TBaseElem> (u, fctGroup[i]))
538 : result = true;
539 : else
540 : UG_LOG ("OK\n");
541 : }
542 0 : return result;
543 0 : }
544 :
545 : /**
546 : * Checks values of a grid function at given elements (like vertices) for nan and inf
547 : *
548 : * This version takes a function group.
549 : *
550 : * \tparam TGridFunction the grid function type
551 : * \tparam TBaseElem the given element type (like 'Vertex')
552 : */
553 : template <typename TGridFunction, typename TBaseElem>
554 : bool CheckGFforNaN
555 : (
556 : const TGridFunction * u, ///< the grid function
557 : const FunctionGroup & fctGroup ///< names of the functions
558 : )
559 : {
560 : // Check the functions
561 : bool result = false;
562 : for (size_t i = 0; i < fctGroup.size (); i++)
563 : {
564 : UG_LOG ("Checking fct #" << i << " ... ");
565 : if (CheckGFforNaN<TGridFunction, TBaseElem> (u, fctGroup[i]))
566 : result = true;
567 : else
568 : UG_LOG ("OK\n");
569 : }
570 : return result;
571 : }
572 :
573 : // //////////////////////////////////////////////////////////////////////////////
574 : // Check values are within bounds
575 : // //////////////////////////////////////////////////////////////////////////////
576 : template <typename TGridFunction, typename TBaseElem>
577 0 : bool CheckGFValuesWithinBounds
578 : (
579 : ConstSmartPtr<TGridFunction> u, ///< the grid function
580 : size_t cmp, ///< the function component to be checked
581 : number lowerBnd,
582 : number upperBnd
583 : )
584 : {
585 : typedef typename TGridFunction::template traits<TBaseElem>::const_iterator elem_it;
586 :
587 : bool ret = true;
588 :
589 : // loop the elements in the subset
590 : std::vector<DoFIndex> vDI;
591 : elem_it it = u->template begin<TBaseElem>();
592 : elem_it itEnd = u->template end<TBaseElem>();
593 0 : for (; it != itEnd; ++it)
594 : {
595 : TBaseElem* elem = *it;
596 :
597 : // loop indices at this element
598 : const size_t nInd = u->inner_dof_indices(elem, cmp, vDI, true);
599 0 : for (size_t i = 0; i < nInd; ++i)
600 : {
601 : const number& val = DoFRef(*u, vDI[i]);
602 0 : if (val < lowerBnd || val > upperBnd)
603 : {
604 0 : UG_LOG_ALL_PROCS("Function value for component " << cmp << " (" << val << ") "
605 : "is outside the specified range [" << lowerBnd << ", " << upperBnd << "] "
606 : "at " << ElementDebugInfo(*u->domain()->grid(), elem) << std::endl);
607 : ret = false;
608 : }
609 : }
610 : }
611 :
612 : #ifdef UG_PARALLEL
613 : if (pcl::NumProcs() > 1)
614 : {
615 : pcl::ProcessCommunicator pc;
616 : int retInt = ret;
617 : ret = pc.allreduce(retInt, PCL_RO_BAND);
618 : }
619 : #endif
620 :
621 0 : return ret;
622 : }
623 :
624 : template <typename TGridFunction>
625 0 : bool CheckGFValuesWithinBounds
626 : (
627 : ConstSmartPtr<TGridFunction> u, ///< the grid function
628 : size_t cmp, ///< the function component to be checked
629 : number lowerBnd,
630 : number upperBnd
631 : )
632 : {
633 : bool ret = true;
634 0 : if (u->max_fct_dofs(cmp, 0))
635 0 : ret = ret && CheckGFValuesWithinBounds<TGridFunction, Vertex>(u, cmp, lowerBnd, upperBnd);
636 0 : if (u->max_fct_dofs(cmp, 1))
637 0 : ret = ret && CheckGFValuesWithinBounds<TGridFunction, Edge>(u, cmp, lowerBnd, upperBnd);
638 0 : if (u->max_fct_dofs(cmp, 2))
639 0 : ret = ret && CheckGFValuesWithinBounds<TGridFunction, Face>(u, cmp, lowerBnd, upperBnd);
640 0 : if (u->max_fct_dofs(cmp, 3))
641 0 : ret = ret && CheckGFValuesWithinBounds<TGridFunction, Volume>(u, cmp, lowerBnd, upperBnd);
642 :
643 0 : return ret;
644 : }
645 :
646 : template <typename TGridFunction>
647 0 : bool CheckGFValuesWithinBounds
648 : (
649 : ConstSmartPtr<TGridFunction> u, ///< the grid function
650 : const char* fctNames, ///< the functions to be checked
651 : number lowerBnd,
652 : number upperBnd
653 : )
654 : {
655 : std::vector<std::string> vFctNames;
656 0 : TokenizeTrimString (fctNames, vFctNames);
657 0 : FunctionGroup fctGroup (u->function_pattern());
658 : const size_t nFct = vFctNames.size();
659 0 : for (size_t f = 0; f < nFct; ++f)
660 : {
661 0 : try {fctGroup.add(vFctNames[f]);}
662 0 : UG_CATCH_THROW("Could not add function " << vFctNames[f] << " to function group.");
663 : }
664 :
665 : bool ret = true;
666 : const size_t fctGrpSz = fctGroup.size();
667 0 : for (size_t f = 0; f < fctGrpSz; ++f)
668 0 : ret = ret && CheckGFValuesWithinBounds(u, fctGroup[f], lowerBnd, upperBnd);
669 :
670 0 : return ret;
671 0 : }
672 :
673 :
674 :
675 :
676 : ////////////////////////////////////////////////////////////////////////////////
677 : // Writing algebra
678 : ////////////////////////////////////////////////////////////////////////////////
679 :
680 : template<typename TDomain>
681 0 : void WriteAlgebraIndices(std::string name, ConstSmartPtr<TDomain> domain, ConstSmartPtr<DoFDistribution> dd)
682 : {
683 : PROFILE_FUNC_GROUP("debug");
684 :
685 : #ifdef UG_PARALLEL
686 : name = ConnectionViewer::GetParallelName(name, dd->layouts()->proc_comm());
687 : #endif
688 :
689 : std::vector<size_t> fctIndex;
690 : std::vector<std::string> fctNames;
691 :
692 0 : ExtractAlgebraIndices<TDomain>(domain, dd, fctIndex);
693 :
694 : size_t numFct = dd->num_fct();
695 0 : if(numFct <= 1) return;
696 :
697 0 : fctNames.resize(numFct);
698 0 : for(size_t i=0; i<numFct; i++)
699 0 : fctNames[i] = dd->name(i);
700 :
701 0 : name.append(".indices");
702 0 : std::fstream file(name.c_str(), std::ios::out);
703 :
704 : //std::cout << "name is " << name << "\n";
705 0 : file << "NUMDOF " << fctNames.size() << "\n";
706 0 : for(size_t i=0; i<numFct; i++)
707 0 : file << fctNames[i] << "\n";
708 :
709 0 : for(size_t i=0; i<fctIndex.size(); i++)
710 0 : file << fctIndex[i] << "\n";
711 0 : }
712 :
713 : template<class TFunction>
714 0 : void WriteMatrixToConnectionViewer(const char *filename,
715 : const typename TFunction::algebra_type::matrix_type &A,
716 : const TFunction &u) {
717 :
718 : PROFILE_FUNC_GROUP("debug");
719 : // check name
720 0 : if ( !FileTypeIs( filename, ".mat") ) {
721 0 : UG_THROW( "Only '.mat' format supported for matrices, but"
722 : " filename is '" << filename << "'." );
723 : }
724 :
725 : // position array
726 : const static int dim = TFunction::domain_type::dim;
727 : std::vector<MathVector<dim> > vPos;
728 0 : ExtractPositions(u, vPos);
729 :
730 : // write matrix
731 0 : if(vPos.empty())
732 0 : ConnectionViewer::WriteMatrixPar( filename, A, (MathVector<dim>*)NULL, dim );
733 : else
734 0 : ConnectionViewer::WriteMatrixPar( filename, A, &vPos[0], dim );
735 :
736 0 : WriteAlgebraIndices(filename, u.domain(),u.dof_distribution());
737 :
738 0 : }
739 :
740 : template<typename TGridFunction>
741 0 : void SaveMatrixForConnectionViewer(
742 : TGridFunction& u,
743 : MatrixOperator<typename TGridFunction::algebra_type::matrix_type,
744 : typename TGridFunction::vector_type>& A,
745 : const char* filename) {
746 : PROFILE_FUNC_GROUP("debug");
747 : // forward
748 0 : WriteMatrixToConnectionViewer(filename, A.get_matrix(), u);
749 0 : }
750 :
751 : /**
752 : * \brief Save the assembled matrix of a matrix operator to MatrixMarket format
753 : *
754 : * \param[in] filename name of the file; must end on '<tt>.mtx</tt>'
755 : * \param[in] A matrix operator of with \c CPUAlgebra matrix
756 : * \param[in] comment optional comment for the header of the MTX file
757 : *
758 : * \note Until now only CPUAlgebra matrices are supported.
759 : */
760 : // TODO extend MatrixIO to support other than CPUAlgebra
761 0 : inline void SaveMatrixToMTX( const char *filename,
762 : MatrixOperator< CPUAlgebra::matrix_type, CPUAlgebra::vector_type >& A,
763 : std::string comment="%Generated with ug4." ) {
764 : PROFILE_FUNC_GROUP("debug");
765 :
766 : // check name
767 0 : if ( !FileTypeIs( filename, ".mtx" ) ) {
768 0 : UG_THROW( "Please use '.mtx' as file extension for MatrixMarket exchange files."
769 : " (Filename is '" << filename << "')" );
770 : }
771 :
772 : // automatically detect the file mode to use
773 0 : MatrixIO::OpenMode openMode = (FileExists( filename )) ? MatrixIO::EXISTING : MatrixIO::NEW;
774 :
775 : // create MatrixIO object for handling the export
776 0 : MatrixIOMtx mtx( filename, openMode );
777 0 : mtx.write_from( A.get_matrix(), comment );
778 0 : }
779 :
780 : template<class TFunction>
781 0 : void WriteVectorToConnectionViewer(const char *filename,
782 : const typename TFunction::algebra_type::vector_type &b,
783 : const TFunction &u,
784 : const typename TFunction::algebra_type::vector_type *pCompareVec = NULL) {
785 : PROFILE_FUNC_GROUP("debug");
786 : // check name
787 0 : if ( !FileTypeIs( filename, ".vec") ) {
788 0 : UG_THROW( "Only '.vec' format supported for vectors, but"
789 : " filename is '" << filename << "'." );
790 : }
791 :
792 : // get positions of vertices
793 : const static int dim = TFunction::domain_type::dim;
794 : std::vector<MathVector<dim> > vPos;
795 0 : ExtractPositions(u, vPos);
796 :
797 : // write vector
798 0 : ConnectionViewer::WriteVectorPar( filename, b, &vPos[0], dim, pCompareVec);
799 0 : WriteAlgebraIndices(filename, u.domain(),u.dof_distribution());
800 0 : }
801 :
802 : template<class TFunction>
803 0 : void WriteVectorToConnectionViewer(
804 : const char *filename,
805 : const typename TFunction::algebra_type::matrix_type &A,
806 : const typename TFunction::algebra_type::vector_type &b,
807 : const TFunction &u,
808 : const typename TFunction::algebra_type::vector_type *pCompareVec = NULL) {
809 : PROFILE_FUNC_GROUP("debug");
810 : // get dimension
811 : const static int dim = TFunction::domain_type::dim;
812 :
813 : // check name
814 0 : if ( !FileTypeIs( filename, ".vec") ) {
815 0 : UG_THROW( "Only '.vec' format supported for vectors." );
816 : }
817 :
818 : // get positions of vertices
819 : std::vector<MathVector<dim> > positions;
820 0 : ExtractPositions(u, positions);
821 :
822 : // write vector
823 0 : ConnectionViewer::WriteVectorPar( filename, A, b, &positions[0], dim, pCompareVec );
824 0 : }
825 :
826 : template<typename TGridFunction>
827 0 : void SaveVectorForConnectionViewer(TGridFunction& b, const char* filename) {
828 0 : WriteVectorToConnectionViewer(filename, b, b);
829 0 : }
830 :
831 : template<typename TGridFunction>
832 0 : void SaveVectorDiffForConnectionViewer(TGridFunction& b, TGridFunction& bCompare, const char* filename) {
833 0 : WriteVectorToConnectionViewer(filename, b, b, &bCompare);
834 0 : }
835 :
836 : template<typename TGridFunction>
837 0 : void SaveVectorForConnectionViewer(
838 : TGridFunction& u,
839 : MatrixOperator<typename TGridFunction::algebra_type::matrix_type,
840 : typename TGridFunction::vector_type>& A,
841 : const char* filename) {
842 0 : WriteVectorToConnectionViewer(filename, A.get_matrix(), u, u);
843 0 : }
844 :
845 : template<typename TGridFunction>
846 0 : void SaveVectorDiffForConnectionViewer(
847 : TGridFunction& u,
848 : TGridFunction& compareVec,
849 : MatrixOperator<typename TGridFunction::algebra_type::matrix_type,
850 : typename TGridFunction::vector_type>& A,
851 : const char* filename) {
852 : // forward
853 0 : WriteVectorToConnectionViewer(filename, A.get_matrix(), u, u, &compareVec);
854 0 : }
855 :
856 : // from connection_viewer_input.h
857 : // with additional checks
858 : template<typename vector_type>
859 0 : bool ReadVector(std::string filename, vector_type &vec,int dim)
860 : {
861 0 : Progress p;
862 0 : std::cout << " Reading std::vector from " << filename << "... ";
863 0 : std::fstream matfile(filename.c_str(), std::ios::in);
864 0 : if(matfile.is_open() == false) { std::cout << "failed.\n"; return false; }
865 :
866 0 : int version=-1, dimension=-1, gridsize;
867 :
868 0 : matfile >> version;
869 0 : matfile >> dimension;
870 0 : matfile >> gridsize;
871 :
872 : assert(version == 1);
873 : assert(dimension == dim);
874 : // todo check positions and not just size
875 : assert(gridsize == (int)vec.size());
876 :
877 :
878 0 : PROGRESS_START(prog, gridsize*2, "ReadVector " << dimension << "d from " << filename << " , " << gridsize << " x " << gridsize);
879 0 : for(int i=0; i<gridsize; i++)
880 : {
881 0 : if(i%100) { PROGRESS_UPDATE(prog, i); }
882 0 : if(matfile.eof())
883 : {
884 0 : std::cout << " failed.\n";
885 : assert(0);
886 0 : return false;
887 : }
888 : double x, y, z;
889 : matfile >> x >> y;
890 0 : if(dimension==3) matfile >> z;
891 : }
892 :
893 : int printStringsInWindow;
894 0 : matfile >> printStringsInWindow;
895 :
896 : // vec.resize(gridsize);
897 : bool bEOF = matfile.eof();
898 0 : while(!bEOF)
899 : {
900 : int from, to; double value;
901 0 : char c = matfile.peek();
902 0 : if(c == -1 || c == 'c' || c == 'v' || matfile.eof())
903 : break;
904 :
905 0 : matfile >> from >> to >> value;
906 : assert(from == to);
907 0 : vec[from] = value;
908 0 : if(from%100) { PROGRESS_UPDATE(prog, from); }
909 : bEOF = matfile.eof();
910 : }
911 : return true;
912 0 : }
913 :
914 : // load vector that has been saved in connection viewer format and write it
915 : // into grid function
916 : template<typename TGridFunction>
917 0 : void LoadVector(TGridFunction& u,const char* filename){
918 : PROFILE_FUNC();
919 : typename TGridFunction::algebra_type::vector_type b;
920 : b.resize(u.num_indices());
921 0 : ReadVector(filename,b,TGridFunction::dim);
922 0 : u.assign(b);
923 0 : }
924 :
925 : // Same as before, but for comma separated value (CSV)
926 : template<class TFunction>
927 0 : void WriteVectorCSV(const char *filename,
928 : const typename TFunction::algebra_type::vector_type &b,
929 : const TFunction &u) {
930 : PROFILE_FUNC_GROUP("debug");
931 : // get dimension
932 : const static int dim = TFunction::domain_type::dim;
933 :
934 : // check name
935 0 : if ( !FileTypeIs( filename, ".csv") ) {
936 0 : UG_THROW( "Only '.csv' format supported for vectors, but"
937 : " filename is '" << filename << "'." );
938 : }
939 :
940 : // extended filename
941 : // add p000X extension in parallel
942 : #ifdef UG_PARALLEL
943 : std::string name(filename);
944 : size_t iExtPos = name.find_last_of(".");
945 : name.resize(iExtPos);
946 : int rank = pcl::ProcRank();
947 : char ext[20]; snprintf(ext, 20, "_p%05d.csv", rank);
948 : name.append(ext);
949 : #endif
950 :
951 : // get positions of vertices
952 : std::vector<MathVector<dim> > positions;
953 0 : ExtractPositions(u, positions);
954 :
955 : // write vector
956 0 : WriteVectorCSV( filename, b, &positions[0], dim );
957 0 : }
958 :
959 : template<typename TGridFunction>
960 0 : void SaveVectorCSV(TGridFunction& b, const char* filename) {
961 : PROFILE_FUNC();
962 0 : WriteVectorCSV(filename, b, b);
963 0 : }
964 :
965 : /**
966 : * \brief Calculates the average of the pointwise difference of two functions on given subset
967 : * \details Iterates over all vertices of given subset, calculates difference
968 : * of \c fct1 and \c fct2 and computes arithmetic mean of the differences at
969 : * all vertices (grid points).
970 : * \note \c fct1 and \c fct2 must both be defined on \c subset !
971 : * \param[in] spGridFct GridFunction holding functions \c fct1 and \c fct2 defined on subset \c subset
972 : * \param[in] subset name of the subset to compare on
973 : * \param[in] fct1 name of the first function
974 : * \param[in] fct2 name of the second function
975 : * \return arithmetic average of the pointwise differences of \c fct1 and \c fct2 on
976 : * all vertices of \c subset
977 : */
978 : template<typename TDomain, typename TAlgebra>
979 0 : number AverageFunctionDifference(
980 : SmartPtr< GridFunction<TDomain, TAlgebra> > spGridFct,
981 : std::string subset, std::string fct1, std::string fct2 )
982 : {
983 : PROFILE_FUNC();
984 : // get subset index
985 : size_t subSetID = spGridFct->subset_id_by_name( subset.c_str() );
986 :
987 : // get function indices
988 : size_t fct1ID = spGridFct->fct_id_by_name( fct1.c_str() );
989 : size_t fct2ID = spGridFct->fct_id_by_name( fct2.c_str() );
990 :
991 : // create space for sum of difference
992 : number sum = 0.0;
993 : size_t numElements = 0;
994 :
995 : // loop over all vertices in given subset and compare values of fct1 and fct2
996 : typedef typename GridFunction<TDomain, TAlgebra>::template traits<Vertex>::const_iterator gridFctIterator;
997 0 : for( gridFctIterator iter = spGridFct->template begin<Vertex>((int)subSetID);
998 0 : iter != spGridFct->template end<Vertex>((int)subSetID); ++iter ) {
999 : // get dof_indices for the two functions on given subset
1000 : std::vector< DoFIndex > indFct1, indFct2;
1001 : spGridFct->template dof_indices<Vertex>( *iter, fct1ID, indFct1 );
1002 : spGridFct->template dof_indices<Vertex>( *iter, fct2ID, indFct2 );
1003 :
1004 : // calculate the difference between the two functions at this grid point
1005 0 : sum += DoFRef( *spGridFct, indFct1[0] ) - DoFRef( *spGridFct, indFct2[0] );
1006 0 : numElements++;
1007 : }
1008 :
1009 : // return overal arithmetic average of differences
1010 0 : return sum / numElements;
1011 : }
1012 :
1013 :
1014 : ////////////////////////////////////////////////////////////////////////////////
1015 : // Grid Function Debug Writer
1016 : ////////////////////////////////////////////////////////////////////////////////
1017 :
1018 : template<typename TDomain, typename TAlgebra>
1019 : class GridFunctionDebugWriter: public IDebugWriter<TAlgebra>
1020 : {
1021 : /// dimension
1022 : static const int dim = TDomain::dim;
1023 :
1024 : public:
1025 : /// type of matrix
1026 : typedef TAlgebra algebra_type;
1027 :
1028 : /// type of vector
1029 : typedef typename algebra_type::vector_type vector_type;
1030 :
1031 : /// type of matrix
1032 : typedef typename algebra_type::matrix_type matrix_type;
1033 :
1034 : /// type of approximation space
1035 : typedef ApproximationSpace<TDomain> approximation_space_type;
1036 :
1037 : typedef IDebugWriter<TAlgebra> super;
1038 : using super::get_base_dir;
1039 :
1040 : public:
1041 : /// Constructor
1042 0 : GridFunctionDebugWriter(
1043 : SmartPtr<ApproximationSpace<TDomain> > spApproxSpace) :
1044 0 : m_spApproxSpace(spApproxSpace), bConnViewerOut(
1045 0 : true), bConnViewerIndices(false), bVTKOut(true), m_printConsistent(true)
1046 : {
1047 0 : IVectorDebugWriter<vector_type>::m_currentDim = dim;
1048 : reset();
1049 0 : }
1050 :
1051 0 : virtual ~GridFunctionDebugWriter() {};
1052 :
1053 : /// sets the grid level
1054 0 : void set_grid_level(const GridLevel& gridLevel) {
1055 0 : m_glFrom = gridLevel; m_glTo = gridLevel;
1056 0 : }
1057 :
1058 : /// sets the grid level
1059 : void set_grid_levels(const GridLevel& glFrom, const GridLevel& glTo) {
1060 0 : m_glFrom = glFrom; m_glTo = glTo;
1061 : }
1062 :
1063 : /// returns current grid level
1064 0 : GridLevel grid_level() const {return m_glFrom;}
1065 :
1066 : /// sets to toplevel on surface
1067 0 : void reset() {
1068 : set_grid_level(GridLevel(GridLevel::TOP, GridLevel::SURFACE));
1069 0 : }
1070 :
1071 :
1072 : /// sets if writing to vtk
1073 0 : void set_vtk_output(bool b) {bVTKOut = b;}
1074 :
1075 : /// sets if writing to conn viewer
1076 0 : void set_conn_viewer_output(bool b) {bConnViewerOut = b;}
1077 :
1078 : /// sets if .indices file is written or conn viewer
1079 0 : void set_conn_viewer_indices(bool b) { bConnViewerIndices = b; }
1080 :
1081 : /// sets if data shall be made consistent before printing
1082 0 : void set_print_consistent(bool b) {m_printConsistent = b;}
1083 :
1084 : /// write vector
1085 0 : virtual void write_vector(const vector_type& vec, const char* filename) {
1086 : // write to conn viewer
1087 0 : if (bConnViewerOut)
1088 0 : write_vector_to_conn_viewer(vec, filename);
1089 :
1090 : // write to vtk
1091 0 : if (bVTKOut)
1092 0 : write_vector_to_vtk(vec, filename);
1093 0 : }
1094 :
1095 0 : virtual void update_positions()
1096 : {
1097 0 : extract_positions(m_glFrom);
1098 0 : }
1099 :
1100 : /// write matrix
1101 0 : virtual void write_matrix(const matrix_type& mat, const char* filename) {
1102 : PROFILE_FUNC_GROUP("debug");
1103 : // something to do ?
1104 0 : if (!bConnViewerOut)
1105 0 : return;
1106 :
1107 0 : update_positions();
1108 :
1109 : std::string name;
1110 0 : this->compose_file_path (name);
1111 : name += filename;
1112 :
1113 0 : if ( !FileTypeIs( filename, ".mat" ) ) {
1114 0 : UG_THROW( "Only '.mat' format supported for matrices, but"
1115 : " filename is '" << filename << "'." );
1116 : }
1117 :
1118 : // write to file
1119 0 : write_algebra_indices_CV(name);
1120 : if(m_glFrom == m_glTo){
1121 0 : if(mat.num_rows() != mat.num_cols())
1122 0 : UG_THROW("DebugWriter: grid level the same, but non-square matrix.");
1123 :
1124 0 : const std::vector<MathVector<dim> >& vPos = this->template get_positions<dim>();
1125 0 : if(vPos.empty())
1126 0 : ConnectionViewer::WriteMatrixPar(name.c_str(), mat,(MathVector<dim>*)NULL, dim);
1127 : else
1128 0 : ConnectionViewer::WriteMatrixPar(name.c_str(), mat, &vPos[0], dim);
1129 : }
1130 : else{
1131 0 : const std::vector<MathVector<dim> >& vFromPos = this->template get_positions<dim>();
1132 :
1133 : std::vector<MathVector<dim> > vToPos;
1134 0 : ExtractPositions<TDomain>(m_spApproxSpace->domain(),
1135 0 : m_spApproxSpace->dof_distribution(m_glTo),
1136 : vToPos);
1137 :
1138 0 : if(mat.num_cols() == vFromPos.size() && mat.num_rows() == vToPos.size())
1139 : {
1140 0 : ConnectionViewer::WriteMatrixPar(name, mat, vFromPos, vToPos, dim);
1141 : }
1142 : else{
1143 0 : UG_THROW("GridFunctionDebugWriter: Wrong size of matrix for writing"
1144 : "matrix ("<<m_glTo<<" x "<<m_glFrom<<"), that would be a ("
1145 : <<vToPos.size()<<" x "<<vFromPos.size()<<") matrix. But "
1146 : "passed matrix of size ("<<mat.num_rows()<<" x "
1147 : <<mat.num_cols()<<").");
1148 : }
1149 0 : }
1150 : }
1151 :
1152 : protected:
1153 :
1154 0 : void write_algebra_indices_CV(std::string name)
1155 : {
1156 0 : if(bConnViewerIndices)
1157 0 : WriteAlgebraIndices<TDomain>(name, m_spApproxSpace->domain(), m_spApproxSpace->dof_distribution(m_glFrom));
1158 0 : }
1159 :
1160 : /// reads the positions
1161 0 : void extract_positions(const GridLevel& gridLevel) {
1162 : // extract positions for this grid function
1163 : std::vector<MathVector<dim> >& vPos =
1164 : this->template positions<dim>();
1165 :
1166 : vPos.clear();
1167 0 : ExtractPositions<TDomain>(
1168 : m_spApproxSpace->domain(),
1169 : m_spApproxSpace->dof_distribution(gridLevel),
1170 : vPos);
1171 0 : }
1172 :
1173 : /// write vector
1174 0 : virtual void write_vector_to_conn_viewer(const vector_type& vec,
1175 : const char* filename)
1176 : {
1177 : PROFILE_FUNC_GROUP("debug");
1178 0 : update_positions();
1179 :
1180 : std::string name;
1181 0 : this->compose_file_path (name);
1182 : name += filename;
1183 :
1184 0 : if ( !FileTypeIs( filename, ".vec" ) ) {
1185 0 : UG_THROW( "Only '.vec' format supported for vectors, but"
1186 : " filename is '" << name << "'.");
1187 : }
1188 :
1189 : // write
1190 0 : write_algebra_indices_CV(filename);
1191 : const std::vector<MathVector<dim> >& vPos =
1192 0 : this->template get_positions<dim>();
1193 0 : if(vPos.empty())
1194 0 : ConnectionViewer::WriteVectorPar(name.c_str(), vec, (MathVector<dim>*)NULL, dim);
1195 : else
1196 0 : ConnectionViewer::WriteVectorPar(name.c_str(), vec, &vPos[0], dim);
1197 0 : }
1198 :
1199 0 : void write_vector_to_vtk(const vector_type& vec, const char* filename) {
1200 : PROFILE_FUNC_GROUP("debug");
1201 : // check name
1202 :
1203 : std::string name;
1204 0 : this->compose_file_path (name);
1205 : name += filename;
1206 :
1207 : typedef GridFunction<TDomain, TAlgebra> TGridFunction;
1208 0 : TGridFunction vtkFunc(
1209 : m_spApproxSpace,
1210 0 : m_spApproxSpace->dof_distribution(m_glFrom));
1211 0 : vtkFunc.resize_values(vec.size());
1212 0 : vtkFunc.assign(vec);
1213 0 : VTKOutput<dim> out;
1214 0 : out.print(name.c_str(), vtkFunc, m_printConsistent);
1215 0 : }
1216 :
1217 : protected:
1218 : // underlying approximation space
1219 : SmartPtr<approximation_space_type> m_spApproxSpace;
1220 :
1221 : // flag if write to conn viewer
1222 : bool bConnViewerOut;
1223 : bool bConnViewerIndices;
1224 :
1225 : // flag if write to vtk
1226 : bool bVTKOut;
1227 :
1228 : // flag if data shall be made consistent before printing
1229 : bool m_printConsistent;
1230 :
1231 : // current grid level
1232 : GridLevel m_glFrom, m_glTo;
1233 : };
1234 :
1235 : ////////////////////////////////////////////////////////////////////////////////
1236 : // Grid Function Position Provider
1237 : ////////////////////////////////////////////////////////////////////////////////
1238 :
1239 : template<typename TGridFunction>
1240 : class GridFunctionPositionProvider: public IPositionProvider<
1241 : TGridFunction::domain_type::dim> {
1242 : public:
1243 : /// Constructor
1244 0 : GridFunctionPositionProvider() :
1245 0 : m_pGridFunc(NULL) {
1246 : }
1247 0 : GridFunctionPositionProvider(const TGridFunction& u) :
1248 0 : m_pGridFunc(&u) {
1249 : }
1250 0 : void set_reference_grid_function(const TGridFunction& u) {
1251 0 : m_pGridFunc = &u;
1252 0 : }
1253 :
1254 0 : virtual bool get_positions(
1255 : std::vector<MathVector<TGridFunction::domain_type::dim> >&vec) {
1256 : UG_ASSERT(m_pGridFunc != NULL,
1257 : "provide a grid function with set_reference_grid_function");
1258 0 : ExtractPositions(*m_pGridFunc, vec);
1259 0 : return true;
1260 : }
1261 :
1262 : protected:
1263 : const TGridFunction *m_pGridFunc;
1264 : };
1265 :
1266 : ////////////////////////////////////////////////////////////////////////////////
1267 : // Grid Function Vector Writer
1268 : ////////////////////////////////////////////////////////////////////////////////
1269 :
1270 : template<typename TGridFunction, typename TVector>
1271 : class GridFunctionVectorWriter: public IVectorWriter<TVector> {
1272 : public:
1273 : typedef typename TVector::value_type value_type;
1274 : typedef typename TGridFunction::domain_type domain_type;
1275 : typedef TVector vector_type;
1276 :
1277 : public:
1278 : /// Constructor
1279 0 : GridFunctionVectorWriter() :
1280 0 : m_pGridFunc(NULL) {
1281 : }
1282 :
1283 0 : void set_user_data(SmartPtr<UserData<number, domain_type::dim> > userData) {
1284 0 : m_userData = userData;
1285 0 : }
1286 :
1287 0 : void set_reference_grid_function(const TGridFunction& u) {
1288 0 : m_pGridFunc = &u;
1289 0 : }
1290 :
1291 : /*virtual double calculate(MathVector<3> &v, double time)
1292 : {
1293 : }*/
1294 :
1295 0 : virtual bool update(vector_type &vec) {
1296 : PROFILE_FUNC_GROUP("debug");
1297 : UG_ASSERT(m_pGridFunc != NULL,
1298 : "provide a grid function with set_reference_grid_function");
1299 : UG_ASSERT(m_userData.valid(), "provide user data with set_user_data");
1300 :
1301 0 : const TGridFunction &u = *m_pGridFunc;
1302 :
1303 : // get position accessor
1304 :
1305 : const typename domain_type::position_accessor_type& aaPos =
1306 0 : u.domain()->position_accessor();
1307 :
1308 : // number of total dofs
1309 0 : int nr = u.num_indices();
1310 :
1311 : // resize positions
1312 0 : vec.resize(nr);
1313 :
1314 : typedef typename TGridFunction::template traits<Vertex>::const_iterator const_iterator;
1315 :
1316 : // loop all subsets
1317 0 : for (int si = 0; si < u.num_subsets(); ++si) {
1318 : // loop all vertices
1319 : const_iterator iter = u.template begin<Vertex>(si);
1320 : const_iterator iterEnd = u.template end<Vertex>(si);
1321 :
1322 0 : for (; iter != iterEnd; ++iter) {
1323 : // get vertex
1324 : Vertex* v = *iter;
1325 :
1326 : // algebra indices vector
1327 : std::vector<size_t> ind;
1328 :
1329 : // load indices associated with vertex
1330 : u.inner_algebra_indices(v, ind);
1331 :
1332 : number t = 0.0;
1333 :
1334 : number d;
1335 0 : (*m_userData)(d, aaPos[v], t, si);
1336 :
1337 : // write
1338 0 : for (size_t i = 0; i < ind.size(); ++i) {
1339 0 : const size_t index = ind[i];
1340 0 : for (size_t alpha = 0; alpha < GetSize(vec[index]); ++alpha)
1341 0 : BlockRef(vec[index], alpha) = d;
1342 : }
1343 : }
1344 : }
1345 0 : return true;
1346 : }
1347 :
1348 : protected:
1349 : const TGridFunction *m_pGridFunc;
1350 : SmartPtr<UserData<number, domain_type::dim> > m_userData;
1351 : };
1352 :
1353 : ////////////////////////////////////////////////////////////////////////////////
1354 : // Grid Function Vector Writer Dirichlet 0
1355 : ////////////////////////////////////////////////////////////////////////////////
1356 :
1357 : template<typename TGridFunction>
1358 : class GridFunctionVectorWriterDirichlet0: public IVectorWriter<
1359 : typename TGridFunction::algebra_type::vector_type> {
1360 : public:
1361 : typedef typename TGridFunction::domain_type domain_type;
1362 :
1363 : typedef typename TGridFunction::approximation_space_type approximation_space_type;
1364 :
1365 : typedef typename TGridFunction::algebra_type algebra_type;
1366 : typedef typename algebra_type::vector_type vector_type;
1367 : typedef typename vector_type::value_type value_type;
1368 :
1369 : public:
1370 : /// Constructor
1371 0 : GridFunctionVectorWriterDirichlet0() :
1372 0 : m_pApproxSpace(NULL), m_spPostProcess(NULL), m_level(-1) {
1373 : }
1374 :
1375 0 : void set_level(size_t level) {
1376 0 : m_level = level;
1377 0 : }
1378 :
1379 0 : void init(SmartPtr<IConstraint<algebra_type> > pp,
1380 : approximation_space_type& approxSpace) {
1381 0 : m_spPostProcess = pp;
1382 0 : m_pApproxSpace = &approxSpace;
1383 0 : }
1384 :
1385 : /*virtual double calculate(MathVector<3> &v, double time)
1386 : {
1387 : }*/
1388 :
1389 0 : virtual bool update(vector_type &vec) {
1390 : PROFILE_FUNC_GROUP("debug");
1391 : UG_ASSERT(m_spPostProcess.valid(), "provide a post process with init");
1392 : UG_ASSERT(m_pApproxSpace != NULL, "provide approximation space init");
1393 :
1394 : size_t numDoFs = 0;
1395 : GridLevel gl;
1396 0 : if (m_level == (size_t) -1) {
1397 0 : numDoFs = m_pApproxSpace->dof_distribution(GridLevel(GridLevel::TOP, GridLevel::SURFACE))->num_indices();
1398 0 : gl = GridLevel();
1399 : } else {
1400 : numDoFs =
1401 0 : m_pApproxSpace->dof_distribution(GridLevel(m_level, GridLevel::LEVEL, true))->num_indices();
1402 0 : gl = GridLevel(m_level, GridLevel::LEVEL);
1403 : }
1404 : vec.resize(numDoFs);
1405 : vec.set(1.0);
1406 :
1407 0 : m_spPostProcess->adjust_defect(vec, vec, m_pApproxSpace->dof_distribution(gl), CT_DIRICHLET);
1408 :
1409 0 : return true;
1410 : }
1411 :
1412 : protected:
1413 : approximation_space_type * m_pApproxSpace;
1414 : SmartPtr<IConstraint<algebra_type> > m_spPostProcess;
1415 : size_t m_level;
1416 : };
1417 :
1418 : } // end namespace ug
1419 :
1420 : #endif /* __H__UG__LIB_DISC__FUNCTION_SPACE__GRID_FUNCTION_UTIL__ */
|
