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