SheafSystem  0.0.0.0
cylindrical_point_locator.cc
Go to the documentation of this file.
1 
2 //
3 // Copyright (c) 2014 Limit Point Systems, Inc.
4 //
5 // Licensed under the Apache License, Version 2.0 (the "License");
6 // you may not use this file except in compliance with the License.
7 // You may obtain a copy of the License at
8 //
9 // http://www.apache.org/licenses/LICENSE-2.0
10 //
11 // Unless required by applicable law or agreed to in writing, software
12 // distributed under the License is distributed on an "AS IS" BASIS,
13 // WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 // See the License for the specific language governing permissions and
15 // limitations under the License.
16 //
17 
20 
21 #include "SheafSystem/cylindrical_point_locator.h"
22 
23 #include "SheafSystem/assert_contract.h"
24 #include "SheafSystem/bilinear_2d.h"
25 #include "SheafSystem/chart_point_3d.h"
26 #include "SheafSystem/d_bin_coordinates.h"
27 #include "SheafSystem/eval_iterator.h"
28 #include "SheafSystem/error_message.h"
29 #include "SheafSystem/linear_2d.h"
30 #include "SheafSystem/preorder_iterator.h"
31 #include "SheafSystem/sec_at1_space.h"
32 #include "SheafSystem/sec_ed.h"
33 #include "SheafSystem/std_bitset.h"
34 
35 using namespace std;
36 using namespace geometry; // Workaround for MS C++ bug.
37 
38 //#undef DIAGNOSTIC_OUTPUT
39 //#define DIAGNOSTIC_OUTPUT
40 
41 
42 
43 // ===========================================================
44 // CYLINDRICAL_POINT_LOCATOR FACET
45 // ===========================================================
46 
47 // PUBLIC MEMBER FUNCTIONS
48 
49 //const int geometry::cylindrical_point_locator::DC;
50 //const int geometry::cylindrical_point_locator::DB;
51 
53 cylindrical_point_locator(sec_ed& xcoords) : point_locator(xcoords),_eval_itr(xcoords.schema(), false)
54  //point_locator::point_locator(xcoords),_eval_itr(xcoords.schema(), false)
55 {
56  // Preconditions:
57 
58  require(precondition_of(point_locator(xcoords)));
59 
60  // Body:
61 
62  _bin_ub.reserve(DB);
63  _bin_ub.set_ct(DB);
64 
65  _bin_size.reserve(DB);
66  _bin_size.set_ct(DB);
67 
70 
71  _box_ct = 0;
72 
73  // The eval_itr owns the evaluator_family which owns the section_evaluators.
74  // The bounding boxes in _boxes keep pointers to these section evaluators.
75  // So the lifetime of the eval_family must be as long as _boxes.
76 
83 
85 
86  // Postconditions:
87 
88 
89  // Exit:
90 
91  return;
92 }
93 
96 {
97  // Preconditions:
98 
99  // Body:
100 
101  // Nothing to do.
102 
103  // Postconditions:
104 
105  // Exit:
106 
107  return;
108 }
109 
112 bin_ub() const
113 {
114  // Preconditions:
115 
116 
117  // Body:
118 
119  const block<size_type>& result = _bin_ub;
120 
121  // Postconditions:
122 
123  ensure_for_all(i, 0, db(), result[i] > 0);
124 
125  // Exit:
126 
127  return result;
128 }
129 
132 bin_size() const
133 {
134 
135  // Preconditions:
136 
137 
138  // Body:
139 
140  const block<sec_vd_value_type>& result = _bin_size;
141 
142  // Postconditions:
143 
144 
145  // Exit:
146 
147  return result;
148 }
149 
152 box_ct() const
153 {
154  size_type result;
155 
156  // Preconditions:
157 
158 
159  // Body:
160 
161  result = _box_ct;
162 
163  // Postconditions:
164 
165 
166  // Exit:
167 
168  return result;
169 }
170 
171 bool
173 is_empty() const
174 {
175  bool result;
176 
177  // Preconditions:
178 
179 
180  // Body:
181 
182  result = (box_ct() == 0);
183 
184  // Postconditions:
185 
186  ensure(result == (box_ct() == 0));
187 
188  // Exit:
189 
190  return result;
191 }
192 
196 {
197  d_bin_coordinates<2, 2>* result;
198 
199  // Preconditions:
200 
201  require(xpt != 0);
202  require(xpt_ub >= db());
203 
204  // Body:
205 
206  result = new d_bin_coordinates<2, 2>();
207  relative_position_pa(xpt, xpt_ub, *result);
208 
209  // Postconditions:
210 
211  ensure(result != 0);
212  ensure(postcondition_of(relative_position_pa(xpt, xpt_ub, result)));
213 
214  // Exit:
215 
216  return result;
217 }
218 
219 void
222  size_type xpt_ub,
223  d_bin_coordinates<2, 2>& xresult) const
224 {
225  // Preconditions:
226 
227  require(xpt != 0);
228  require(xpt_ub >= db());
229  require((-90.0 < xpt[1]) && (xpt[1] < 90.0));
230 
231 
232  // Body:
233 
234  sec_vd_value_type lrel_pos[DB];
235 
236  lrel_pos[0] = (xpt[0] + 180.0)*_one_over_min_bin_size[0]; // - (-180)
237 
238  // Put lon values outside (-180, 180) into the min or max bins.
239 
240  if(lrel_pos[0] < 0.0)
241  {
242  lrel_pos[0] = 0.01; // Truncates to 0
243  }
244  else if(lrel_pos[0] >= _bin_ub[0])
245  {
246  lrel_pos[0] = _bin_0_max; // Truncates to _bin_ub[0] - 1
247  }
248 
249  lrel_pos[1] = (xpt[1] + 90.0)*_one_over_min_bin_size[1]; // - (-90)
250 
251 #ifdef DIAGNOSTIC_OUTPUT
252 
253  cout << "rel pos: " << setw(12) << lrel_pos[0] << setw(12) << lrel_pos[1] << endl;
254 #endif
255 
256  xresult = lrel_pos;
257 
258  // Postconditions:
259 
260  // Exit:
261 
262  return;
263 }
264 
265 
266 void
269 {
270  // Preconditions:
271 
272  require(xbox != 0);
273 
274  // Body:
275 
276  define_old_variable(size_type old_box_ct = _box_ct);
277 
278  is_abstract();
279 
280  // Postconditions:
281 
282  ensure(contains_box(xbox));
283  ensure(box_ct() == old_box_ct + 1);
284 
285  // Exit:
286 
287  return;
288 }
289 
290 void
293 {
294  // Preconditions:
295 
296  require(xbox != 0);
297  require(contains_box(xbox));
298 
299  // Body:
300 
301  define_old_variable(size_type old_box_ct = _box_ct);
302 
303  is_abstract();
304 
305  // Postconditions:
306 
307  ensure(!contains_box(xbox));
308  ensure(box_ct() == old_box_ct - 1);
309 
310  // Exit:
311 
312  return;
313 }
314 
318 {
319  // Preconditions:
320 
321  require(xpt != 0);
322  require(xpt_ub >= db());
323 
324  // Body:
325 
326 
327  is_abstract();
328 
329  static box_list_type result; // just to avoid compiler warnings
330 
331  // Postconditions:
332 
333 
334  // Exit:
335 
336  return result;
337 }
338 
339 bool
342 {
343  bool result = true;
344 
345  // Preconditions:
346 
347  require(xbox != 0);
348 
349  // Body:
350 
351  is_abstract();
352 
353  // Postconditions:
354 
355 
356  // Exit:
357 
358  return result;
359 }
360 
361 void
364 {
365  // Preconditions:
366 
367  // Body:
368 
369  is_abstract();
370 
371  // Postconditions:
372 
373  ensure(is_empty());
374 
375  // Exit:
376 
377  return;
378 }
379 
380 
381 // PROTECTED MEMBER FUNCTIONS
382 
383 void
386 {
387  // Preconditions:
388 
389  require(is_empty());
390 
391  // Body:
392 
393  is_abstract();
394 
395  // Postconditions:
396 
397  // Exit:
398 
399  return;
400 }
401 
402 void
405  const scoped_index& xbranch_id,
406  const scoped_index& xeval_id)
407 {
408  // Preconditions:
409 
410 
411  // Body:
412 
413  // Vertex 0:
414 
415  bitset<4> neg_x, neg_y;
416  block<sec_vd_dof_type> lat_lon(8);
417  lat_lon.set_ct(8);
418 
419  for(int i=0, j=0, k=0; i<4; i++, j+=2, k+=3)
420  {
421  neg_x[i] = (xdofs[k] < 0.0);
422  neg_y[i] = (xdofs[k+1] < 0.0);
423  xyz_to_lat_lon(&xdofs[k], &lat_lon[j]);
424  }
425 
426 #ifdef DIAGNOSTIC_OUTPUT
427  cout << "xdofs " << xdofs << endl;
428  cout << "neg_x " << neg_x << endl;
429  cout << "neg_y " << neg_y << endl << endl;
430 #endif
431 
432  int neg_x_ct = neg_x.count();
433  int neg_y_ct = neg_y.count();
434  bool on_dateline = (neg_x_ct == 4) && (neg_y_ct != 4) && (neg_y_ct != 0);
435 
436  if(on_dateline)
437  {
438  // Create 2 bounding boxes, one at -180, one at +180.
439 
440  // Create the box at 180:
441 
442  for(int i=0; i<4; ++i)
443  {
444  if(neg_y[i])
445  {
446  lat_lon[i*2] += 360.0;
447  }
448  }
449 
450 #ifdef DIAGNOSTIC_OUTPUT
451  cout << "creating bbox at 180" << endl;
452 #endif
453 
454  make_quad_bounding_box(lat_lon, xbranch_id, xeval_id);
455 
456  // Create the box at -180:
457 
458  for(int i=0; i<4; ++i)
459  {
460  lat_lon[i*2] -= 360.0;
461  }
462 
463 #ifdef DIAGNOSTIC_OUTPUT
464  cout << "creating bbox at -180" << endl;
465 #endif
466 
467  make_quad_bounding_box(lat_lon, xbranch_id, xeval_id);
468 
469  }
470  else
471  {
472  // Create a single bounding box.
473 
474  make_quad_bounding_box(lat_lon, xbranch_id, xeval_id);
475  }
476 
477  // Postconditions:
478 
479 
480  // Exit:
481 
482  return;
483 }
484 
485 void
488  const scoped_index& xbranch_id,
489  const scoped_index& xeval_id)
490 {
491  // Preconditions:
492 
493 
494  // Body:
495 
496  // Vertex 0:
497 
498  bitset<3> neg_x, neg_y;
499  block<sec_vd_dof_type> lat_lon(6);
500  lat_lon.set_ct(6);
501 
502  for(int i=0, j=0, k=0; i<3; i++, j+=2, k+=3)
503  {
504  neg_x[i] = (xdofs[k] < 0.0);
505  neg_y[i] = (xdofs[k+1] < 0.0);
506  xyz_to_lat_lon(&xdofs[k], &lat_lon[j]);
507  }
508 
509 #ifdef DIAGNOSTIC_OUTPUT
510  cout << "xdofs " << xdofs << endl;
511  cout << "neg_x " << neg_x << endl;
512  cout << "neg_y " << neg_y << endl << endl;
513 #endif
514 
515  int neg_x_ct = neg_x.count();
516  int neg_y_ct = neg_y.count();
517  bool on_dateline = (neg_x_ct == 3) && (neg_y_ct != 3) && (neg_y_ct != 0);
518 
519  if(on_dateline)
520  {
521  // Create 2 bounding boxes, one at -180, one at +180.
522 
523  // Create the box at 180:
524 
525  for(int i=0; i<3; ++i)
526  {
527  if(neg_y[i])
528  {
529  lat_lon[i*2] += 360.0;
530  }
531  }
532 
533 #ifdef DIAGNOSTIC_OUTPUT
534  cout << "creating bbox at 180" << endl;
535 #endif
536 
537  make_triangle_bounding_box(lat_lon, xbranch_id, xeval_id);
538 
539  // Create the box at -180:
540 
541  for(int i=0; i<3; ++i)
542  {
543  lat_lon[i*2] -= 360.0;
544  }
545 
546 #ifdef DIAGNOSTIC_OUTPUT
547  cout << "creating bbox at -180" << endl;
548 #endif
549 
550  make_triangle_bounding_box(lat_lon, xbranch_id, xeval_id);
551 
552  }
553  else
554  {
555  // Create a single bounding box.
556 
557  make_triangle_bounding_box(lat_lon, xbranch_id, xeval_id);
558  }
559 
560  // Postconditions:
561 
562 
563  // Exit:
564 
565  return;
566 }
567 
568 void
570 xyz_to_lat_lon(const sec_vd_dof_type xcartesian[3], sec_vd_dof_type xresult[2])
571 {
572  // Preconditions:
573 
574  // Body:
575 
576  static const double RAD2DEG = 180.0/M_PI;
577 
578  double x = xcartesian[0];
579  double y = xcartesian[1];
580  double z = xcartesian[2];
581 
582  // Convert (x,y,z) to phi and theta in radians.
583 
584  double radius = sqrt(x*x + y*y + z*z);
585  double theta = atan2(y, x);
586  double phi = acos(z/radius);
587 
588  double longitude = RAD2DEG*theta;
589  double latitude = 90.0 - RAD2DEG*phi;
590 
591  // Note that we have longitude before latitude.
592 
593  xresult[0] = longitude;
594  xresult[1] = latitude;
595 
596  // Postconditions:
597 
598  // Exit:
599 
600 }
601 
602 void
605  const scoped_index& xbranch_id,
606  const scoped_index& xeval_id)
607 {
608  // Preconditions:
609 
610 
611  // Body:
612 
613 #ifdef DIAGNOSTIC_OUTPUT
614  cout << "lat lon dofs: " << xlat_lon_dofs << endl;
615 #endif
616 
617  static bilinear_2d levaluator;
618 
619  // Get a bounding box.
620 
621  _boxes.set_ct(_boxes.ct() + 1);
622  d_bounding_box<2, 2>& bbox = _boxes.back();
623 
624  // Store the dofs.
625 
627  _gathered_dofs.push_back(xlat_lon_dofs);
628 
629  // Get the bounds for this eval member.
630 
631  sec_vd_value_type lat_lon_min[2];
632  levaluator.min(xlat_lon_dofs.base(), xlat_lon_dofs.ct(), lat_lon_min, 2);
633 
634 #ifdef DIAGNOSTIC_OUTPUT
635 
636  cout << "lat_lon_min: " << setw(12) << lat_lon_min[0] << setw(12) << lat_lon_min[1] << endl;
637 #endif
638 
639  sec_vd_value_type lat_lon_max[2];
640  levaluator.max(xlat_lon_dofs.base(), xlat_lon_dofs.ct(), lat_lon_max, 2);
641 
642 #ifdef DIAGNOSTIC_OUTPUT
643 
644  cout << "lat_lon_max: " << setw(12) << lat_lon_max[0] << setw(12) << lat_lon_max[1] << endl;
645 #endif
646 
647  // Put the bounds in the bounding box.
648 
649  d_bin_coordinates<2, 2> rel_pos;
650  relative_position_pa(lat_lon_min, DB, rel_pos);
651  bbox.put_lb(rel_pos);
652 
653  relative_position_pa(lat_lon_max, DB, rel_pos);
654  bbox.put_ub(rel_pos);
655 
656  // Set the other attributes of the bounding box.
657 
658  bbox.put_member_id(xeval_id);
659  bbox.put_branch_id(xbranch_id);
660  bbox.put_evaluator(&levaluator);
661  bbox.put_dof_ct(xlat_lon_dofs.ct());
662 
663 #ifdef DIAGNOSTIC_OUTPUT
664 
665  cout << "bounding box: " << bbox << endl;
666 #endif
667 
668  // Insert the bounding box into the search structure.
669 
670  insert_box(&bbox);
671 
672  // Postconditions:
673 
674 
675  // Exit:
676 
677  return;
678 }
679 
680 void
683  const scoped_index& xbranch_id,
684  const scoped_index& xeval_id)
685 {
686  // Preconditions:
687 
688 
689  // Body:
690 
691 #ifdef DIAGNOSTIC_OUTPUT
692  cout << "lat lon dofs: " << xlat_lon_dofs << endl;
693 #endif
694 
695  static linear_2d levaluator;
696 
697  // Get a bounding box.
698 
699  _boxes.set_ct(_boxes.ct() + 1);
700  d_bounding_box<2, 2>& bbox = _boxes.back();
701 
702  // Store the dofs.
703 
705  _gathered_dofs.push_back(xlat_lon_dofs);
706 
707  // Get the bounds for this eval member.
708 
709  sec_vd_value_type lat_lon_min[2];
710  levaluator.min(xlat_lon_dofs.base(), xlat_lon_dofs.ct(), lat_lon_min, 2);
711 
712 #ifdef DIAGNOSTIC_OUTPUT
713 
714  cout << "lat_lon_min: " << setw(12) << lat_lon_min[0] << setw(12) << lat_lon_min[1] << endl;
715 #endif
716 
717  sec_vd_value_type lat_lon_max[2];
718  levaluator.max(xlat_lon_dofs.base(), xlat_lon_dofs.ct(), lat_lon_max, 2);
719 
720 #ifdef DIAGNOSTIC_OUTPUT
721 
722  cout << "lat_lon_max: " << setw(12) << lat_lon_max[0] << setw(12) << lat_lon_max[1] << endl;
723 #endif
724 
725  // Put the bounds in the bounding box.
726 
727  d_bin_coordinates<2, 2> rel_pos;
728  relative_position_pa(lat_lon_min, DB, rel_pos);
729  bbox.put_lb(rel_pos);
730 
731  relative_position_pa(lat_lon_max, DB, rel_pos);
732  bbox.put_ub(rel_pos);
733 
734  // Set the other attributes of the bounding box.
735 
736  bbox.put_member_id(xeval_id);
737  bbox.put_branch_id(xbranch_id);
738  bbox.put_evaluator(&levaluator);
739  bbox.put_dof_ct(xlat_lon_dofs.ct());
740 
741 #ifdef DIAGNOSTIC_OUTPUT
742 
743  cout << "bounding box: " << bbox << endl;
744 #endif
745 
746  // Insert the bounding box into the search structure.
747 
748  insert_box(&bbox);
749 
750  // Postconditions:
751 
752 
753  // Exit:
754 
755  return;
756 }
757 
758 
759 // ===========================================================
760 // POINT_LOCATOR FACET
761 // ===========================================================
762 
763 // PUBLIC MEMBER FUNCTIONS
764 
768 operator=(const point_locator& xother)
769 {
770 
771  // Preconditions:
772 
773  // Body:
774 
775  not_implemented();
776 
777  // Postconditions:
778 
779  ensure(invariant());
780 
781  // Exit:
782 
783  return *this;
784 }
785 
790 {
791 
792  // Preconditions:
793 
794  // Body:
795 
796  not_implemented();
797 
798  // Postconditions:
799 
800  ensure(invariant());
801 
802  // Exit:
803 
804  return *this;
805 }
806 
807 bool
809 invariant() const
810 {
811  bool result = true;
812 
813  invariance(coordinates().schema().df() == DC);
814 
815  return result;
816 }
817 
818 void
821 {
822 #ifdef DIAGNOSTIC_OUTPUT
823  post_information_message("Entering update");
824 #endif
825 
826  // Preconditions:
827 
828  require(coordinates().state_is_read_accessible());
829 
830  // Body:
831 
832  update_domain();
833  update_bins();
834 
835  // For each coordinate field evaluator, we need to
836  // create a bounding box and insert the bounding box into the search structure.
837  // It's more efficient if we allocate all the bounding boxes at once.
838  // It's also more efficient if we gather all the dofs for each evaluator
839  // in advance. So find out how many evaluators we have:
840 
841  size_type leval_ct = coordinates().schema().evaluation_ct();
842 
843  // We need at least one bounding box for each evaluation_member,
844  // two for any that straddle the dateline.
845 
846  _boxes.reserve(2*leval_ct);
847  _boxes.set_ct(0);
848 
849  // Start the dof buffer off at a reasonable size.
850  // We don't know how many dofs there are for each element,
851  // but the buffer will resize as necessary.
852 
853  _gathered_dofs.reserve(leval_ct*8);
855 
856  block<sec_vd_dof_type> xyz_dofs(12);
857 
858  // Iterate over the branches.
860 
861  sec_vd lbranch;
862  preorder_iterator lbranch_itr(coordinates(), "jims", DOWN, NOT_STRICT);
863  while(!lbranch_itr.is_done())
864  {
865  // Get a handle to this branch.
866 
867  const scoped_index& lbranch_id = lbranch_itr.index();
868  lbranch.attach_to_state(coordinates().host(), lbranch_id);
869 
870  // Re-anchor the eval itr to this branch.
871 
873  _eval_itr.reset(false);
874 
875  while(!_eval_itr.is_done())
876  {
877  // Gather the dofs for this eval member.
878 
879  xyz_dofs.set_ct(0);
880  _eval_itr.gather_dofs(lbranch, xyz_dofs);
881 
882  // Insert this eval member into the search structure.
883 
885 
886  if(ldisc_ct == 4)
887  {
888  insert_quad(xyz_dofs, lbranch_id, _eval_itr.index());
889  }
890  else if(ldisc_ct == 3)
891  {
892  insert_triangle(xyz_dofs, lbranch_id, _eval_itr.index());
893  }
894  else
895  {
896  post_warning_message("Not a bilinear quad or linear triangle.");
897  }
898 
899  // Move on to the next evaluator.
900 
901  _eval_itr.next();
902  }
903 
904  // Move on to the next branch.
905 
906  lbranch_itr.truncate();
907  }
908 
909  // Clean up.
910 
911  lbranch.detach_from_state();
912 
913  // Postconditions:
914 
915 
916  // Exit:
917 
918 #ifdef DIAGNOSTIC_OUTPUT
919 
920  post_information_message("Leaving update");
921 #endif
922 
923  return;
924 }
925 
926 void
929  size_type xvalue_ub,
930  chart_point& xresult)
931 {
932  require(xvalue != 0);
933  require(xvalue_ub >= dc());
934  require(xresult.db() >= db());
935 
936  // Body:
937 
938  sec_vd_value_type global_coords[DB];
939  xyz_to_lat_lon(xvalue, global_coords);
940 
941  // Initialize the chart id to invalid.
942  // If the inversion succeeds, it will be set to valid.
943 
944  xresult.invalidate();
945 
946  // Find the bounding boxes that may contain the value.
947 
948  d_bin_coordinates<DB, DB> lrel_pos;
949  relative_position_pa(global_coords, DB, lrel_pos);
950 
951  const box_list_type& blist = box_list(global_coords, DB);
952 
953  // Iterate over the bounding boxes.
954 
955  box_list_type::const_iterator iter;
956 
957  for(iter = blist.begin(); !xresult.is_valid() && (iter != blist.end()); ++iter)
958  {
959  const d_bounding_box<DB, DB>* bbox = *iter;
960 
961  if(bbox->contains_point(lrel_pos))
962  {
963  // Found a cell with a bounding box that contains the value;
964  // try inverting the value.
965 
966  section_evaluator* levaluator = bbox->evaluator();
967 
968  size_type dof_ct = bbox->dof_ct();
969  sec_vd_dof_type* dofs = _gathered_dofs.base() + bbox->dofs_index();
970 
971  levaluator->coord_at_value(dofs,
972  dof_ct,
973  global_coords,
974  DB,
975  xresult.local_coords(),
976  xresult.db());
977 
978  if(levaluator->in_standard_domain(xresult.local_coords(), xresult.db()))
979  {
980  // Inversion succeeded.
981  // The coordinates are already stored in xresult_coords;
982  // so just set the chart id.
983 
984  xresult.put_chart_id(bbox->member_id());
985  }
986  else
987  {
988  // Inversion failed.
989 
990  xresult.invalidate();
991  }
992  }
993  }
994 
995 
996  // Postconditions:
997 
998  ensure(unexecutable(xresult.is_valid() ?
999  xresult.chart() contains xvalue at xresult.local_coord() :
1000  no chart contains xvalue));
1001 
1002 
1003  // Exit:
1004 
1005  return;
1006 }
1007 
1008 void
1011  size_type xvalue_ub,
1012  block<chart_point_3d>& xresult)
1013 {
1014 #ifdef DIAGNOSTIC_OUTPUT
1015  post_information_message("Entering all_points_at_value");
1016 #endif
1017 
1018  require(xvalue != 0);
1019  require(xvalue_ub >= dc());
1020  require(db() <= 3);
1021 
1022  // Body:
1023 
1024  sec_vd_value_type global_coords[DB];
1025  xyz_to_lat_lon(xvalue, global_coords);
1026 
1027  define_old_variable(int old_xresult_ct = xresult.ct());
1028 
1029  chart_point_3d lchart_pt(invalid_pod_index(), 0.0, 0.0, 0.0);
1030 
1031  // Find the bounding boxes that may contain the value.
1032 
1033  d_bin_coordinates<DB, DB> lrel_pos;
1034  relative_position_pa(global_coords, DB, lrel_pos);
1035 
1036  const box_list_type& blist = box_list(global_coords, DB);
1037 
1038  // Iterate over the bounding boxes.
1039 
1040  box_list_type::const_iterator iter;
1041 
1042  for(iter = blist.begin(); iter != blist.end(); ++iter)
1043  {
1044  const d_bounding_box<DB, DB>* bbox = *iter;
1045 
1046  if(bbox->contains_point(lrel_pos))
1047  {
1048  // Found a cell with a bounding box that contains the value;
1049  // try inverting the value.
1050 
1051  section_evaluator* levaluator = bbox->evaluator();
1052 
1053  size_type dof_ct = bbox->dof_ct();
1054  sec_vd_dof_type* dofs = _gathered_dofs.base() + bbox->dofs_index();
1055 
1056  levaluator->coord_at_value(dofs,
1057  dof_ct,
1058  global_coords,
1059  DB,
1060  lchart_pt.local_coords(),
1061  lchart_pt.db());
1062 
1063  if(levaluator->in_standard_domain(lchart_pt.local_coords(), lchart_pt.db()))
1064  {
1065  // Inversion succeeded.
1066  // The coordinates are already stored in lchart_pt coords;
1067  // so just set the chart id.
1068 
1069  lchart_pt.put_chart_id(bbox->member_id());
1070  xresult.push_back(lchart_pt);
1071  }
1072  else
1073  {
1074  // Inversion failed; do nothing.
1075  }
1076  }
1077  }
1078 
1079 
1080 #ifdef DIAGNOSTIC_OUTPUT
1081  cout << "result: " << xresult << endl;
1082 #endif
1083 
1084  // Postconditions:
1085 
1086  ensure_for_all(i, 0, xresult.ct(), xresult[i].is_valid());
1087 
1088 
1089  // Exit:
1090 
1091 #ifdef DIAGNOSTIC_OUTPUT
1092 
1093  post_information_message("Leaving all_points_at_value");
1094 #endif
1095 
1096  return;
1097 }
1098 
1099 void
1102  size_type xvalue_ub,
1103  block<branch_point_pair>& xresult)
1104 {
1105 #ifdef DIAGNOSTIC_OUTPUT
1106  post_information_message("Entering branch_points_at_value");
1107 #endif
1108 
1109  require(xvalue != 0);
1110  require(xvalue_ub >= dc());
1111  require(dc() <= 3);
1112 
1113  // Body:
1114 
1115  define_old_variable(int old_xresult_ct = xresult.ct());
1116 
1117  sec_vd_value_type global_coords[DB];
1118  xyz_to_lat_lon(xvalue, global_coords);
1119 
1120  branch_point_pair lbranch_pt;
1121  scoped_index& lbranch_id = lbranch_pt.first;
1122  chart_point_3d& lchart_pt = lbranch_pt.second;
1123 
1124  // Find the bounding boxes that may contain the value.
1125 
1126  d_bin_coordinates<DB, DB> lrel_pos;
1127  relative_position_pa(global_coords, DB, lrel_pos);
1128 
1129  const box_list_type& blist = box_list(global_coords, DB);
1130 
1131  // Iterate over the bounding boxes.
1132 
1133  box_list_type::const_iterator iter;
1134 
1135  for(iter = blist.begin(); iter != blist.end(); ++iter)
1136  {
1137  const d_bounding_box<DB, DB>* bbox = *iter;
1138 
1139 #ifdef DIAGNOSTIC_OUTPUT
1140 
1141  cout << "considering " << *bbox << endl;
1142 #endif
1143 
1144  lbranch_id = bbox->branch_id();
1145  if((_branches.find(lbranch_id) == _branches.end()) &&
1146  (bbox->contains_point(lrel_pos)))
1147  {
1148 
1149 #ifdef DIAGNOSTIC_OUTPUT
1150  cout << "Inverting it ... ";
1151 #endif
1152  // Haven't found a point for this branch yet and we've
1153  // found a cell with a bounding box that contains the value;
1154  // try inverting the value.
1155 
1156  section_evaluator* levaluator = bbox->evaluator();
1157 
1158 
1159  size_type dof_ct = bbox->dof_ct();
1160  sec_vd_dof_type* dofs = _gathered_dofs.base() + bbox->dofs_index();
1161 
1162  levaluator->coord_at_value(dofs,
1163  dof_ct,
1164  global_coords,
1165  DB,
1166  lchart_pt.local_coords(),
1167  lchart_pt.db());
1168 
1169  if(levaluator->in_standard_domain(lchart_pt.local_coords(), lchart_pt.db()))
1170  {
1171  // Inversion succeeded.
1172  // The coordinates are already stored in lchart_pt coords;
1173  // so just set the chart id.
1174 
1175 #ifdef DIAGNOSTIC_OUTPUT
1176  cout << "valid" << endl;
1177 #endif
1178 
1179  lchart_pt.put_chart_id(bbox->member_id());
1180  xresult.push_back(lbranch_pt);
1181 
1182  // We've found a point for this branch.
1183 
1184  _branches.insert(lbranch_id);
1185  }
1186  else
1187  {
1188  // Inversion failed; do nothing.
1189 
1190 #ifdef DIAGNOSTIC_OUTPUT
1191  cout << "invalid" << endl;
1192 #endif
1193 
1194  }
1195  }
1196  else
1197  {
1198 #ifdef DIAGNOSTIC_OUTPUT
1199  cout << "Ignoring it" << endl;
1200 #endif
1201 
1202  }
1203 
1204  }
1205 
1206  // Clean up.
1207 
1208  _branches.clear();
1209 
1210 #ifdef DIAGNOSTIC_OUTPUT
1211 
1212  cout << "result: " << xresult << endl;
1213 #endif
1214 
1215  // Postconditions:
1216 
1217  ensure(xresult.ct() >= old_xresult_ct);
1218  ensure_for_all(i, old_xresult_ct, xresult.ct(),
1219  coordinates().host()->contains_member(xresult[i].first, false));
1220  ensure_for_all(i, old_xresult_ct, xresult.ct(), xresult[i].second.is_valid());
1221 
1222 
1223  // Exit:
1224 
1225 #ifdef DIAGNOSTIC_OUTPUT
1226 
1227  post_information_message("Leaving branch_points_at_value");
1228 #endif
1229 
1230  return;
1231 }
1232 
1233 
1234 // ===========================================================
1235 // NON-MEMBER FUNCTIONS
1236 // ===========================================================
SHEAF_DLL_SPEC void sqrt(const sec_at0 &x0, sec_at0 &xresult, bool xauto_access)
Compute sqrt of x0 (sqrt(x0)) (pre-allocated version).
Definition: sec_at0.cc:1556
virtual void update()
Updates the search structure to the current values of coordinates().
void truncate()
Makes this the next member of the subset which is not less than old this, i.e. the depth-first descen...
block< sec_vd_value_type > _one_over_min_bin_size
Reciprocal of the dimensions of the smallest bins.
size_type ct() const
The number of items currently in use.
sec_ed & coordinates() const
The coordinate section this inverts.
const scoped_index & member_id() const
Index of the evaluation member this bounds.
const block< sec_vd_value_type > & bin_size() const
The dimensions of the smallest bins.
void put_lb(const d_bin_coordinates< DC, DB > &xlb)
Copies the contents of xlb to lb().
size_type box_ct() const
The number of bounding boxes stored in the search structure.
A point in a 3D chart space.
virtual void force_is_done()
Force the iterator to be done.
const block< size_type > & bin_ub() const
The upper bound for the bin index.
virtual int db() const =0
The dimension of this chart.
Definition: chart_point.cc:141
void put_dof_ct(size_type xct)
Sets dof_ct() to xct.
virtual int db() const
The dimension of this chart.
virtual cylindrical_point_locator & operator=(const point_locator &xother)
Assignment operator.
virtual void branch_points_at_value(const sec_vd_value_type *xvalue, size_type xvalue_ub, block< branch_point_pair > &xresult)
Finds one chart point in each branch at which coordinates() has value xvalue and appends them to xres...
std::set< stl_scoped_index<> > _branches
The branches for which a point has already been found for the current evaluation member. Used in all_points_at_value, allocated her to avoid reallocation for each query.
STL namespace.
void reserve(index_type xub)
Makes ub() at least xub; if new storage is allocated, it is uninitialized.
eval_iterator _eval_itr
The evaluator iterator used to populate the search structure; must have same life time as the search ...
point_locator()
Default constructor.
void put_chart_id(pod_index_type xchart)
Sets chart_id() to xchart_id.
Definition: chart_point.cc:103
A point in chart space.
Definition: chart_point.h:52
void update_domain()
Initializes the domain bounds and dimension.
cylindrical_point_locator()
Default constructor; disabled.
host_type * host() const
The poset this is a member of.
Definition: sec_at1.cc:525
block< size_type > _bin_ub
The upper bound for the bin index.
const scoped_index & index() const
The index of the component state this handle is attached to.
virtual void point_at_value(const sec_vd_value_type *xvalue, size_type xvalue_ub, chart_point &xresult)
Finds a chart point at which coordinates() has value xvalue.
virtual bool contains_box(d_bounding_box< 2, 2 > *xbox) const =0
True if xbox is in the box list of some bin.
iterator begin()
Returns an iterator to the first element of the container.
Fixed point relative coordinates for a tree domain.
sec_vd_value_type _bin_0_max
A sec_vd_value_type that truncates to the maximum index for bin 0.
virtual bool in_standard_domain(const dof_type xlocal_coords[], size_type xlocal_coords_ub) const =0
True if the specified local coordinates are in the "standard" domain.
bool is_valid() const
True if this ia a valid point in a chart.
Definition: chart_point.cc:303
virtual const box_list_type & box_list(sec_vd_value_type *xpt, size_type xpt_ub) const =0
The list of bounding boxes which may contain xpt.
A section of a fiber bundle with a d-dimensional Euclidean vector space fiber.
Definition: sec_ed.h:47
void relative_position_pa(sec_vd_value_type *xpt, size_type xpt_ub, d_bin_coordinates< 2, 2 > &xresult) const
The position of xpt relative to the lower bound, in integer coordinates; pre-allocated version...
std::pair< scoped_index, chart_point_3d > branch_point_pair
A point in base space paired with a branch in section space.
Definition: geometry.h:79
void invalidate()
Makes this invalid.
Definition: chart_point.cc:322
const bool NOT_STRICT
Iteration strictness control.
Definition: sheaf.h:102
const bool DOWN
Iteration directions.
Definition: sheaf.h:77
block< sec_vd_value_type > _bin_size
The dimensions of the smallest bins.
bool contains_point(const d_bin_coordinates< DC, DB > &xpt) const
True if this contains point xpt.
virtual coord_type * local_coords()=0
The array of local coordinates.
Definition: chart_point.cc:209
size_type dof_ct() const
The number of dofs in the gathered dofs array for the evaluation member this bounds.
void push_back(const_reference_type item)
Insert item at the end of the items in the auto_block.
virtual dof_type max(const dof_type xdofs[], size_type xdofs_ub) const
The maximum value of the scalar or component section defined by xdofs.
size_type dofs_index() const
Index into the gathered dofs array for the evaluation member this bounds.
pointer_type base() const
The underlying storage array.
void put_schema_anchor(const section_space_schema_member &xschema_anchor)
Set schema_anchor() to the same state as xschema_anchor.
void set_ct(size_type xct)
Sets ct() == xct.
void put_member_id(const scoped_index &xid)
Sets member_id() to xid.
An index within the external ("client") scope of a given id space.
Definition: scoped_index.h:116
block< sec_vd_dof_type > _gathered_dofs
The dofs of gathered by evaluation member.
virtual void remove_box(d_bounding_box< 2, 2 > *xbox)=0
Remove xbox from the search structure.
unsigned long size_type
An unsigned integral type used to represent sizes and capacities.
Definition: sheaf.h:52
A bounding box that can be strung together into a list.
SHEAF_DLL_SPEC void atan2(const sec_at0 &x0, const sec_at0 &x1, sec_at0 &xresult, bool xauto_access)
Compute atan2 of x0/x1 (atan2(x0, x1)) (pre-allocated version).
Definition: sec_at0.cc:1228
virtual void all_points_at_value(const sec_vd_value_type *xvalue, size_type xvalue_ub, block< chart_point_3d > &xresult)
Finds all chart points at which coordinates() has value xvalue and appends them to xresult...
void insert_triangle(block< sec_vd_dof_type > &xdofs, const scoped_index &xbranch_id, const scoped_index &xeval_id)
Creates the bounding box or boxes for the triangle defined by the dofs starting at xdofs and inserts ...
virtual void coord_at_value(const dof_type xdofs[], size_type xdofs_ub, const dof_type xvalue[], size_type xvalue_ub, dof_type xlocal_coords[], size_type xlocal_coords_ub) const =0
The local coordinates of a point at which the field has the value xvalue. The dofs are assumed to be ...
virtual bool contains_member(pod_index_type xmbr_hub_id, bool xauto_access=true) const
True if some version of this poset contains poset member with hub id xmbr_hub_id. ...
d_bin_coordinates< 2, 2 > * relative_position(sec_vd_value_type *xpt, size_type xpt_ub) const
The position of xpt relative to the lower bound, in integer coordinates; auto-allocated version...
void put_dofs_index(size_type xindex)
Sets dofs_index() to xindex.
void make_triangle_bounding_box(block< sec_vd_dof_type > &xlat_lon_dofs, const scoped_index &xbranch_id, const scoped_index &xeval_id)
Makes a bounding box for a triangle.
void next()
Makes this the next member of the subset.
void xyz_to_lat_lon(const sec_vd_dof_type xcartesian[3], sec_vd_dof_type xresult[2])
Converts cartesian coordinates xcartesian to (longitude, latitude) xresult.
iterator end()
Returns an iterator to the element following the last element of the container.
An abstract local section evaluator; a map from {local coordinates x dofs} to section value...
coord_type * local_coords()
The array of local coordinates.
A section of a fiber bundle with a d-dimensional vector space fiber.
Definition: sec_vd.h:54
section_evaluator * evaluator() const
Evaluator for the evaluation member this bounds.
virtual void insert_box(d_bounding_box< 2, 2 > *xbox)=0
Insert xbox into the search structure.
bool is_done() const
True if iteration finished.
An abstract point location query in domains with global coordinate dimension dc and local coordinate ...
void detach_from_state()
Detaches field from state it is currently attached to.
Definition: sec_tuple.cc:603
size_type _box_ct
The number of bounding boxes stored in the search structure.
void put_evaluator(section_evaluator *xevaluator)
Sets evaluator() to xevaluator.
int db() const
The intrinsic dimension of the domain; the dimension of the local coordinates.
virtual void reset(bool xreset_markers=true)
Restarts the iteration over the down set of anchor().
virtual dof_type min(const dof_type xdofs[], size_type xdofs_ub) const
The minimum value of the scalar or component section defined by xdofs.
void put_ub(const d_bin_coordinates< DC, DB > &xub)
Copies the contents of xub to ub().
virtual section_space_schema_member & schema()
The restricted schema for this (mutable version).
A section evaluator using linear interpolation over a triangular 2D domain.
Definition: linear_2d.h:39
SHEAF_DLL_SPEC void acos(const sec_at0 &x0, sec_at0 &xresult, bool xauto_access)
Compute acos of x0 (acos(x0)) (pre-allocated version).
Definition: sec_at0.cc:1154
int dc() const
The spatial dimension of the domain; the dimension of the global coordinates.
const scoped_index & branch_id() const
Index of the branch that contains the evaluation member this bounds.
int evaluation_ct() const
The number of members in the intersection of the evaluation subposet and the down set of the base spa...
void gather_dofs(const sec_vd &xsec, block< sec_vd::dof_type > &xdofs)
Gathers the dofs for the current evalaution member from section xsec and appends them to the back of ...
virtual void update_bins()=0
Updates the bin parameters.
block< d_bounding_box< 2, 2 > > _boxes
Bounding boxes for the evaluation members.
bool is_empty() const
True if this contains no bounding boxes.
void attach_to_state(const namespace_poset *xns, const poset_path &xpath, bool xauto_access=true)
Attach to the state specified by path xpath in the namespace xns.
void put_branch_id(const scoped_index &xid)
Sets branch_id() to xid.
double sec_vd_dof_type
The type of degree of freedom in the section space.
Definition: fiber_bundle.h:78
block< scoped_index > & discretization_members()
The discretization members in the downset of the current evaluation member (mutable version)...
An auto_block with a no-initialization initialization policy.
Namespace for geometry component of sheaf system.
Definition: field_vd.h:54
A section evaluator using bilinear interpolation over a square 2D domain.
Definition: bilinear_2d.h:42
vd_value_type sec_vd_value_type
The type of component in the value of a section at a point.
Definition: fiber_bundle.h:73
SHEAF_DLL_SPEC pod_index_type invalid_pod_index()
The invalid pod index value.
Definition: pod_types.cc:31
virtual bool invariant() const
Class invariant.
virtual coord_type local_coord(int xi) const =0
The xi-th local coordinate of this point.
Definition: chart_point.cc:163
void insert_quad(block< sec_vd_dof_type > &xdofs, const scoped_index &xbranch_id, const scoped_index &xeval_id)
Creates the bounding box or boxes for the quad defined by the dofs starting at xdofs and inserts them...
void make_quad_bounding_box(block< sec_vd_dof_type > &xlat_lon_dofs, const scoped_index &xbranch_id, const scoped_index &xeval_id)
Makes a bounding box for a quad.
Wrapper class for forward_list or slist depending on compiler. The class replicates the minimum subse...
An abstract point location query in domains with global coordinate dimension dc and local coordinate ...
Definition: point_locator.h:52
virtual void clear()=0
Clear the search structure of all bounding boxes.
const scoped_index & index() const
The index of the current member of the iteration.