Line data Source code
1 : #include <common/util/string_util.h>
2 : #include "lib_grid/algorithms/space_partitioning/lg_ntree.h"
3 : #include <lib_grid/grid_objects/grid_dim_traits.h>
4 : #ifdef UG_PARALLEL
5 : #include "lib_grid/parallelization/gather_grid.h"
6 : #include "lib_grid/parallelization/broadcast.h"
7 : #include "lib_grid/parallelization/distributed_grid.h"
8 : #endif
9 : #include "lib_grid/file_io/file_io.h"
10 :
11 : namespace ug {
12 :
13 : template <typename TDomain>
14 : class OverlyingSubsetFinder
15 : {
16 : static const int dim = TDomain::dim;
17 : typedef TDomain domain_type;
18 :
19 : typedef typename grid_dim_traits<dim-1>::element_type side_t;
20 :
21 : typedef lg_ntree<dim-1, dim, side_t> top_tracer_tree_t;
22 :
23 : typedef RayElemIntersectionRecord<side_t*> top_intersection_record_t;
24 :
25 : typedef typename domain_type::position_attachment_type position_attachment_type;
26 : typedef typename domain_type::position_accessor_type position_accessor_type;
27 :
28 : typedef typename Grid::traits<side_t>::iterator SideIterator;
29 :
30 : public:
31 :
32 0 : OverlyingSubsetFinder(
33 : SmartPtr<TDomain> sp_domain, ///< the domain
34 : const std::string & top_ss_names) ///< top surface subset names
35 : : m_sp_domain (sp_domain),
36 0 : m_top_ss_grp (sp_domain->subset_handler (), TokenizeString (top_ss_names))
37 : {
38 :
39 0 : MultiGrid & mg = * m_sp_domain->grid ();
40 0 : MGSubsetHandler & sh = * m_sp_domain->subset_handler ();
41 : std::vector<side_t*> topSides;
42 :
43 0 : m_top_tracer_tree.set_grid(m_top_grid, m_sp_domain->position_attachment ());
44 :
45 0 : mg.attach_to<side_t>(subset_index_attachment);
46 0 : m_top_grid.attach_to<side_t>(subset_index_attachment);
47 :
48 : Grid::AttachmentAccessor<side_t, Attachment<int> > subset_index_accessor;
49 0 : subset_index_accessor.access(mg, subset_index_attachment);
50 :
51 : // select the Sides on the top
52 0 : Selector sel (mg);
53 0 : for (size_t i = 0; i < m_top_ss_grp.size (); i++)
54 : {
55 : int si = m_top_ss_grp [i];
56 :
57 0 : for (SideIterator it = sh.begin<side_t> (si, 0); it != sh.end<side_t> (si, 0); ++it)
58 : {
59 : sel.select (*it);
60 0 : subset_index_accessor[*it] = si;
61 : }
62 : }
63 :
64 : position_accessor_type accessor;
65 : accessor.access(mg, m_sp_domain->position_attachment());
66 :
67 :
68 : //UG_LOG_ALL_PROCS("Total count of Edges on this node before distributing: " << sel.num<side_t>() << std::endl);
69 :
70 : // for (SideIterator it = sel.begin<side_t>(); it != sel.end<side_t>(); ++it)
71 : // {
72 : // UG_LOG_ALL_PROCS(subset_index_accessor[(*it)] << ": " << accessor[(*it)->vertex(0)] << "-" << accessor[(*it)->vertex(1)] << std::endl);
73 : // }
74 :
75 : // copy the top Sides into a new grid
76 :
77 : #ifdef UG_PARALLEL
78 :
79 : GridDataSerializationHandler serializer;
80 : serializer.add
81 : (GeomObjAttachmentSerializer<Vertex, position_attachment_type>::create (mg, m_sp_domain->position_attachment ()));
82 : serializer.add
83 : (GeomObjAttachmentSerializer<side_t, Attachment<int> >::create (mg, subset_index_attachment));
84 :
85 : GridDataSerializationHandler deserializer;
86 : deserializer.add
87 : (GeomObjAttachmentSerializer<Vertex, position_attachment_type>::create (m_top_grid, m_sp_domain->position_attachment ()));
88 : deserializer.add
89 : (GeomObjAttachmentSerializer<side_t, Attachment<int> >::create (m_top_grid, subset_index_attachment));
90 :
91 : BroadcastGrid (m_top_grid, sel, serializer, deserializer, 0);
92 :
93 : #endif
94 :
95 0 : for (SideIterator it = m_top_grid.begin<side_t> (); it != m_top_grid.end<side_t> (); ++it)
96 0 : topSides.push_back (*it);
97 :
98 0 : m_top_tracer_tree.create_tree (topSides.begin (), topSides.end ());
99 :
100 : // debug();
101 :
102 : // SaveGridToFile(m_top_grid, mkstr("top_grid_p" << pcl::ProcRank() << ".ugx").c_str(),
103 : // m_sp_domain->position_attachment());
104 0 : }
105 :
106 0 : std::string findOverlyingSubset(const std::vector<number> & point)
107 : {
108 :
109 0 : if (point.size () != dim)
110 0 : UG_THROW ("OverlyingSubsetFinder::findOverlyingSubset: Exactly " << dim << " coordinates have to be specified!");
111 :
112 : MathVector<dim> measure_point;
113 :
114 0 : for(size_t d = 0; d < dim; d++)
115 : {
116 0 : measure_point[d] = point[d];
117 : }
118 :
119 : // find all the intersections
120 : MathVector<dim> up_dir;
121 : up_dir = 0;
122 0 : up_dir [dim - 1] = 1;
123 : m_top_intersection_records.clear ();
124 0 : RayElementIntersections (m_top_intersection_records, m_top_tracer_tree, measure_point, up_dir, 0.001);
125 :
126 : // check if there are intersections at all
127 0 : if (m_top_intersection_records.size () == 0)
128 0 : return "";
129 :
130 : Grid::AttachmentAccessor<side_t, Attachment<int> > subset_index_accessor;
131 0 : subset_index_accessor.access(m_top_grid, subset_index_attachment);
132 :
133 : //debug();
134 :
135 : // find the lowest point
136 0 : MathVector<dim> x = PointOnRay (measure_point, up_dir, m_top_intersection_records[0].smin);
137 0 : number z_min = x [dim - 1];
138 0 : int si = subset_index_accessor[m_top_intersection_records[0].elem];
139 :
140 : // UG_LOG("<dbg> x: " << x << std::endl);
141 0 : for (size_t i = 1; i < m_top_intersection_records.size (); i++)
142 : {
143 : top_intersection_record_t & r = m_top_intersection_records [i];
144 0 : x = PointOnRay (measure_point, up_dir, r.smin);
145 0 : if (x [dim - 1] < z_min)
146 : {
147 : z_min = x [dim - 1];
148 0 : si = subset_index_accessor[r.elem];
149 : }
150 : }
151 :
152 0 : return m_sp_domain->subset_handler ()->get_subset_name(si);
153 : }
154 :
155 : void debug()
156 : {
157 : position_accessor_type accessor;
158 : accessor.access(m_top_grid, m_sp_domain->position_attachment());
159 :
160 : Grid::AttachmentAccessor<side_t, Attachment<int> > subset_index_accessor;
161 : subset_index_accessor.access(m_top_grid, subset_index_attachment);
162 :
163 : UG_LOG_ALL_PROCS("Total count of Edges: " << m_top_grid.num<side_t>() << std::endl);
164 :
165 : for (SideIterator it = m_top_grid.begin<side_t>(); it != m_top_grid.end<side_t>(); ++it)
166 : {
167 : UG_LOG_ALL_PROCS(subset_index_accessor[(*it)] << ": ");
168 : for (size_t i = 0; i < dim; i++)
169 : {
170 : UG_LOG_ALL_PROCS(accessor[(*it)->vertex(i)] << "-");
171 : }
172 : UG_LOG_ALL_PROCS(std::endl);
173 : }
174 :
175 : }
176 :
177 : private:
178 : SmartPtr<domain_type> m_sp_domain;
179 :
180 : SubsetGroup m_top_ss_grp;
181 : Grid m_top_grid;
182 : Attachment<int> subset_index_attachment;
183 :
184 : top_tracer_tree_t m_top_tracer_tree;
185 :
186 : /// array to store all the intersections
187 : std::vector<top_intersection_record_t> m_top_intersection_records;
188 :
189 : };
190 :
191 : }
|
