Flow123d  release_3.0.0-688-g6b683cf
observe.cc
Go to the documentation of this file.
1 /*
2  * observe.cc
3  *
4  * Created on: Jun 28, 2016
5  * Author: jb
6  */
7 
8 #include <string>
9 #include <numeric>
10 #include <cmath>
11 #include <algorithm>
12 #include <unordered_set>
13 #include <queue>
14 
15 #include "system/global_defs.h"
16 #include "input/accessors.hh"
17 #include "input/input_type.hh"
19 
20 #include "mesh/side_impl.hh"
21 #include "mesh/mesh.h"
22 #include "mesh/bih_tree.hh"
23 #include "mesh/region.hh"
24 #include "mesh/accessors.hh"
25 #include "io/observe.hh"
26 #include "io/element_data_cache.hh"
27 #include "fem/mapping_p1.hh"
28 #include "tools/unit_si.hh"
29 
30 
31 namespace IT = Input::Type;
32 
33 
34 /**
35  * Helper class allows to work with ObservePoint above elements with different dimensions
36  *
37  * Allows:
38  * - calculate projection of points by dimension
39  * - snap to subelement
40  */
41 template<unsigned int dim>
43 public:
44  /// Constructor
46 
47  ObservePointData projection(arma::vec3 input_point, unsigned int i_elm, ElementAccessor<3> elm) {
48  arma::mat::fixed<3, dim+1> elm_map = mapping_.element_map(elm);
49  arma::vec::fixed<dim+1> projection = mapping_.project_real_to_unit(input_point, elm_map);
50  projection = mapping_.clip_to_element(projection);
51 
52  ObservePointData data;
53  data.element_idx_ = i_elm;
54  data.local_coords_ = projection.rows(1, elm.dim());
55  data.global_coords_ = elm_map*projection;
56  data.distance_ = arma::norm(data.global_coords_ - input_point, 2);
57 
58  return data;
59  }
60 
61  /**
62  * Snap local coords to the subelement. Called by the ObservePoint::snap method.
63  */
64  void snap_to_subelement(ObservePointData & observe_data, ElementAccessor<3> elm, unsigned int snap_dim)
65  {
66  if (snap_dim <= dim) {
67  double min_dist = 2.0; // on the ref element the max distance should be about 1.0, smaler then 2.0
68  arma::vec min_center;
69  for(auto &center : RefElement<dim>::centers_of_subelements(snap_dim))
70  {
71  double dist = arma::norm(center - observe_data.local_coords_, 2);
72  if ( dist < min_dist) {
73  min_dist = dist;
74  min_center = center;
75  }
76  }
77  observe_data.local_coords_ = min_center;
78  }
79 
80  arma::mat::fixed<3, dim+1> elm_map = mapping_.element_map(elm);
81  observe_data.global_coords_ = elm_map * RefElement<dim>::local_to_bary(observe_data.local_coords_);
82  }
83 
84 private:
85  /// Mapping object.
87 };
88 
89 template class ProjectionHandler<1>;
90 template class ProjectionHandler<2>;
91 template class ProjectionHandler<3>;
92 
93 
94 /**
95  * Helper struct, used as comparator of priority queue in ObservePoint::find_observe_point.
96  */
98 {
99  bool operator()(const ObservePointData& lhs, const ObservePointData& rhs) const
100  {
101  return lhs.distance_ > rhs.distance_;
102  }
103 };
104 
105 
106 /*******************************************************************
107  * implementation of ObservePoint
108  */
109 
111  return IT::Record("ObservePoint", "Specification of the observation point.\n"
112  "The actual observation element and the observation point on it is determined as follows:\n\n"
113  "1. Find an initial element containing the initial point. If no such element exists, we report an error.\n"
114  "2. Use BFS (Breadth-first search) starting from the inital element to find the 'observe element'. The observe element is the closest element.\n"
115  "3. Find the closest projection of the inital point on the observe element and snap this projection according to the ``snap_dim``.\n")
116  .allow_auto_conversion("point")
117  .declare_key("name", IT::String(),
119  "Default name have the form 'obs_<id>', where 'id' "
120  "is the rank of the point on the input."),
121  "Optional point name, which has to be unique.\n"
122  "Any string that is a valid YAML key in record without any quoting can be used, however, "
123  "using just alpha-numerical characters, and underscore instead of the space, is recommended."
124  )
126  "Initial point for the observe point search.")
127  .declare_key("snap_dim", IT::Integer(0, 4), IT::Default("4"),
128  "The dimension of the sub-element to which center we snap. For value 4 no snapping is done. "
129  "For values 0 up to 3 the element containing the initial point is found and then the observe"
130  "point is snapped to the nearest center of the sub-element of the given dimension. "
131  "E.g. for dimension 2 we snap to the nearest center of the face of the initial element."
132  )
133  .declare_key("snap_region", IT::String(), IT::Default("\"ALL\""),
134  "The region of the initial element for snapping. Without snapping we make a projection to the initial element.")
135  .declare_key("search_radius", IT::Double(0.0),
136  IT::Default::read_time("Maximal distance of observe point from Mesh relative to its size (bounding box). "),
137  "Global value is define in Mesh by the key global_snap_radius.")
138  .close();
139 }
140 
142 {}
143 
144 
145 ObservePoint::ObservePoint(Input::Record in_rec, Mesh &mesh, unsigned int point_idx)
146 {
147  in_rec_ = in_rec;
148 
149  string default_label = string("obs_") + std::to_string(point_idx);
150  name_ = in_rec.val<string>("name", default_label );
151 
152  vector<double> tmp_coords;
153  in_rec.val<Input::Array>("point").copy_to(tmp_coords);
154  input_point_= arma::vec(tmp_coords);
155 
156  snap_dim_ = in_rec.val<unsigned int>("snap_dim");
157 
158  snap_region_name_ = in_rec.val<string>("snap_region");
159 
160  const BoundingBox &main_box = mesh.get_bih_tree().tree_box();
161  double max_mesh_size = arma::max(main_box.max() - main_box.min());
162  max_search_radius_ = in_rec_.val<double>("search_radius", mesh.global_snap_radius()) * max_mesh_size;
163 }
164 
165 
166 
168  return observe_data_.distance_ < numeric_limits<double>::infinity();
169 }
170 
171 
172 
174 {
175  ElementAccessor<3> elm = mesh.element_accessor(observe_data_.element_idx_);
176  switch (elm.dim()) {
177  case 1:
178  {
180  ph.snap_to_subelement(observe_data_, elm, snap_dim_);
181  break;
182  }
183  case 2:
184  {
186  ph.snap_to_subelement(observe_data_, elm, snap_dim_);
187  break;
188  }
189  case 3:
190  {
192  ph.snap_to_subelement(observe_data_, elm, snap_dim_);
193  break;
194  }
195  default: ASSERT(false).error("Clipping supported only for dim=1,2,3.");
196  }
197 }
198 
199 
200 
202  RegionSet region_set = mesh.region_db().get_region_set(snap_region_name_);
203  if (region_set.size() == 0)
204  THROW( RegionDB::ExcUnknownSet() << RegionDB::EI_Label(snap_region_name_) << in_rec_.ei_address() );
205 
206 
207  const BIHTree &bih_tree=mesh.get_bih_tree();
208  vector<unsigned int> candidate_list;
209  std::unordered_set<unsigned int> closed_elements(1023);
210  std::priority_queue< ObservePointData, std::vector<ObservePointData>, CompareByDist > candidate_queue;
211 
212  // search for the initial element
213  auto projected_point = bih_tree.tree_box().project_point(input_point_);
214  bih_tree.find_point(projected_point, candidate_list, true);
215 
216  for (unsigned int i_candidate=0; i_candidate<candidate_list.size(); ++i_candidate) {
217  unsigned int i_elm=candidate_list[i_candidate];
218  ElementAccessor<3> elm = mesh.element_accessor(i_elm);
219 
220  // project point, add candidate to queue
221  auto observe_data = point_projection( i_elm, elm );
222  if (observe_data.distance_ <= max_search_radius_)
223  candidate_queue.push(observe_data);
224  closed_elements.insert(i_elm);
225  }
226 
227  while (!candidate_queue.empty())
228  {
229  auto candidate_data = candidate_queue.top();
230  candidate_queue.pop();
231 
232  unsigned int i_elm=candidate_data.element_idx_;
233  ElementAccessor<3> elm = mesh.element_accessor(i_elm);
234 
235  // test if candidate is in region and update projection
236  if (elm.region().is_in_region_set(region_set)) {
237  ASSERT_LE(candidate_data.distance_, observe_data_.distance_).error();
238 
239  observe_data_.distance_ = candidate_data.distance_;
240  observe_data_.element_idx_ = candidate_data.element_idx_;
241  observe_data_.local_coords_ = candidate_data.local_coords_;
242  observe_data_.global_coords_ = candidate_data.global_coords_;
243  break;
244  }
245 
246  // add candidates to queue
247  for (unsigned int n=0; n < elm->n_nodes(); n++)
248  for(unsigned int i_node_ele : mesh.node_elements()[elm.node_accessor(n).idx()]) {
249  if (closed_elements.find(i_node_ele) == closed_elements.end()) {
250  ElementAccessor<3> neighbor_elm = mesh.element_accessor(i_node_ele);
251  auto observe_data = point_projection( i_node_ele, neighbor_elm );
252  if (observe_data.distance_ <= max_search_radius_)
253  candidate_queue.push(observe_data);
254  closed_elements.insert(i_node_ele);
255  }
256  }
257  }
258 
259  if (! have_observe_element()) {
260  THROW(ExcNoObserveElement() << EI_RegionName(snap_region_name_) );
261  }
262  snap( mesh );
263  ElementAccessor<3> elm = mesh.element_accessor(observe_data_.element_idx_);
264  double dist = arma::norm(elm.centre() - input_point_, 2);
265  double elm_norm = arma::norm(elm.bounding_box().max() - elm.bounding_box().min(), 2);
266  if (dist > 2*elm_norm)
267  WarningOut().fmt("Observe point ({}) is too distant from the mesh.\n", name_);
268 }
269 
270 
271 
272 
273 void ObservePoint::output(ostream &out, unsigned int indent_spaces, unsigned int precision)
274 {
275  out << setw(indent_spaces) << "" << "- name: " << name_ << endl;
276  out << setw(indent_spaces) << "" << " init_point: " << field_value_to_yaml(input_point_, precision) << endl;
277  out << setw(indent_spaces) << "" << " snap_dim: " << snap_dim_ << endl;
278  out << setw(indent_spaces) << "" << " snap_region: " << snap_region_name_ << endl;
279  out << setw(indent_spaces) << "" << " observe_point: " << field_value_to_yaml(observe_data_.global_coords_, precision) << endl;
280 }
281 
282 
283 
285  switch (elm.dim()) {
286  case 1:
287  {
289  return ph.projection(input_point_, i_elm, elm);
290  break;
291  }
292  case 2:
293  {
295  return ph.projection(input_point_, i_elm, elm);
296  break;
297  }
298  case 3:
299  {
301  return ph.projection(input_point_, i_elm, elm);
302  break;
303  }
304  default:
305  ASSERT(false).error("Invalid element dimension!");
306  }
307 
308  return ObservePointData(); // Should not happen.
309 }
310 
311 
312 
313 
314 /*******************************************************************
315  * implementation of Observe
316  */
317 
318 Observe::Observe(string observe_name, Mesh &mesh, Input::Array in_array, unsigned int precision, std::string unit_str)
319 : observe_values_time_(numeric_limits<double>::signaling_NaN()),
320  observe_name_(observe_name),
321  precision_(precision)
322 {
323  // in_rec is Output input record.
324 
325  for(auto it = in_array.begin<Input::Record>(); it != in_array.end(); ++it) {
326  ObservePoint point(*it, mesh, points_.size());
327  point.find_observe_point(mesh);
328  points_.push_back( point );
329  observed_element_indices_.push_back(point.observe_data_.element_idx_);
330  }
331  // make indices unique
332  std::sort(observed_element_indices_.begin(), observed_element_indices_.end());
333  auto last = std::unique(observed_element_indices_.begin(), observed_element_indices_.end());
335 
336  time_unit_str_ = unit_str;
338 
339  if (points_.size() == 0) return;
341  if (rank_==0) {
342  FilePath observe_file_path(observe_name_ + "_observe.yaml", FilePath::output_file);
343  try {
344  observe_file_path.open_stream(observe_file_);
345  //observe_file_.setf(std::ios::scientific);
346  observe_file_.precision(this->precision_);
347 
348  } INPUT_CATCH(FilePath::ExcFileOpen, FilePath::EI_Address_String, in_array)
349  output_header();
350  }
351 }
352 
354  observe_file_.close();
355 }
356 
357 
358 template <typename T>
359 ElementDataCache<T> & Observe::prepare_compute_data(std::string field_name, double field_time, unsigned int n_rows,
360  unsigned int n_cols)
361 {
364  else
365  ASSERT(fabs(field_time / time_unit_seconds_ - observe_values_time_) < 2*numeric_limits<double>::epsilon())
366  (field_time)(observe_values_time_);
367 
368  OutputDataFieldMap::iterator it=observe_field_values_.find(field_name);
369  if (it == observe_field_values_.end()) {
370  observe_field_values_[field_name]
371  = std::make_shared< ElementDataCache<T> >(field_name, n_rows, n_cols, points_.size());
372  it=observe_field_values_.find(field_name);
373  }
374  return dynamic_cast<ElementDataCache<T> &>(*(it->second));
375 }
376 
377 // explicit instantiation of template method
378 #define OBSERVE_PREPARE_COMPUTE_DATA(TYPE) \
379 template ElementDataCache<TYPE> & Observe::prepare_compute_data<TYPE>(std::string field_name, double field_time, \
380  unsigned int n_rows, unsigned int n_cols)
381 
383 OBSERVE_PREPARE_COMPUTE_DATA(unsigned int);
385 
386 
388  unsigned int indent = 2;
389  observe_file_ << "# Observation file: " << observe_name_ << endl;
390  observe_file_ << "time_unit: " << time_unit_str_ << endl;
391  observe_file_ << "time_unit_in_seconds: " << time_unit_seconds_ << endl;
392  observe_file_ << "points:" << endl;
393  for(auto &point : points_)
394  point.output(observe_file_, indent, precision_);
395  observe_file_ << "data:" << endl;
396 
397 }
398 
399 void Observe::output_time_frame(double time) {
400  if (points_.size() == 0) return;
401 
402  if ( ! no_fields_warning ) {
403  no_fields_warning=true;
404  // check that observe fields are set
406  // first call and no fields
407  ASSERT(observe_field_values_.size() == 0);
408  WarningOut() << "No observe fields for the observation stream: " << observe_name_ << endl;
409  }
410  }
411 
413  ASSERT(observe_field_values_.size() == 0);
414  return;
415  }
416 
417  if (rank_ == 0) {
418  unsigned int indent = 2;
419  observe_file_ << setw(indent) << "" << "- time: " << observe_values_time_ << endl;
420  for(auto &field_data : observe_field_values_) {
421  observe_file_ << setw(indent) << "" << " " << field_data.second->field_input_name() << ": ";
422  field_data.second->print_all_yaml(observe_file_, precision_);
423  observe_file_ << endl;
424  }
425  }
426 
427  observe_values_time_ = numeric_limits<double>::signaling_NaN();
428 
429 }
ObservePointData point_projection(unsigned int i_elm, ElementAccessor< 3 > elm)
Project point to given element by dimension of this element.
Definition: observe.cc:284
std::vector< ObservePoint > points_
Full information about observe points.
Definition: observe.hh:234
Class MappingP1 implements the affine transformation of the unit cell onto the actual cell...
Iterator< ValueType > begin() const
Bounding box in 3d ambient space.
Definition: bounding_box.hh:54
#define OBSERVE_PREPARE_COMPUTE_DATA(TYPE)
Definition: observe.cc:378
OutputDataFieldMap observe_field_values_
Stored field values.
Definition: observe.hh:239
unsigned int n_nodes() const
Definition: elements.h:129
Observe(string observe_name, Mesh &mesh, Input::Array in_array, unsigned int precision, std::string unit_str)
Definition: observe.cc:318
vector< vector< unsigned int > > const & node_elements()
Definition: mesh.cc:991
std::ofstream observe_file_
Output file stream.
Definition: observe.hh:249
RegionSet get_region_set(const std::string &set_name) const
Definition: region.cc:329
Accessor to input data conforming to declared Array.
Definition: accessors.hh:567
string field_value_to_yaml(const T &mat, unsigned int prec)
NodeAccessor< 3 > node_accessor(unsigned int ni) const
Definition: accessors.hh:149
BaryPoint project_real_to_unit(const RealPoint &point, const ElementMap &map) const
Definition: mapping_p1.cc:275
Class Input::Type::Default specifies default value of keys of a Input::Type::Record.
Definition: type_record.hh:61
ProjectionHandler()
Constructor.
Definition: observe.cc:45
double convert_unit_from(std::string actual_unit) const
Convert and check user-defined unit.
Definition: unit_si.cc:217
static const Input::Type::Record & get_input_type()
Definition: observe.cc:110
#define INPUT_CATCH(ExceptionType, AddressEITag, input_accessor)
Definition: accessors.hh:64
bool have_observe_element()
Definition: observe.cc:167
static Default obligatory()
The factory function to make an empty default value which is obligatory.
Definition: type_record.hh:110
double global_snap_radius() const
Maximal distance of observe point from Mesh relative to its size.
Definition: mesh.cc:933
Definition: mesh.h:80
void output(ostream &out, unsigned int indent_spaces, unsigned int precision)
Definition: observe.cc:273
#define ASSERT_LE(a, b)
Definition of comparative assert macro (Less or Equal)
Definition: asserts.hh:303
unsigned int element_idx_
Final element of the observe point. The index in the mesh.
Definition: observe.hh:38
std::string time_unit_str_
String representation of the time unit.
Definition: observe.hh:252
int rank_
Definition: observe.hh:231
const RegionDB & region_db() const
Definition: mesh.h:147
#define ASSERT(expr)
Allow use shorter versions of macro names if these names is not used with external library...
Definition: asserts.hh:346
std::string to_string(const T &value)
Definition: string.h:29
std::string observe_name_
Definition: observe.hh:246
Class for declaration of the integral input data.
Definition: type_base.hh:490
arma::vec local_coords_
Local (barycentric) coordinates on the element.
Definition: observe.hh:47
Record & close() const
Close the Record for further declarations of keys.
Definition: type_record.cc:303
Class for declaration of inputs sequences.
Definition: type_base.hh:346
virtual ElementAccessor< 3 > element_accessor(unsigned int idx) const
Create and return ElementAccessor to element of given idx.
Definition: mesh.cc:719
IteratorBase end() const
arma::vec::fixed< spacedim > centre() const
Computes the barycenter.
Definition: accessors.hh:285
void open_stream(Stream &stream) const
Definition: file_path.cc:211
Global macros to enhance readability and debugging, general constants.
virtual Record & allow_auto_conversion(const string &from_key)
Allows shorter input of the Record providing only value of the from_key given as the parameter...
Definition: type_record.cc:132
Class for declaration of the input data that are floating point numbers.
Definition: type_base.hh:541
arma::vec3 global_coords_
Global coordinates of the observation point.
Definition: observe.hh:44
unsigned int precision_
Precision of float output.
Definition: observe.hh:256
void find_observe_point(Mesh &mesh)
Definition: observe.cc:201
DummyInt isnan(...)
Definition: format.h:306
bool operator()(const ObservePointData &lhs, const ObservePointData &rhs) const
Definition: observe.cc:99
const BoundingBox & tree_box() const
Definition: bih_tree.cc:229
Accessor to the data with type Type::Record.
Definition: accessors.hh:292
const Ret val(const string &key) const
UnitSI & s(int exp=1)
Definition: unit_si.cc:76
void output_header()
Definition: observe.cc:387
std::vector< unsigned int > observed_element_indices_
Elements of the o_points.
Definition: observe.hh:236
const Point & max() const
double distance_
Definition: observe.hh:51
Class for O(log N) lookup for intersections with a set of bounding boxes.
Definition: bih_tree.hh:38
BoundingBox bounding_box() const
Definition: accessors.hh:157
Record & declare_key(const string &key, std::shared_ptr< TypeBase > type, const Default &default_value, const string &description, TypeBase::attribute_map key_attributes=TypeBase::attribute_map())
Declares a new key of the Record.
Definition: type_record.cc:501
void snap_to_subelement(ObservePointData &observe_data, ElementAccessor< 3 > elm, unsigned int snap_dim)
Definition: observe.cc:64
ElementDataCache< T > & prepare_compute_data(std::string field_name, double field_time, unsigned int n_rows, unsigned int n_cols)
Definition: observe.cc:359
double observe_values_time_
Common evaluation time of the fields for single time frame.
Definition: observe.hh:243
Region region() const
Definition: accessors.hh:95
ObservePointData projection(arma::vec3 input_point, unsigned int i_elm, ElementAccessor< 3 > elm)
Definition: observe.cc:47
const BIHTree & get_bih_tree()
Getter for BIH. Creates and compute BIH at first call.
Definition: mesh.cc:924
#define MPI_Comm_rank
Definition: mpi.h:236
Dedicated class for storing path to input and output files.
Definition: file_path.hh:54
static Default read_time(const std::string &description)
The factory function to make an default value that will be specified at the time when a key will be r...
Definition: type_record.hh:97
void find_point(const Space< 3 >::Point &point, std::vector< unsigned int > &result_list, bool full_list=false) const
Definition: bih_tree.cc:287
double time_unit_seconds_
Time unit in seconds.
Definition: observe.hh:254
bool is_in_region_set(const RegionSet &set) const
Definition: region.cc:47
MappingP1< dim, 3 > mapping_
Mapping object.
Definition: observe.cc:86
Point project_point(const Point &point) const
#define WarningOut()
Macro defining &#39;warning&#39; record of log.
Definition: logger.hh:246
#define MPI_COMM_WORLD
Definition: mpi.h:123
static BaryPoint local_to_bary(const LocalPoint &lp)
Converts from local to barycentric coordinates.
Definition: ref_element.cc:286
const Point & min() const
Record type proxy class.
Definition: type_record.hh:182
unsigned int dim() const
Definition: accessors.hh:87
Class for representation SI units of Fields.
Definition: unit_si.hh:40
BaryPoint clip_to_element(BaryPoint &barycentric)
Definition: mapping_p1.cc:297
ElementMap element_map(ElementAccessor< 3 > elm) const
Definition: mapping_p1.cc:265
void snap(Mesh &mesh)
Definition: observe.cc:173
bool no_fields_warning
Definition: observe.hh:259
Class for declaration of the input data that are in string format.
Definition: type_base.hh:589
#define THROW(whole_exception_expr)
Wrapper for throw. Saves the throwing point.
Definition: exceptions.hh:53
~Observe()
Destructor, must close the file.
Definition: observe.cc:353
void output_time_frame(double time)
Definition: observe.cc:399