Line data Source code
1 : /*
2 : * Copyright (c) 2013-2015: G-CSC, Goethe University Frankfurt
3 : * Author: Ivo Muha
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 : // extern headers
34 : #include <iostream>
35 : #include <sstream>
36 : #include <string>
37 : #include <limits>
38 : #include <set>
39 :
40 : // include bridge
41 : #include "bridge/bridge.h"
42 : #include "bridge/util.h"
43 : #include "bridge/util_domain_algebra_dependent.h"
44 :
45 : // lib_disc includes
46 : #include "lib_disc/dof_manager/dof_distribution.h"
47 : #include "lib_disc/function_spaces/grid_function.h"
48 :
49 : // user data
50 : #include "lib_disc/spatial_disc/user_data/user_data.h"
51 :
52 : #include "lib_disc/reference_element/reference_element_traits.h"
53 : #include "lib_disc/reference_element/reference_mapping_provider.h"
54 : #include "lib_grid/algorithms/space_partitioning/lg_ntree.h"
55 : #include "common/util/provider.h"
56 : #include "lib_disc/domain_util.h"
57 : #include "lib_disc/time_disc/time_integrator_observers/time_integrator_observer_interface.h"
58 :
59 : // pcl includes
60 : #ifdef UG_PARALLEL
61 : #include "pcl/pcl_util.h"
62 : #include "pcl/pcl_process_communicator.h"
63 : #endif
64 :
65 : #ifdef UG_PARALLEL
66 : #include "lib_grid/parallelization/distributed_grid.h"
67 : #endif
68 :
69 : using namespace std;
70 :
71 : namespace ug{
72 : namespace bridge{
73 : namespace Evaluate{
74 :
75 :
76 : template <typename TDomain, typename TAlgebra>
77 0 : class NumberValuedUserDataEvaluator
78 : {
79 : static const int dim = TDomain::dim;
80 : typedef GridFunction<TDomain, TAlgebra> TGridFunction;
81 : typedef typename TGridFunction::template dim_traits<dim>::grid_base_object TElem;
82 :
83 : typedef lg_ntree<dim, dim, TElem> tree_t;
84 :
85 : public:
86 :
87 0 : NumberValuedUserDataEvaluator(SmartPtr<UserData<number, dim> > userData) : m_userData(userData)
88 : {
89 :
90 : }
91 :
92 0 : number evaluateLua(const std::vector<number>& pos, SmartPtr<TGridFunction> u, number time)
93 : {
94 0 : number result = 0;
95 :
96 0 : this->evaluate(pos, result, u, time);
97 :
98 0 : return result;
99 : }
100 :
101 :
102 0 : bool evaluate(const std::vector<number>& pos,
103 : number& result,
104 : SmartPtr<TGridFunction> u,
105 : number time)
106 : {
107 0 : if(!m_initialized)
108 0 : initialize(u);
109 :
110 0 : bool found = evaluateOnThisProcess(pos, result, u, time);
111 :
112 :
113 : #ifndef UG_PARALLEL
114 0 : return found;
115 : #endif
116 :
117 : #ifdef UG_PARALLEL
118 : // share value between all procs
119 : pcl::ProcessCommunicator com;
120 : int numFound = (found ? 1 : 0);
121 : numFound = com.allreduce(numFound, PCL_RO_SUM);
122 :
123 : if(numFound == 0)
124 : return false;
125 :
126 : // get overall value
127 : // if found on more than one processor, the data will be the same
128 : number globalResult = com.allreduce(result, PCL_RO_SUM);
129 :
130 : // set as result
131 : result = globalResult / numFound;
132 :
133 :
134 : return true;
135 : #endif
136 : }
137 :
138 :
139 : private:
140 :
141 0 : bool evaluateOnThisProcess(const std::vector<number>& pos,
142 : number& result,
143 : SmartPtr<TGridFunction> u,
144 : number time)
145 : {
146 0 : TElem* elem = NULL;
147 :
148 : MathVector<dim> globalPosition;
149 :
150 0 : for(int i = 0; i < dim; i++)
151 : {
152 0 : globalPosition[i] = pos[i];
153 : }
154 :
155 0 : if(!FindContainingElement(elem, *m_tree, globalPosition))
156 : {
157 0 : result = 0;
158 0 : return false;
159 : }
160 :
161 : // get corners of element
162 : std::vector<MathVector<dim> > vCornerCoords;
163 0 : CollectCornerCoordinates(vCornerCoords, *elem, *u->domain());
164 :
165 : // get subset
166 0 : int si = u->domain()->subset_handler()->get_subset_index(elem);
167 :
168 : // reference object id
169 0 : const ReferenceObjectID roid = elem->reference_object_id();
170 :
171 : // get local position of DoF
172 : DimReferenceMapping<dim, dim>& map
173 : = ReferenceMappingProvider::get<dim, dim>(roid, vCornerCoords);
174 : MathVector<dim> locPos;
175 : VecSet(locPos, 0.5);
176 0 : map.global_to_local(locPos, globalPosition);
177 :
178 : // storage for the result
179 : number value;
180 :
181 : // get local solution if needed
182 0 : if(m_userData->requires_grid_fct())
183 : {
184 : // create storage
185 : LocalIndices ind;
186 : LocalVector localU;
187 :
188 : // get global indices
189 : u->indices(elem, ind);
190 :
191 : // adapt local algebra
192 0 : localU.resize(ind);
193 :
194 : // read local values of u
195 0 : GetLocalVector(localU, *u);
196 :
197 : try
198 : {
199 0 : (*m_userData)(&value, &globalPosition, time, si, elem,
200 : &vCornerCoords[0], &locPos, 1, &localU, NULL);
201 : }
202 0 : UG_CATCH_THROW("NumberValuedUserDataEvaluator: Cannot evaluate data.");
203 : }
204 : else
205 : {
206 : try
207 : {
208 0 : (*m_userData)(&value, &globalPosition, time, si, elem,
209 : &vCornerCoords[0], &locPos, 1, NULL, NULL);
210 : }
211 0 : UG_CATCH_THROW("NumberValuedUserDataEvaluator: Cannot evaluate data.");
212 : }
213 :
214 0 : result = value;
215 :
216 : return true;
217 0 : }
218 :
219 0 : void initialize(SmartPtr<TGridFunction> u)
220 : {
221 0 : m_initialized = true;
222 :
223 0 : m_tree = make_sp(new tree_t(*u->domain()->grid(), u->domain()->position_attachment()));
224 0 : m_tree->create_tree(u->template begin<TElem>(), u->template end<TElem>());
225 0 : }
226 :
227 : bool m_initialized = false;
228 : SmartPtr<tree_t> m_tree;
229 : SmartPtr<UserData<number, dim> > m_userData;
230 :
231 : };
232 :
233 : template <typename TDomain, typename TAlgebra>
234 0 : class VectorValuedUserDataEvaluator
235 : {
236 : static const int dim = TDomain::dim;
237 : typedef GridFunction<TDomain, TAlgebra> TGridFunction;
238 : typedef typename TGridFunction::template dim_traits<dim>::grid_base_object TElem;
239 :
240 : typedef lg_ntree<dim, dim, TElem> tree_t;
241 :
242 : public:
243 :
244 0 : VectorValuedUserDataEvaluator(SmartPtr<UserData<MathVector<dim>, dim> > userData) : m_userData(userData)
245 : {
246 :
247 : }
248 :
249 0 : std::vector<number> evaluateLua(const std::vector<number>& pos, SmartPtr<TGridFunction> u, number time, bool interpolateOverNeighbouringTriangles)
250 : {
251 : std::vector<number> result;
252 :
253 0 : this->evaluate(pos, result, u, time, interpolateOverNeighbouringTriangles);
254 :
255 0 : return result;
256 0 : }
257 :
258 :
259 0 : bool evaluate(const std::vector<number>& pos,
260 : std::vector<number>& result,
261 : SmartPtr<TGridFunction> u,
262 : number time,
263 : bool interpolateOverNeighbouringTriangles)
264 : {
265 0 : if(!m_initialized)
266 0 : initialize(u);
267 :
268 : bool found = false;
269 :
270 0 : if(interpolateOverNeighbouringTriangles)
271 0 : found = evaluateOnThisProcessNeighbouring(pos, result, u, time);
272 : else
273 0 : found = evaluateOnThisProcess(pos, result, u, time);
274 :
275 :
276 : #ifndef UG_PARALLEL
277 0 : return found;
278 : #else
279 : if (! found)
280 : for (int i = 0; i < dim; i++)
281 : result[i] = 0;
282 : #endif
283 :
284 : #ifdef UG_PARALLEL
285 : // share value between all procs
286 : pcl::ProcessCommunicator com;
287 : int numFound = (found ? 1 : 0);
288 : numFound = com.allreduce(numFound, PCL_RO_SUM);
289 :
290 : if(numFound == 0)
291 : return false;
292 :
293 : // get overall value
294 : // if found on more than one processor, the data will be the same
295 : std::vector<number> globalResult;
296 : com.allreduce(result, globalResult, PCL_RO_SUM);
297 :
298 : // set as result
299 : for(int i = 0; i < dim; i++)
300 : {
301 : result[i] = globalResult[i] / numFound;
302 : }
303 :
304 : return true;
305 : #endif
306 : }
307 :
308 :
309 : private:
310 :
311 0 : bool evaluateOnThisProcessNeighbouring(const std::vector<number>& pos,
312 : std::vector<number>& result,
313 : SmartPtr<TGridFunction> u,
314 : number time)
315 : {
316 : typedef typename TDomain::grid_type TGrid;
317 : typedef typename Grid::traits<TElem>::secure_container TElemContainer;
318 :
319 0 : result.resize(TDomain::dim);
320 :
321 0 : TElem* elem = NULL;
322 :
323 : MathVector<dim> globalPosition;
324 :
325 0 : for(int i = 0; i < dim; i++)
326 : {
327 0 : globalPosition[i] = pos[i];
328 : }
329 :
330 0 : if(!FindContainingElement(elem, *m_tree, globalPosition))
331 : {
332 0 : return false;
333 : }
334 :
335 : // find all adjacent faces/volumes/edges
336 : std::set<TElem* > elements;
337 : elements.insert(elem);
338 :
339 0 : TGrid* grid = u->domain()->grid().get();
340 :
341 0 : for(size_t i = 0; i < elem->num_vertices(); i++)
342 : {
343 : // for each vertex of the found element: get associated faces/volumes/edges
344 0 : Vertex* vrt = elem->vertex(i);
345 :
346 : TElemContainer associatedElements;
347 0 : grid->associated_elements(associatedElements, vrt);
348 :
349 0 : for(size_t j = 0; j < associatedElements.size(); j++)
350 : {
351 0 : elements.insert(associatedElements[j]);
352 : }
353 : }
354 :
355 : double sumWeights = 0;
356 :
357 0 : for(TElem * element : elements)
358 : {
359 : // get corners of element
360 : std::vector<MathVector<dim> > vCornerCoords;
361 0 : CollectCornerCoordinates(vCornerCoords, *element, *u->domain());
362 :
363 : // find center of element
364 : MathVector<dim> centerGlobal;
365 0 : for(size_t i = 0; i < vCornerCoords.size(); i++)
366 : {
367 0 : for(int d = 0; d < dim; d++)
368 0 : centerGlobal[d] += vCornerCoords[i][d] / vCornerCoords.size();
369 : }
370 :
371 : // calculate distance to evaluated position
372 : double distance = 0;
373 0 : for(int d = 0; d < dim; d++)
374 0 : distance += (pos[d]-centerGlobal[d])*(pos[d]-centerGlobal[d]);
375 0 : distance = sqrt(distance);
376 :
377 : // calculate weight
378 0 : double weight = exp(-distance);
379 0 : sumWeights += weight;
380 :
381 : // reference object id
382 0 : const ReferenceObjectID roid = element->reference_object_id();
383 0 : DimReferenceElement<dim> refElem = ReferenceElementProvider::get<dim>(roid);
384 :
385 : // calculate the elements center in local coordinates
386 : MathVector<dim> centerLocal;
387 0 : for(size_t i = 0; i < refElem.corners()->size(); i++)
388 : {
389 0 : for(int d = 0; d < dim; d++)
390 0 : centerLocal[d] += refElem.corners()[i][d] / refElem.corners()->size();
391 : }
392 :
393 : // get subset
394 0 : int si = u->domain()->subset_handler()->get_subset_index(element);
395 :
396 : // storage for the result
397 : MathVector<dim> value;
398 :
399 : // get local solution if needed
400 0 : if(m_userData->requires_grid_fct())
401 : {
402 : // create storage
403 : LocalIndices ind;
404 : LocalVector localU;
405 :
406 : // get global indices
407 : u->indices(element, ind);
408 :
409 : // adapt local algebra
410 0 : localU.resize(ind);
411 :
412 : // read local values of u
413 0 : GetLocalVector(localU, *u);
414 :
415 : try
416 : {
417 0 : (*m_userData)(&value, ¢erGlobal, time, si, element,
418 : &vCornerCoords[0], ¢erLocal, 1, &localU, NULL);
419 : }
420 0 : UG_CATCH_THROW("VectorValuedUserDataEvaluator: Cannot evaluate data.");
421 : }
422 : else
423 : {
424 : try
425 : {
426 0 : (*m_userData)(&value, ¢erGlobal, time, si, element,
427 : &vCornerCoords[0], ¢erLocal, 1, NULL, NULL);
428 : }
429 0 : UG_CATCH_THROW("VectorValuedUserDataEvaluator: Cannot evaluate data.");
430 : }
431 :
432 : // add the vector with the calculated weight
433 0 : for(int i = 0; i < dim; i++)
434 : {
435 0 : result[i] = value[i] * weight;
436 : }
437 : }
438 :
439 : // scale vector by sumWeights
440 0 : for(int i = 0; i < dim; i++)
441 : {
442 0 : result[i] = result[i] / sumWeights;
443 : }
444 :
445 : return true;
446 : }
447 :
448 0 : bool evaluateOnThisProcess(const std::vector<number>& pos,
449 : std::vector<number>& result,
450 : SmartPtr<TGridFunction> u,
451 : number time)
452 : {
453 0 : result.resize(TDomain::dim);
454 :
455 0 : TElem* elem = NULL;
456 :
457 : MathVector<dim> globalPosition;
458 :
459 0 : for(int i = 0; i < dim; i++)
460 : {
461 0 : globalPosition[i] = pos[i];
462 : }
463 :
464 0 : if(!FindContainingElement(elem, *m_tree, globalPosition))
465 : {
466 : return false;
467 : }
468 :
469 : // get corners of element
470 : std::vector<MathVector<dim> > vCornerCoords;
471 0 : CollectCornerCoordinates(vCornerCoords, *elem, *u->domain());
472 :
473 : // get subset
474 0 : int si = u->domain()->subset_handler()->get_subset_index(elem);
475 :
476 : // reference object id
477 0 : const ReferenceObjectID roid = elem->reference_object_id();
478 :
479 : // get local position of DoF
480 : DimReferenceMapping<dim, dim>& map
481 : = ReferenceMappingProvider::get<dim, dim>(roid, vCornerCoords);
482 : MathVector<dim> locPos;
483 : VecSet(locPos, 0.5);
484 0 : map.global_to_local(locPos, globalPosition);
485 :
486 : // storage for the result
487 : MathVector<dim> value;
488 :
489 : // get local solution if needed
490 0 : if(m_userData->requires_grid_fct())
491 : {
492 : // create storage
493 : LocalIndices ind;
494 : LocalVector localU;
495 :
496 : // get global indices
497 : u->indices(elem, ind);
498 :
499 : // adapt local algebra
500 0 : localU.resize(ind);
501 :
502 : // read local values of u
503 0 : GetLocalVector(localU, *u);
504 :
505 : try
506 : {
507 0 : (*m_userData)(&value, &globalPosition, time, si, elem,
508 : &vCornerCoords[0], &locPos, 1, &localU, NULL);
509 : }
510 0 : UG_CATCH_THROW("VectorValuedUserDataEvaluator: Cannot evaluate data.");
511 : }
512 : else
513 : {
514 : try
515 : {
516 0 : (*m_userData)(&value, &globalPosition, time, si, elem,
517 : &vCornerCoords[0], &locPos, 1, NULL, NULL);
518 : }
519 0 : UG_CATCH_THROW("VectorValuedUserDataEvaluator: Cannot evaluate data.");
520 : }
521 :
522 0 : for(int i = 0; i < dim; i++)
523 : {
524 0 : result[i] = value[i];
525 : }
526 :
527 : return true;
528 0 : }
529 :
530 0 : void initialize(SmartPtr<TGridFunction> u)
531 : {
532 0 : m_initialized = true;
533 :
534 0 : m_tree = make_sp(new tree_t(*u->domain()->grid(), u->domain()->position_attachment()));
535 0 : m_tree->create_tree(u->template begin<TElem>(), u->template end<TElem>());
536 0 : }
537 :
538 : bool m_initialized = false;
539 : SmartPtr<tree_t> m_tree;
540 : SmartPtr<UserData<MathVector<dim>, dim> > m_userData;
541 :
542 : };
543 :
544 : template <typename TDomain, typename TAlgebra>
545 : class PointEvaluatorBase : public ITimeIntegratorObserver<TDomain, TAlgebra>
546 : {
547 : typedef std::vector<number> TPoint;
548 : static const int dim = TDomain::dim;
549 : typedef GridFunction<TDomain, TAlgebra> TGridFunction;
550 :
551 : public:
552 0 : void add_evaluation_point(TPoint point)
553 : {
554 0 : UG_COND_THROW(m_initialized, "PointEvaluator: Can't add point, Evaluator in use.");
555 :
556 0 : m_evaluationPoints.push_back(point);
557 0 : }
558 :
559 0 : void set_filename(std::string filename)
560 : {
561 0 : UG_COND_THROW(m_initialized, "PointEvaluator: Can't change filename, Evaluator in use.");
562 :
563 0 : m_filename = filename;
564 0 : }
565 :
566 0 : void set_separator(std::string separator)
567 : {
568 0 : UG_COND_THROW(m_initialized, "PointEvaluator: Can't change separator, Evaluator in use.");
569 :
570 0 : m_separator = separator;
571 0 : }
572 :
573 0 : virtual bool step_process(SmartPtr<TGridFunction> uNew, int step, number time, number dt) override
574 : {
575 0 : if(!m_initialized)
576 0 : initialize();
577 :
578 0 : std::fstream output_file;
579 :
580 : #ifdef UG_PARALLEL
581 : if(pcl::ProcRank() != 0)
582 : {
583 : // passing an unopenend file object. nothing will be written!
584 : write_evaluations(output_file, uNew, step, time, dt);
585 : return true;
586 : }
587 : #endif
588 :
589 0 : output_file.open(m_filename.c_str(), std::fstream::app | std::fstream::out);
590 :
591 0 : UG_COND_THROW(output_file.fail(), "PointEvaluator: Can't open output file.");
592 :
593 : output_file.precision(15);
594 0 : write_evaluations(output_file, uNew, step, time, dt);
595 :
596 0 : output_file.close();
597 :
598 0 : return true;
599 0 : }
600 :
601 : protected:
602 :
603 0 : void initialize()
604 : {
605 0 : m_initialized = true;
606 :
607 : #ifdef UG_PARALLEL
608 : if(pcl::ProcRank() != 0)
609 : {
610 : return;
611 : }
612 : #endif
613 :
614 0 : std::fstream output_file;
615 0 : output_file.open(m_filename.c_str(), std::fstream::trunc | std::fstream::out);
616 :
617 0 : UG_COND_THROW(output_file.fail(), "PointEvaluator: Can't open output file.");
618 :
619 0 : write_header(output_file);
620 :
621 0 : output_file.close();
622 0 : }
623 :
624 0 : virtual void write_header(std::ostream& output) {}
625 0 : virtual void write_evaluations(std::ostream& output, SmartPtr<TGridFunction> uNew, int step, number time, number dt) {}
626 :
627 0 : void write_point(std::ostream& output, const TPoint& point)
628 : {
629 0 : output << "{";
630 0 : for(int i = 0; i < dim; i++)
631 : {
632 0 : if(i > 0)
633 0 : output << ",";
634 0 : output << point[i];
635 : }
636 0 : output << "}";
637 0 : }
638 :
639 : bool m_initialized = false;
640 :
641 : std::vector<TPoint> m_evaluationPoints;
642 : std::string m_filename;
643 : std::string m_separator{"\t"};
644 : };
645 :
646 : template <typename TDomain, typename TAlgebra>
647 : class VectorValuedUserDataPointEvaluator : public PointEvaluatorBase<TDomain, TAlgebra>
648 : {
649 : static const int dim = TDomain::dim;
650 : typedef GridFunction<TDomain, TAlgebra> TGridFunction;
651 :
652 : public:
653 0 : VectorValuedUserDataPointEvaluator(SmartPtr<UserData<MathVector<dim>, dim> > userData)
654 0 : : m_evaluator(userData)
655 0 : {}
656 :
657 : void set_interpolate_over_neighbouring_elements(bool interpolate)
658 : {
659 : m_interpolateOverNeighbouringTriangles = interpolate;
660 : }
661 :
662 : protected:
663 :
664 0 : virtual void write_evaluations(std::ostream& output, SmartPtr<TGridFunction> uNew, int step, number time, number dt) override
665 : {
666 : UG_LOG(" * Write Vector-valued Position Data to '" << this->m_filename << "' ... \n");
667 : output << time << this->m_separator;
668 : std::vector<number> result;
669 0 : for (auto point : this->m_evaluationPoints)
670 : {
671 : result.clear();
672 0 : if(m_evaluator.evaluate(point, result, uNew, time, m_interpolateOverNeighbouringTriangles))
673 : {
674 0 : for(int d = 0; d < dim; d++)
675 : {
676 0 : output << result[d] << this->m_separator;
677 : }
678 : }
679 : else
680 : {
681 0 : for(int d = 0; d < dim; d++)
682 : {
683 : output << "NaN" << this->m_separator;
684 : }
685 : }
686 :
687 : }
688 0 : output << "\n";
689 0 : }
690 :
691 0 : virtual void write_header(std::ostream& output) override
692 : {
693 0 : output << "# Position Evaluating file - vector valued\n";
694 0 : char axis[3] = { 'x', 'y', 'z'};
695 :
696 : output << "time" << this->m_separator;
697 :
698 0 : for (auto point : this->m_evaluationPoints)
699 : {
700 0 : for(int d = 0; d < dim; d++)
701 : {
702 0 : this->write_point(output, point);
703 0 : output << "-" << axis[d] << this->m_separator;
704 : }
705 : }
706 :
707 0 : output << "\n";
708 0 : }
709 :
710 : private:
711 : VectorValuedUserDataEvaluator<TDomain, TAlgebra> m_evaluator;
712 : bool m_interpolateOverNeighbouringTriangles = false;
713 :
714 : };
715 :
716 : template <typename TDomain, typename TAlgebra>
717 : class NumberValuedUserDataPointEvaluator : public PointEvaluatorBase<TDomain, TAlgebra>
718 : {
719 : static const int dim = TDomain::dim;
720 : typedef GridFunction<TDomain, TAlgebra> TGridFunction;
721 :
722 : public:
723 0 : NumberValuedUserDataPointEvaluator(SmartPtr<UserData<number, dim> > userData) : m_evaluator(userData)
724 0 : {}
725 :
726 : protected:
727 :
728 0 : virtual void write_evaluations(std::ostream& output, SmartPtr<TGridFunction> uNew, int step, number time, number dt) override
729 : {
730 : UG_LOG(" * Write Number-valued Position Data to '" << this->m_filename << "' ... \n");
731 : output << time << this->m_separator;
732 : number result;
733 0 : for (auto point : this->m_evaluationPoints)
734 : {
735 0 : if(m_evaluator.evaluate(point, result, uNew, time))
736 : {
737 0 : output << result << this->m_separator;
738 : }
739 : else
740 : {
741 : output << "NaN" << this->m_separator;
742 : }
743 :
744 : }
745 0 : output << "\n";
746 0 : }
747 :
748 0 : virtual void write_header(std::ostream& output) override
749 : {
750 0 : output << "# Position Evaluating file - number valued\n";
751 : output << "time" << this->m_separator;
752 :
753 0 : for (auto point : this->m_evaluationPoints)
754 : {
755 0 : this->write_point(output, point);
756 : }
757 :
758 0 : output << "\n";
759 0 : }
760 :
761 : private:
762 : NumberValuedUserDataEvaluator<TDomain, TAlgebra> m_evaluator;
763 : };
764 :
765 : //! This is a factory for creating a 'PointEvaluatorBase' object from user data.
766 : /*! The class carries TDomain/TAlgebra info required when creating objects.*/
767 : template <typename TDomain, typename TAlgebra>
768 : struct PointEvaluatorFactory
769 : {
770 :
771 : PointEvaluatorFactory(){}
772 :
773 : typedef PointEvaluatorBase<TDomain,TAlgebra> return_type;
774 : typedef UserData<MathVector<TDomain::dim>, TDomain::dim> input_vector_data;
775 : typedef UserData<number, TDomain::dim> input_number_data;
776 :
777 0 : SmartPtr<return_type> create(SmartPtr<input_vector_data> userData) const
778 0 : { return make_sp(new VectorValuedUserDataPointEvaluator<TDomain, TAlgebra>(userData)); }
779 :
780 0 : SmartPtr<return_type> create(SmartPtr<input_number_data> userData) const
781 0 : { return make_sp(new NumberValuedUserDataPointEvaluator<TDomain, TAlgebra>(userData)); }
782 :
783 : };
784 :
785 : template <typename TDomain>
786 0 : bool CloseVertexExists(const MathVector<TDomain::dim>& globPos,
787 : TDomain* dom,
788 : const char* subsets,
789 : SmartPtr<typename TDomain::subset_handler_type> sh,
790 : number maxDist)
791 : {
792 : // domain type
793 : typedef TDomain domain_type;
794 : typedef typename domain_type::grid_type grid_type;
795 : typedef typename domain_type::subset_handler_type subset_handler_type;
796 : // get position accessor
797 0 : grid_type* grid = dom->grid().get();
798 :
799 : const typename domain_type::position_accessor_type& aaPos
800 : = dom->position_accessor();
801 :
802 : typename subset_handler_type::template traits<Vertex>::const_iterator iterEnd, iter;
803 : number minDistanceSq = numeric_limits<number>::max();
804 :
805 : #ifdef UG_PARALLEL
806 : DistributedGridManager* dgm = grid->distributed_grid_manager();
807 : #endif
808 :
809 0 : SubsetGroup ssGrp(sh);
810 0 : if(subsets != NULL)
811 0 : ssGrp.add(TokenizeString(subsets));
812 : else
813 0 : ssGrp.add_all();
814 :
815 0 : for(size_t i = 0; i < ssGrp.size(); ++i)
816 : {
817 : // get subset index
818 : const int si = ssGrp[i];
819 : // iterate over all elements
820 0 : for(size_t lvl = 0; lvl < sh->num_levels(); ++lvl){
821 0 : iterEnd = sh->template end<Vertex>(si, lvl);
822 0 : iter = sh->template begin<Vertex>(si, lvl);
823 0 : for(; iter != iterEnd; ++iter)
824 : {
825 : // get element
826 : //todo: replace most of the following checks by a spGridFct->contains(...)
827 :
828 : Vertex* vrt = *iter;
829 0 : if(grid->has_children(vrt)) continue;
830 :
831 : #ifdef UG_PARALLEL
832 : if(dgm->is_ghost(vrt)) continue;
833 : if(dgm->contains_status(vrt, INT_H_SLAVE)) continue;
834 : #endif
835 : // global position
836 0 : number buffer = VecDistanceSq(globPos, aaPos[vrt]);
837 0 : if(buffer < minDistanceSq)
838 : {
839 : minDistanceSq = buffer;
840 : }
841 : }
842 : }
843 : }
844 0 : return minDistanceSq < sq(maxDist);
845 : }
846 :
847 : /**
848 : * \defgroup interpolate_bridge Interpolation Bridge
849 : * \ingroup disc_bridge
850 : * \{
851 : */
852 :
853 : template <typename TGridFunction>
854 0 : number EvaluateAtVertex(const MathVector<TGridFunction::dim>& globPos,
855 : SmartPtr<TGridFunction> spGridFct,
856 : size_t fct,
857 : const SubsetGroup& ssGrp,
858 : typename TGridFunction::domain_type::subset_handler_type* sh,
859 : bool minimizeOverAllProcs = false)
860 : {
861 :
862 : // domain type
863 : typedef typename TGridFunction::domain_type domain_type;
864 : typedef typename domain_type::grid_type grid_type;
865 : typedef typename domain_type::subset_handler_type subset_handler_type;
866 : // get position accessor
867 0 : domain_type* dom = spGridFct->domain().get();
868 0 : grid_type* grid = dom->grid().get();
869 :
870 0 : subset_handler_type* domSH = dom->subset_handler().get();
871 :
872 : const typename domain_type::position_accessor_type& aaPos
873 : = dom->position_accessor();
874 :
875 : std::vector<DoFIndex> ind;
876 : Vertex* chosen = NULL;
877 :
878 : typename subset_handler_type::template traits<Vertex>::const_iterator iterEnd, iter;
879 : number minDistanceSq = std::numeric_limits<double>::max();
880 :
881 : #ifdef UG_PARALLEL
882 : DistributedGridManager* dgm = grid->distributed_grid_manager();
883 : #endif
884 :
885 0 : for(size_t i = 0; i < ssGrp.size(); ++i)
886 : {
887 : // get subset index
888 : const int si = ssGrp[i];
889 :
890 : // iterate over all elements
891 0 : for(size_t lvl = 0; lvl < sh->num_levels(); ++lvl){
892 0 : iterEnd = sh->template end<Vertex>(si, lvl);
893 0 : iter = sh->template begin<Vertex>(si, lvl);
894 0 : for(; iter != iterEnd; ++iter)
895 : {
896 : // get element
897 : //todo: replace most of the following checks by a spGridFct->contains(...)
898 : Vertex* vrt = *iter;
899 0 : if(grid->has_children(vrt)) continue;
900 :
901 : #ifdef UG_PARALLEL
902 : if(dgm->is_ghost(vrt)) continue;
903 : if(dgm->contains_status(vrt, INT_H_SLAVE)) continue;
904 : #endif
905 :
906 0 : int domSI = domSH->get_subset_index(vrt);
907 :
908 : // skip if function is not defined in subset
909 0 : if(!spGridFct->is_def_in_subset(fct, domSI)) continue;
910 :
911 : // global position
912 0 : number buffer = VecDistanceSq(globPos, aaPos[vrt]);
913 0 : if(buffer < minDistanceSq)
914 : {
915 : minDistanceSq = buffer;
916 : chosen = *iter;
917 : }
918 : }
919 : }
920 : }
921 :
922 : // get corresponding value (if vertex found, otherwise take 0)
923 : number value = 0.0;
924 0 : if (chosen)
925 : {
926 : spGridFct->inner_dof_indices(chosen, fct, ind);
927 0 : value = DoFRef(*spGridFct, ind[0]);
928 : }
929 :
930 : // in parallel environment, find global minimal distance and corresponding value
931 : #ifdef UG_PARALLEL
932 : if (minimizeOverAllProcs)
933 : {
934 : pcl::MinimalKeyValuePairAcrossAllProcs<number, number, std::less<number> >(minDistanceSq, value);
935 :
936 : // check that a vertex has been found
937 : UG_COND_THROW(minDistanceSq == std::numeric_limits<double>::max(),
938 : "No vertex of given subsets could be located on any process.");
939 :
940 : return value;
941 : }
942 : #endif
943 :
944 : // check that a vertex has been found
945 0 : UG_COND_THROW(!chosen, "No vertex of given subsets could be located.");
946 :
947 0 : return value;
948 : }
949 :
950 :
951 : template <typename TGridFunction>
952 0 : number EvaluateAtClosestVertex(const MathVector<TGridFunction::dim>& pos,
953 : SmartPtr<TGridFunction> spGridFct, const char* cmp,
954 : const char* subsets,
955 : SmartPtr<typename TGridFunction::domain_type::subset_handler_type> sh)
956 : {
957 : // get function id of name
958 : const size_t fct = spGridFct->fct_id_by_name(cmp);
959 :
960 : // check that function found
961 0 : if(fct > spGridFct->num_fct())
962 0 : UG_THROW("Evaluate: Name of component '"<<cmp<<"' not found.");
963 :
964 :
965 : // create subset group
966 0 : SubsetGroup ssGrp(sh);
967 0 : if(subsets != NULL)
968 : {
969 0 : ssGrp.add(TokenizeString(subsets));
970 : }
971 : else
972 : {
973 : // add all subsets and remove lower dim subsets afterwards
974 0 : ssGrp.add_all();
975 : }
976 :
977 0 : return EvaluateAtVertex<TGridFunction>(pos, spGridFct, fct, ssGrp, sh.get());
978 : }
979 :
980 :
981 : template <typename TGridFunction>
982 0 : number EvaluateAtClosestVertexAllProcs
983 : (
984 : const MathVector<TGridFunction::dim>& pos,
985 : SmartPtr<TGridFunction> spGridFct,
986 : const char* cmp,
987 : const char* subsets,
988 : SmartPtr<typename TGridFunction::domain_type::subset_handler_type> sh
989 : )
990 : {
991 : // get function id of name
992 : const size_t fct = spGridFct->fct_id_by_name(cmp);
993 :
994 : // check that function found
995 0 : if (fct > spGridFct->num_fct())
996 0 : UG_THROW("Evaluate: Name of component '"<<cmp<<"' not found.");
997 :
998 :
999 : // create subset group
1000 0 : SubsetGroup ssGrp(sh);
1001 0 : if (subsets != NULL)
1002 0 : ssGrp.add(TokenizeString(subsets));
1003 : else
1004 : {
1005 : // add all subsets and remove lower dim subsets afterwards
1006 0 : ssGrp.add_all();
1007 : }
1008 :
1009 0 : return EvaluateAtVertex<TGridFunction>(pos, spGridFct, fct, ssGrp, sh.get(), true);
1010 : }
1011 :
1012 :
1013 : /**
1014 : * Class exporting the functionality. All functionality that is to
1015 : * be used in scripts or visualization must be registered here.
1016 : */
1017 : struct Functionality
1018 : {
1019 :
1020 : /**
1021 : * Function called for the registration of Domain and Algebra dependent parts.
1022 : * All Functions and Classes depending on both Domain and Algebra
1023 : * are to be placed here when registering. The method is called for all
1024 : * available Domain and Algebra types, based on the current build options.
1025 : *
1026 : * @param reg registry
1027 : * @param parentGroup group for sorting of functionality
1028 : */
1029 : template <typename TDomain, typename TAlgebra>
1030 9 : static void DomainAlgebra(Registry& reg, string grp)
1031 : {
1032 : string suffix = GetDomainAlgebraSuffix<TDomain,TAlgebra>();
1033 : string tag = GetDomainAlgebraTag<TDomain,TAlgebra>();
1034 :
1035 : // typedef
1036 : typedef ug::GridFunction<TDomain, TAlgebra> TFct;
1037 :
1038 : {
1039 27 : reg.add_function("EvaluateAtClosestVertex",
1040 : &EvaluateAtClosestVertex<TFct>,
1041 : grp, "Evaluate_at_closest_vertex", "Position#GridFunction#Component#Subsets#SubsetHandler");
1042 27 : reg.add_function("EvaluateAtClosestVertexAllProcs",
1043 : &EvaluateAtClosestVertexAllProcs<TFct>,
1044 : grp, "Evaluate_at_closest_vertex", "Position#GridFunction#Component#Subsets#SubsetHandler");
1045 :
1046 : }
1047 :
1048 :
1049 : {
1050 : typedef PointEvaluatorBase<TDomain, TAlgebra> T;
1051 : typedef ITimeIntegratorObserver<TDomain, TAlgebra> TBase;
1052 :
1053 9 : string name = string("PointEvaluatorBase").append(suffix);
1054 :
1055 27 : reg.add_class_<T, TBase>(name, grp)
1056 27 : .add_method("add_evaluation_point", &T::add_evaluation_point, "point", "")
1057 36 : .add_method("set_filename", &T::set_filename, "filename", "")
1058 27 : .add_method("set_separator", &T::set_separator, "separator", "")
1059 9 : .set_construct_as_smart_pointer(true);
1060 27 : reg.add_class_to_group(name, "PointEvaluatorBase", tag);
1061 : }
1062 :
1063 : {
1064 : typedef VectorValuedUserDataPointEvaluator<TDomain, TAlgebra> T;
1065 : typedef PointEvaluatorBase<TDomain, TAlgebra> TBase;
1066 :
1067 9 : string name = string("VectorValuedUserDataPointEvaluator").append(suffix);
1068 :
1069 27 : reg.add_class_<T, TBase>(name, grp)
1070 18 : .template add_constructor<void (*)(SmartPtr<UserData<MathVector<TDomain::dim>, TDomain::dim> >) >("")
1071 9 : .set_construct_as_smart_pointer(true);
1072 27 : reg.add_class_to_group(name, "VectorValuedUserDataPointEvaluator", tag);
1073 : }
1074 :
1075 : {
1076 : typedef VectorValuedUserDataEvaluator<TDomain, TAlgebra> T;
1077 9 : string name = string("VectorValuedUserDataEvaluator").append(suffix);
1078 :
1079 27 : reg.add_class_<T>(name, grp)
1080 18 : .template add_constructor<void (*)(SmartPtr<UserData<MathVector<TDomain::dim>, TDomain::dim> >) >("")
1081 18 : .add_method("evaluate", &T::evaluateLua, "point#result#solution#time", "")
1082 9 : .set_construct_as_smart_pointer(true);
1083 27 : reg.add_class_to_group(name, "VectorValuedUserDataEvaluator", tag);
1084 : }
1085 :
1086 : {
1087 : typedef NumberValuedUserDataPointEvaluator<TDomain, TAlgebra> T;
1088 : typedef PointEvaluatorBase<TDomain, TAlgebra> TBase;
1089 :
1090 9 : string name = string("NumberValuedUserDataPointEvaluator").append(suffix);
1091 :
1092 27 : reg.add_class_<T, TBase>(name, grp)
1093 18 : .template add_constructor<void (*)(SmartPtr<UserData<number, TDomain::dim> >) >("")
1094 9 : .set_construct_as_smart_pointer(true);
1095 27 : reg.add_class_to_group(name, "NumberValuedUserDataPointEvaluator", tag);
1096 : }
1097 :
1098 : {
1099 : typedef NumberValuedUserDataEvaluator<TDomain, TAlgebra> T;
1100 9 : string name = string("NumberValuedUserDataEvaluator").append(suffix);
1101 :
1102 27 : reg.add_class_<T>(name, grp)
1103 18 : .template add_constructor<void (*)(SmartPtr<UserData<number, TDomain::dim> >) >("")
1104 18 : .add_method("evaluate", &T::evaluateLua, "point#result#solution#time", "")
1105 9 : .set_construct_as_smart_pointer(true);
1106 27 : reg.add_class_to_group(name, "NumberValuedUserDataEvaluator", tag);
1107 : }
1108 : {
1109 : typedef PointEvaluatorFactory<TDomain, TAlgebra> T;
1110 9 : string name = string("PointEvaluatorFactory").append(suffix);
1111 :
1112 27 : reg.add_class_<T>(name, grp)
1113 18 : .template add_constructor<void (*)() >("")
1114 18 : .add_method("create", static_cast<SmartPtr<typename T::return_type> (T::*)(SmartPtr<typename T::input_vector_data>) const>(&T::create), "point#result#solution#time", "")
1115 18 : .add_method("create", static_cast<SmartPtr<typename T::return_type> (T::*)(SmartPtr<typename T::input_number_data>) const>(&T::create), "point#result#solution#time", "")
1116 9 : .set_construct_as_smart_pointer(true);
1117 27 : reg.add_class_to_group(name, "PointEvaluatorFactory", tag);
1118 : }
1119 :
1120 : /*{
1121 : string name = string("GetUserDataPointEvaluator").append(suffix);
1122 : reg.add_function(name,
1123 : static_cast<SmartPtr<PointEvaluatorBase<TDomain,TAlgebra> > (*)(SmartPtr<UserData<MathVector<TDomain::dim>, TDomain::dim> >)>(&GetUserDataPointEvaluator<TDomain,TAlgebra>),
1124 : grp,
1125 : "GetUserDataPointEvaluator",
1126 : "UserDataObject");
1127 :
1128 : reg.add_function(name,
1129 : static_cast<SmartPtr<PointEvaluatorBase<TDomain,TAlgebra> > (*)(SmartPtr<UserData<number, TDomain::dim> >)>(&GetUserDataPointEvaluator<TDomain,TAlgebra>),
1130 : grp,
1131 : "GetUserDataPointEvaluator",
1132 : "UserDataObject");
1133 :
1134 : }*/
1135 9 : }
1136 :
1137 : /**
1138 : * Function called for the registration of Domain dependent parts.
1139 : * All Functions and Classes depending on the Domain
1140 : * are to be placed here when registering. The method is called for all
1141 : * available Domain types, based on the current build options.
1142 : *
1143 : * @param reg registry
1144 : * @param parentGroup group for sorting of functionality
1145 : */
1146 : template <typename TDomain>
1147 3 : static void Domain(Registry& reg, string grp)
1148 : {
1149 : string suffix = GetDomainSuffix<TDomain>();
1150 : string tag = GetDomainTag<TDomain>();
1151 9 : reg.add_function("CloseVertexExists", &CloseVertexExists<TDomain>, grp);
1152 3 : }
1153 :
1154 : /**
1155 : * Function called for the registration of Dimension dependent parts.
1156 : * All Functions and Classes depending on the Dimension
1157 : * are to be placed here when registering. The method is called for all
1158 : * available Dimension types, based on the current build options.
1159 : *
1160 : * @param reg registry
1161 : * @param parentGroup group for sorting of functionality
1162 : */
1163 : template <int dim>
1164 : static void Dimension(Registry& reg, string grp)
1165 : {
1166 : string suffix = GetDimensionSuffix<dim>();
1167 : string tag = GetDimensionTag<dim>();
1168 :
1169 : }
1170 :
1171 : /**
1172 : * Function called for the registration of Algebra dependent parts.
1173 : * All Functions and Classes depending on Algebra
1174 : * are to be placed here when registering. The method is called for all
1175 : * available Algebra types, based on the current build options.
1176 : *
1177 : * @param reg registry
1178 : * @param parentGroup group for sorting of functionality
1179 : */
1180 : template <typename TAlgebra>
1181 : static void Algebra(Registry& reg, string grp)
1182 : {
1183 : string suffix = GetAlgebraSuffix<TAlgebra>();
1184 : string tag = GetAlgebraTag<TAlgebra>();
1185 :
1186 : }
1187 :
1188 : /**
1189 : * Function called for the registration of Domain and Algebra independent parts.
1190 : * All Functions and Classes not depending on Domain and Algebra
1191 : * are to be placed here when registering.
1192 : *
1193 : * @param reg registry
1194 : * @param parentGroup group for sorting of functionality
1195 : */
1196 : static void Common(Registry& reg, string grp)
1197 : {
1198 : }
1199 :
1200 : }; // end Functionality
1201 :
1202 : // end group
1203 : /// \}
1204 :
1205 : }// namespace Evaluate
1206 :
1207 : ///
1208 1 : void RegisterBridge_Evaluate(Registry& reg, string grp)
1209 : {
1210 1 : grp.append("/Evaluate");
1211 : typedef Evaluate::Functionality Functionality;
1212 :
1213 : try{
1214 : // RegisterCommon<Functionality>(reg,grp);
1215 : // RegisterDimensionDependent<Functionality>(reg,grp);
1216 2 : RegisterDomainDependent<Functionality>(reg,grp);
1217 : // RegisterAlgebraDependent<Functionality>(reg,grp);
1218 1 : RegisterDomainAlgebraDependent<Functionality>(reg,grp);
1219 : }
1220 0 : UG_REGISTRY_CATCH_THROW(grp);
1221 1 : }
1222 :
1223 : }// end of namespace bridge
1224 : }// end of namespace ug
|
