Flow123d  master-f44eb46
observe.hh
Go to the documentation of this file.
1  /*
2  * observe.hh
3  *
4  * Created on: Jun 28, 2016
5  * Author: jb
6  */
7 
8 #ifndef SRC_IO_OBSERVE_HH_
9 #define SRC_IO_OBSERVE_HH_
10 
11 
12 #include <iosfwd> // for ofstream, ostream
13 #include <map> // for map, map<>::value_compare
14 #include <memory> // for shared_ptr
15 #include <new> // for operator new[]
16 #include <string> // for string, operator<<
17 #include <vector> // for vector
18 #include <armadillo>
19 #include "input/accessors.hh" // for Array (ptr only), Record
20 #include "input/input_exception.hh" // for DECLARE_INPUT_EXCEPTION
21 #include "system/exceptions.hh" // for operator<<, ExcStream, EI
22 #include "system/armadillo_tools.hh" // for Armadillo vec string
23 #include "system/index_types.hh" // for LongIdx
24 #include "mesh/range_wrapper.hh"
26 #include "la/distribution.hh"
27 
29 class Mesh;
30 class TimeUnitConversion;
31 namespace Input { namespace Type { class Record; } }
32 template <typename T> class ElementDataCache;
33 
34 
35 
36 /// Holds data of one eval point on patch (index of element and local coordinations).
38  /// Default constructor
40 
41  /// Constructor with data mebers initialization
42  PatchPointData(unsigned int elm_idx, arma::vec loc_coords)
43  : elem_idx(elm_idx), local_coords(loc_coords), i_quad(0), i_quad_point(0) {}
44 
45  /// Copy constructor
47  : elem_idx(other.elem_idx), local_coords(other.local_coords),
48  i_quad(other.i_quad), i_quad_point(other.i_quad_point) {}
49 
50  unsigned int elem_idx; ///< Index of element
51  arma::vec local_coords; ///< Local coords of point
52  unsigned int i_reg; ///< Index of region (use during patch creating)
53  unsigned int i_quad; ///< Index of quadrature (use during patch creating), i_quad = dim-1
54  unsigned int i_quad_point; ///< Index of point in quadrature (use during patch creating)
55 };
57 
58 
59 
60 /**
61  * Helper class stores base data of ObservePoint and allows to evaluate
62  * the nearest point to input_point_.
63  */
65 public:
66  /// Constructor
68  : distance_(numeric_limits<double>::infinity()) {};
69 
70  /// Final element of the observe point. The index in the mesh.
71  unsigned int element_idx_;
72 
73  /// Global coordinates of the observation point.
75 
76  /// Local (barycentric) coordinates on the element.
78 
79  /// Distance of found projection from the initial point.
80  /// If we find more candidates we pass in the closest one.
81  double distance_;
82 
83  /// Actual process of the observe point.
84  unsigned int proc_;
85 
86  /// Global index of the observe point.
88 
89  /// Local index on actual process of the observe point.
91 };
92 
93 
94 /**
95  * Class representing single observe point, used internally by the class Observe.
96  * Members: input_pos_, snap_dim_, snap_region_name_ are set in constructor. Should be checked before passed in.
97  * Members: element_idx_, global_coords_, local_coords_ are derived, set in Observe::find_observe_points.
98  */
99 class ObservePoint {
100 public:
101  DECLARE_INPUT_EXCEPTION(ExcNoInitialPoint,
102  << "Failed to find the element containing the initial observe point.\n");
103  TYPEDEF_ERR_INFO(EI_RegionName, std::string);
104  TYPEDEF_ERR_INFO(EI_PointName, std::string);
107  DECLARE_INPUT_EXCEPTION(ExcNoObserveElementCandidates,
108  << "Failed to find any element in the search radius of the observe point " << EI_PointName::qval
109  << " with given coordinates " << field_value_to_yaml(EI_Point::ref(*this)) << ".\n"
110  << "The closest element has index " << EI_ClosestEle::ref(*this).element_idx_ << ", its distance is " << EI_ClosestEle::ref(*this).distance_ << ".\n"
111  << "Solution: check the position of the observe point, possibly increase the maximal snapping distance "
112  << "(keys: observe_points:search_radius, mesh:global_snap_radius)"<< "\n");
113  DECLARE_INPUT_EXCEPTION(ExcNoObserveElement,
114  << "Failed to find any element in the search radius of the observe point" << EI_PointName::qval
115  << " inside the snap region: " << EI_RegionName::qval << ".\n"
116  << "The observe point coordinates are " << field_value_to_yaml(EI_Point::ref(*this)) << ".\n"
117  << "The closest element (outside the snap region) has index " << EI_ClosestEle::ref(*this).element_idx_
118  << ", its distance is " << EI_ClosestEle::ref(*this).distance_ << ".\n"
119  << "Solution: check the position/region of the observe point, possibly increase the maximal snapping distance "
120  << "(keys: observe_points:search_radius, mesh:global_snap_radius)"<< "\n");
121 
122  static const Input::Type::Record & get_input_type();
123 
124  /**
125  * Return index of observation point in the mesh.
126  */
127  inline unsigned int element_idx() const
128  { return observe_data_.element_idx_; }
129 
130  /**
131  * Return Local coordinates of the observation point.
132  */
133  inline arma::vec local_coords() const
134  { return observe_data_.local_coords_; }
135 
136  /**
137  * Return global coordinates of the observation point.
138  */
139  inline arma::vec3 global_coords() const
140  { return observe_data_.global_coords_; }
141 
142 protected:
143  /**
144  * Default constructor just for testing.
145  */
146  ObservePoint();
147 
148  /**
149  * Constructor. Read from input.
150  */
151  ObservePoint(Input::Record in_rec, Mesh &mesh, unsigned int point_idx);
152 
153  /**
154  * Returns true if we have already found any observe element.
155  */
156  bool have_observe_element();
157 
158  /**
159  * Snap to the center of closest subelement with dimension snap_dim_.
160  * This makes final adjustment of global_coords_ and local_coords_.
161  */
162  void snap(Mesh &mesh);
163 
164 
165  /**
166  * Find the observe element and the definitive observe point.
167  *
168  * Algorithm:
169  * 1. find element containing the point (initial element)
170  * 2. check initial element for region match possibly set it as (observe element)
171  * 3. add neighbours into processing_list for the next level
172  * 4. while we have no observe element: pass through the processing list of the current level
173  * 5. if element match the region, project and clip the init point, update observe element.
174  * 6. snapping on the observe element
175  *
176  */
177  void find_observe_point(Mesh &mesh);
178 
179  /**
180  * Output the observe point information into a YAML formated stream, indent by
181  * given number of spaces + "- ".
182  */
183  void output(ostream &out, unsigned int indent_spaces, unsigned int precision);
184 
185  /// Project point to given element by dimension of this element.
186  ObservePointData point_projection(unsigned int i_elm, ElementAccessor<3> elm);
187 
188  /// Index in the input array.
190 
191  /// Observation point name.
192  std::string name_;
193 
194  /**
195  * Snap to the center of the object of given dimension.
196  * Value 4 and greater means no snapping.
197  */
198  unsigned int snap_dim_;
199 
200  /**
201  * Region of the snapping element.
202  */
204 
205  /**
206  * Maximal distance of observe element from input point.
207  */
209 
210  /// Input coordinates of the initial position of the observation point.
212 
213  /// Helper object stored projection data
215 
216  /// Only Observe should use this class directly.
217  friend class Observe;
218 
219 };
220 
221 
223 
224 /**
225  * This class takes care about the observe points in the output stream, storing observe values of the fields and
226  * their output in the YAML format.
227  */
228 class Observe {
229 public:
230 
231  typedef std::shared_ptr<ElementDataCacheBase> OutputDataPtr;
233 
234  /**
235  * Construct the observation object.
236  *
237  * observe_name - base name of the output file, the equation name.
238  * mesh - the mesh used for search for the observe points
239  * in_array - the array of observe points
240  */
241  Observe(string observe_name, Mesh &mesh, Input::Array in_array,
242  unsigned int precision, const std::shared_ptr<TimeUnitConversion>& time_unit_conv);
243 
244  /// Destructor, must close the file.
245  ~Observe();
246 
247 
248  /**
249  * Provides a vector of element indices on which the observation values are evaluated.
250  * This can be used to evaluate derived fields only on these elements in the times not selected to
251  * full output.
252  */
254  { return observed_element_indices_;}
255 
256  /**
257  * Output file header.
258  */
259  void output_header();
260 
261  /**
262  * Sets next output time frame of observe. If the table is full, writes field values to the output
263  * file. Using the YAML format. Argument flush starts writing to output file explicitly.
264  */
265  void output_time_frame(bool flush);
266 
267  /**
268  * Return \p points_ vector
269  */
270  inline const std::vector<ObservePoint> & points() const
271  { return points_; }
272 
273  /**
274  * Return point distribution
275  */
276  inline const Distribution * point_ds() const
277  { return point_ds_; }
278 
279  /// Returns local range of observe points
281 
282  /**
283  * Prepare data for computing observe values.
284  *
285  * Method:
286  * - check that all fields of one time frame are evaluated at the same time
287  * - find and return ElementDataCache of given field_name, create its if doesn't exist.
288  *
289  * @param field_name Quantity name of founding ElementDataCache
290  * @param field_time Actual computing time
291  * @param n_shape Count of rows * cols of data cache (used only if new cache is created)
292  */
293  OutputDataPtr prepare_compute_data(std::string field_name, double field_time, unsigned int n_shape);
294 
295  /**
296  * Return pointer to ElementDataCache of given field.
297  *
298  * @param field_name Quantity name of founding ElementDataCache
299  */
300  inline OutputDataPtr get_output_cache(std::string field_name) {
301  OutputDataFieldMap::iterator it=observe_field_values_.find(field_name);
302  ASSERT(it != observe_field_values_.end())(field_name).error("Element data cache of given field doesn't exist. Did you call prepare_compute_data?\n");
303  return it->second;
304  }
305 
306  /// Getter of patch_point_data.
308  return patch_point_data_;
309  }
310 
311 
312 
313 protected:
314  /// Effectively writes the data into the observe stream.
315  void flush_values();
316 
317  /// Maximal size of observe values times vector
318  static const unsigned int max_observe_value_time;
319 
320  // MPI rank.
321  int rank_;
322 
323  /// Full information about observe points.
325  /// Elements of the o_points.
327 
328  /// Stored field values.
330 
331 
332  /// Common evaluation time of the fields for single time frame.
334 
335  // Name of the observation stream. Base for the output filename.
336  std::string observe_name_;
337 
338  /// Output file stream.
339  std::ofstream observe_file_;
340 
341  /// Precision of float output
342  unsigned int precision_;
343  /// Time unit conversion object.
344  std::shared_ptr<TimeUnitConversion> time_unit_conversion_;
345 
346  // Warn for no observe fields only once.
347  bool no_fields_warning=false;
348 
349  /// Index set assigning to local point index its global index.
351 
352  /// Parallel distribution of observe points.
354 
355  /// Index of actual (last) time in \p observe_values_time_ vector
356  unsigned int observe_time_idx_;
357 
358  /// Holds observe data of eval points on patch
360 
361  friend class ObservePointAccessor;
362 };
363 
364 
365 /**
366  * @brief Point accessor allow iterate over local Observe points.
367  *
368  * Iterator is defined by:
369  * - Observe object
370  * - local index of observe point (iterated value)
371  */
373 public:
374  /// Default invalid accessor.
376  : observe_(NULL)
377  {}
378 
379  /**
380  * Observe point accessor.
381  */
382  ObservePointAccessor(const Observe *observe, unsigned int loc_idx)
383  : observe_(observe), loc_point_idx_(loc_idx)
384  {}
385 
386  /// Return local index to point.
387  inline unsigned int local_idx() const {
388  return loc_point_idx_;
389  }
390 
391  /// Return global index to point.
392  inline unsigned int global_idx() const {
394  }
395 
396  /// Return ElementAccessor to element of loc_ele_idx_.
397  inline const ObservePoint observe_point() const {
398  return observe_->points_[ this->global_idx() ];
399  }
400 
401  /// Return local index in data cache (combination of local point index and index of stored time)
402  inline unsigned int loc_point_time_index() const {
404  }
405 
406  /// Check validity of accessor (see default constructor)
407  inline bool is_valid() const {
408  return observe_ != NULL;
409  }
410 
411  /// Iterates to next local point.
412  inline void inc() {
413  loc_point_idx_++;
414  }
415 
416  /// Comparison of accessors.
417  bool operator==(const ObservePointAccessor& other) {
418  return (loc_point_idx_ == other.loc_point_idx_);
419  }
420 
421 protected:
422  /// Pointer to the Observe owning the point.
423  const Observe * observe_;
424  /// Index into Observe::point_4_loc_ array.
425  unsigned int loc_point_idx_;
426 
427 };
428 
429 
430 
431 
432 #endif /* SRC_IO_OBSERVE_HH_ */
string field_value_to_yaml(const T &mat, unsigned int prec)
#define ASSERT(expr)
Definition: asserts.hh:351
Accessor to input data conforming to declared Array.
Definition: accessors.hh:566
Accessor to the data with type Type::Record.
Definition: accessors.hh:291
Record type proxy class.
Definition: type_record.hh:182
Definition: mesh.h:362
Point accessor allow iterate over local Observe points.
Definition: observe.hh:372
void inc()
Iterates to next local point.
Definition: observe.hh:412
ObservePointAccessor()
Default invalid accessor.
Definition: observe.hh:375
const ObservePoint observe_point() const
Return ElementAccessor to element of loc_ele_idx_.
Definition: observe.hh:397
bool operator==(const ObservePointAccessor &other)
Comparison of accessors.
Definition: observe.hh:417
bool is_valid() const
Check validity of accessor (see default constructor)
Definition: observe.hh:407
unsigned int local_idx() const
Return local index to point.
Definition: observe.hh:387
unsigned int loc_point_idx_
Index into Observe::point_4_loc_ array.
Definition: observe.hh:425
ObservePointAccessor(const Observe *observe, unsigned int loc_idx)
Definition: observe.hh:382
const Observe * observe_
Pointer to the Observe owning the point.
Definition: observe.hh:423
unsigned int global_idx() const
Return global index to point.
Definition: observe.hh:392
unsigned int loc_point_time_index() const
Return local index in data cache (combination of local point index and index of stored time)
Definition: observe.hh:402
double distance_
Definition: observe.hh:81
ObservePointData()
Constructor.
Definition: observe.hh:67
unsigned int element_idx_
Final element of the observe point. The index in the mesh.
Definition: observe.hh:68
unsigned int proc_
Actual process of the observe point.
Definition: observe.hh:84
LongIdx global_idx_
Global index of the observe point.
Definition: observe.hh:87
LongIdx local_idx_
Local index on actual process of the observe point.
Definition: observe.hh:90
arma::vec local_coords_
Local (barycentric) coordinates on the element.
Definition: observe.hh:77
arma::vec3 global_coords_
Global coordinates of the observation point.
Definition: observe.hh:74
arma::vec3 global_coords() const
Definition: observe.hh:139
arma::vec3 input_point_
Input coordinates of the initial position of the observation point.
Definition: observe.hh:211
Input::Record in_rec_
Index in the input array.
Definition: observe.hh:189
DECLARE_INPUT_EXCEPTION(ExcNoObserveElement,<< "Failed to find any element in the search radius of the observe point"<< EI_PointName::qval<< " inside the snap region: "<< EI_RegionName::qval<< ".\n"<< "The observe point coordinates are "<< field_value_to_yaml(EI_Point::ref(*this))<< ".\n"<< "The closest element (outside the snap region) has index "<< EI_ClosestEle::ref(*this).element_idx_<< ", its distance is "<< EI_ClosestEle::ref(*this).distance_<< ".\n"<< "Solution: check the position/region of the observe point, possibly increase the maximal snapping distance "<< "(keys: observe_points:search_radius, mesh:global_snap_radius)"<< "\n")
TYPEDEF_ERR_INFO(EI_ClosestEle, ObservePointData)
unsigned int element_idx() const
Definition: observe.hh:127
void output(ostream &out, unsigned int indent_spaces, unsigned int precision)
Definition: observe.cc:292
DECLARE_INPUT_EXCEPTION(ExcNoObserveElementCandidates,<< "Failed to find any element in the search radius of the observe point "<< EI_PointName::qval<< " with given coordinates "<< field_value_to_yaml(EI_Point::ref(*this))<< ".\n"<< "The closest element has index "<< EI_ClosestEle::ref(*this).element_idx_<< ", its distance is "<< EI_ClosestEle::ref(*this).distance_<< ".\n"<< "Solution: check the position of the observe point, possibly increase the maximal snapping distance "<< "(keys: observe_points:search_radius, mesh:global_snap_radius)"<< "\n")
TYPEDEF_ERR_INFO(EI_RegionName, std::string)
string snap_region_name_
Definition: observe.hh:203
bool have_observe_element()
Definition: observe.cc:164
TYPEDEF_ERR_INFO(EI_PointName, std::string)
unsigned int snap_dim_
Definition: observe.hh:198
arma::vec local_coords() const
Definition: observe.hh:133
std::string name_
Observation point name.
Definition: observe.hh:192
ObservePointData point_projection(unsigned int i_elm, ElementAccessor< 3 > elm)
Project point to given element by dimension of this element.
Definition: observe.cc:303
static const Input::Type::Record & get_input_type()
Definition: observe.cc:107
ObservePointData observe_data_
Helper object stored projection data.
Definition: observe.hh:214
void snap(Mesh &mesh)
Definition: observe.cc:170
DECLARE_INPUT_EXCEPTION(ExcNoInitialPoint,<< "Failed to find the element containing the initial observe point.\n")
double max_search_radius_
Definition: observe.hh:208
TYPEDEF_ERR_INFO(EI_Point, arma::vec3)
void find_observe_point(Mesh &mesh)
Definition: observe.cc:198
const Distribution * point_ds() const
Definition: observe.hh:276
PatchPointVec patch_point_data_
Holds observe data of eval points on patch.
Definition: observe.hh:359
static const unsigned int max_observe_value_time
Maximal size of observe values times vector.
Definition: observe.hh:318
std::vector< ObservePoint > points_
Full information about observe points.
Definition: observe.hh:324
unsigned int observe_time_idx_
Index of actual (last) time in observe_values_time_ vector.
Definition: observe.hh:356
Observe(string observe_name, Mesh &mesh, Input::Array in_array, unsigned int precision, const std::shared_ptr< TimeUnitConversion > &time_unit_conv)
Definition: observe.cc:340
const std::vector< unsigned int > & observed_elements() const
Definition: observe.hh:253
Range< ObservePointAccessor > local_range() const
Returns local range of observe points.
Definition: observe.cc:488
void flush_values()
Effectively writes the data into the observe stream.
Definition: observe.cc:432
bool no_fields_warning
Definition: observe.hh:347
std::vector< unsigned int > observed_element_indices_
Elements of the o_points.
Definition: observe.hh:326
std::ofstream observe_file_
Output file stream.
Definition: observe.hh:339
std::vector< double > observe_values_time_
Common evaluation time of the fields for single time frame.
Definition: observe.hh:333
std::string observe_name_
Definition: observe.hh:336
OutputDataPtr get_output_cache(std::string field_name)
Definition: observe.hh:300
~Observe()
Destructor, must close the file.
Definition: observe.cc:394
int rank_
Definition: observe.hh:321
OutputDataFieldMap observe_field_values_
Stored field values.
Definition: observe.hh:329
std::map< string, OutputDataPtr > OutputDataFieldMap
Definition: observe.hh:232
PatchPointVec & patch_point_data()
Getter of patch_point_data.
Definition: observe.hh:307
std::shared_ptr< TimeUnitConversion > time_unit_conversion_
Time unit conversion object.
Definition: observe.hh:344
unsigned int precision_
Precision of float output.
Definition: observe.hh:342
void output_header()
Definition: observe.cc:420
const std::vector< ObservePoint > & points() const
Definition: observe.hh:270
std::vector< LongIdx > point_4_loc_
Index set assigning to local point index its global index.
Definition: observe.hh:350
Distribution * point_ds_
Parallel distribution of observe points.
Definition: observe.hh:353
std::shared_ptr< ElementDataCacheBase > OutputDataPtr
Definition: observe.hh:231
void output_time_frame(bool flush)
Definition: observe.cc:463
OutputDataPtr prepare_compute_data(std::string field_name, double field_time, unsigned int n_shape)
Definition: observe.cc:401
Range helper class.
Helper class storing unit conversion coefficient and functionality for conversion of units.
Support classes for parallel programing.
Template Iter serves as general template for internal iterators.
int LongIdx
Define type that represents indices of large arrays (elements, nodes, dofs etc.)
Definition: index_types.hh:24
ArmaVec< double, N > vec
Definition: armor.hh:885
Abstract linear system class.
Definition: balance.hh:40
std::vector< PatchPointData > PatchPointVec
Definition: observe.hh:56
Implementation of range helper class.
Holds data of one eval point on patch (index of element and local coordinations).
Definition: observe.hh:37
unsigned int i_reg
Index of region (use during patch creating)
Definition: observe.hh:52
PatchPointData(const PatchPointData &other)
Copy constructor.
Definition: observe.hh:46
unsigned int i_quad
Index of quadrature (use during patch creating), i_quad = dim-1.
Definition: observe.hh:53
PatchPointData(unsigned int elm_idx, arma::vec loc_coords)
Constructor with data mebers initialization.
Definition: observe.hh:42
PatchPointData()
Default constructor.
Definition: observe.hh:39
unsigned int elem_idx
Index of element.
Definition: observe.hh:50
arma::vec local_coords
Local coords of point.
Definition: observe.hh:51
unsigned int i_quad_point
Index of point in quadrature (use during patch creating)
Definition: observe.hh:54