Line data Source code
1 : /*
2 : * Copyright (c) 2013-2015: G-CSC, Goethe University Frankfurt
3 : * Author: Sebastian Reiter
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 : #ifndef __H__UG__lg_ntree__
34 : #define __H__UG__lg_ntree__
35 :
36 : #include "common/space_partitioning/ntree.h"
37 : #include "common/space_partitioning/ntree_traverser.h"
38 : #include "common/math/misc/shapes.h"
39 : #include "lib_grid/algorithms/ray_element_intersection_util.h"
40 : #include "lib_grid/algorithms/geom_obj_util/geom_obj_util.h"
41 :
42 : namespace ug{
43 :
44 :
45 : template <int world_dim>
46 : class NTreeGridData
47 : {
48 : public:
49 : typedef MathVector<world_dim > position_t;
50 : typedef Attachment<position_t> position_attachment_t;
51 : typedef Grid::VertexAttachmentAccessor<position_attachment_t> position_accessor_t;
52 :
53 0 : NTreeGridData() : m_pGrid(NULL) {}
54 :
55 0 : NTreeGridData(Grid& grid, position_attachment_t aPos)
56 : {
57 0 : m_pGrid = &grid;
58 0 : if(!grid.has_vertex_attachment(aPos))
59 0 : grid.attach_to_vertices(aPos);
60 : m_aaPos.access(grid, aPos);
61 0 : }
62 :
63 : const position_t& position(Vertex* v) const
64 : {
65 : UG_ASSERT(m_aaPos.valid(),
66 : "Make sure to pass an instance of NTreeGridData to lg_ntree::set_common_data");
67 : return m_aaPos[v];
68 : }
69 :
70 : position_accessor_t position_accessor() const {return m_aaPos;}
71 :
72 0 : Grid* grid_ptr() const {return m_pGrid;}
73 :
74 : private:
75 : position_accessor_t m_aaPos;
76 : Grid* m_pGrid;
77 :
78 : };
79 :
80 :
81 : template <int tree_dim, int world_dim, class elem_t_, class common_data_t_>
82 : struct lg_ntree_traits_base
83 : {
84 : typedef number real_t;
85 : typedef MathVector<world_dim> vector_t;
86 : typedef AABox<vector_t> box_t;
87 : typedef common_data_t_ common_data_t;
88 : typedef elem_t_ elem_t;
89 :
90 0 : static void calculate_center(vector_t& centerOut, const elem_t& e,
91 : const common_data_t& commonData)
92 : {
93 : Grid::vertex_traits::secure_container vrts;
94 0 : commonData.grid_ptr()->associated_elements(vrts, e);
95 :
96 : assert(vrts.size() > 0);
97 : centerOut = commonData.position(vrts[0]);
98 0 : for(size_t i = 1; i < vrts.size(); ++i)
99 : VecAdd(centerOut, centerOut, commonData.position(vrts[i]));
100 0 : VecScale(centerOut, centerOut, 1. / (real_t)vrts.size());
101 0 : }
102 :
103 0 : static void calculate_bounding_box(box_t& boxOut, const elem_t& e,
104 : const common_data_t& commonData)
105 : {
106 : Grid::vertex_traits::secure_container vrts;
107 0 : commonData.grid_ptr()->associated_elements(vrts, e);
108 :
109 : assert(vrts.size() > 0);
110 : boxOut.min = boxOut.max = commonData.position(vrts[0]);
111 0 : for(size_t i = 1; i < vrts.size(); ++i)
112 0 : boxOut = box_t(boxOut, commonData.position(vrts[i]));
113 0 : }
114 :
115 : static void grow_box(box_t& boxOut, const box_t& box,
116 : const vector_t& offset)
117 : {
118 : VecSubtract(boxOut.min, box.min, offset);
119 : VecAdd(boxOut.max, box.max, offset);
120 : }
121 :
122 : static vector_t box_diagonal(const box_t& box)
123 : {
124 : vector_t d;
125 : VecSubtract(d, box.max, box.min);
126 : return d;
127 : }
128 :
129 : static bool box_contains_point(const box_t& box, const vector_t& point)
130 : {
131 : return box.contains_point(point);
132 : }
133 :
134 : /// returns true if the given boxes intersect
135 : static bool box_box_intersection(const box_t& box1, const box_t& box2)
136 : {
137 : for(int i = 0; i < world_dim; ++i){
138 : if(box1.min[i] > box2.max[i] || box1.max[i] < box2.min[i])
139 : return false;
140 : }
141 : return true;
142 : }
143 :
144 : static bool ray_box_intersection(const vector_t& from,
145 : const vector_t& dir,
146 : const box_t& box)
147 : {
148 0 : return RayBoxIntersection(from, dir, box.min, box.max);
149 : }
150 :
151 : /// returns the smallest box that contains both box1 and box2
152 0 : static void merge_boxes(box_t& boxOut, const box_t& box1, const box_t& box2)
153 : {
154 0 : boxOut = box_t(box1, box2);
155 0 : }
156 :
157 0 : static bool contains_point(const elem_t& e, const vector_t& point,
158 : const common_data_t& commonData)
159 : {
160 0 : return ContainsPoint(e, point, commonData.position_accessor());
161 :
162 : // todo: think about some fast-checks like the following (but faster!!!)
163 : // first check whether the bounding box contains the point, then perform
164 : // the exact check (slow e.g. for tetrahedrons...)
165 : // typename common_data_t::position_accessor_t aaPos = commonData.position_accessor();
166 : // box_t box(aaPos[e->vertex(0)], aaPos[e->vertex(1)]);
167 : // for(size_t i = 2; i < e->num_vertices(); ++i)
168 : // box = box_t(box, aaPos[e->vertex(i)]);
169 : //
170 : // if(box.contains_point(point))
171 : // return ContainsPoint(e, point, aaPos);
172 : // return false;
173 : }
174 :
175 :
176 : static bool intersects_ray( Vertex* e,
177 : const vector_t& rayFrom,
178 : const vector_t& rayDir,
179 : const common_data_t& cd,
180 : number& s0out,
181 : number& s1out)
182 : {
183 : UG_THROW("intersects_ray not yet implemented for vertices");
184 : return false;
185 : }
186 :
187 0 : static bool intersects_ray( Edge* e,
188 : const vector_t& rayFrom,
189 : const vector_t& rayDir,
190 : const common_data_t& cd,
191 : number& s0out,
192 : number& s1out,
193 : number sml = SMALL)
194 : {
195 0 : UG_COND_THROW(!cd.grid_ptr(), "No grid assigned to ntree::common_data.");
196 0 : return RayElementIntersection(s0out, s1out, rayFrom, rayDir,
197 : e, *cd.grid_ptr(), cd.position_accessor(),
198 0 : sml);
199 : }
200 :
201 0 : static bool intersects_ray( Face* e,
202 : const vector_t& rayFrom,
203 : const vector_t& rayDir,
204 : const common_data_t& cd,
205 : number& s0out,
206 : number& s1out,
207 : number sml = SMALL)
208 : {
209 0 : UG_COND_THROW(!cd.grid_ptr(), "No grid assigned to ntree::common_data.");
210 0 : return RayElementIntersection(s0out, s1out, rayFrom, rayDir,
211 : e, *cd.grid_ptr(), cd.position_accessor(),
212 0 : sml);
213 : }
214 :
215 : static bool intersects_ray( Volume* e,
216 : const vector_t& rayFrom,
217 : const vector_t& rayDir,
218 : const common_data_t& cd,
219 : number& s0out,
220 : number& s1out,
221 : number sml = SMALL)
222 : {
223 : UG_COND_THROW(!cd.grid_ptr(), "No grid assigned to ntree::common_data.");
224 : return RayElementIntersection(s0out, s1out, rayFrom, rayDir,
225 : e, *cd.grid_ptr(), cd.position_accessor(),
226 : sml);
227 : }
228 : };
229 :
230 :
231 : template <class elem_t>
232 : struct ntree_traits<1, 1, elem_t, NTreeGridData<1> > :
233 : public lg_ntree_traits_base<1, 1, elem_t, NTreeGridData<1> >
234 : {
235 : typedef MathVector<1> vector_t;
236 : typedef AABox<vector_t> box_t;
237 :
238 : static void split_box(box_t* boxesOut, const box_t& box, const vector_t& splitPoint)
239 : {
240 : boxesOut[0] = box_t(box.min, splitPoint);
241 : boxesOut[1] = box_t(splitPoint, box.max);
242 : }
243 : };
244 :
245 :
246 : template <class elem_t>
247 : struct ntree_traits<1, 2, elem_t, NTreeGridData<2> > :
248 : public lg_ntree_traits_base<1, 2, elem_t, NTreeGridData<2> >
249 : {
250 : typedef MathVector<2> vector_t;
251 : typedef AABox<vector_t> box_t;
252 :
253 : static void split_box(box_t* boxesOut, const box_t& box, const vector_t& splitPoint)
254 : {
255 : vector_t splitPointYMin = splitPoint;
256 0 : splitPointYMin.y() = box.min.y();
257 : vector_t splitPointYMax = splitPoint;
258 0 : splitPointYMax.y() = box.max.y();
259 : boxesOut[0] = box_t(box.min, splitPointYMax);
260 : boxesOut[1] = box_t(splitPointYMin, box.max);
261 : }
262 : };
263 :
264 :
265 : template <class elem_t>
266 : struct ntree_traits<2, 2, elem_t, NTreeGridData<2> > :
267 : public lg_ntree_traits_base<2, 2, elem_t, NTreeGridData<2> >
268 : {
269 : typedef MathVector<2> vector_t;
270 : typedef AABox<vector_t> box_t;
271 :
272 : static void split_box(box_t* boxesOut, const box_t& box, const vector_t& splitPoint)
273 : {
274 : boxesOut[0] = box_t(box.min, splitPoint);
275 : boxesOut[1] = box_t(vector_t(splitPoint.x(), box.min.y()),
276 0 : vector_t(box.max.x(), splitPoint.y()));
277 : boxesOut[2] = box_t(vector_t(box.min.x(), splitPoint.y()),
278 0 : vector_t(splitPoint.x(), box.max.y()));
279 : boxesOut[3] = box_t(splitPoint, box.max);
280 : }
281 : };
282 :
283 :
284 : template <class elem_t>
285 : struct ntree_traits<2, 3, elem_t, NTreeGridData<3> > :
286 : public lg_ntree_traits_base<2, 3, elem_t, NTreeGridData<3> >
287 : {
288 : typedef MathVector<3> vector_t;
289 : typedef AABox<vector_t> box_t;
290 :
291 0 : static void split_box(box_t* boxesOut, const box_t& box, const vector_t& splitPoint)
292 : {
293 : vector_t splitPointZMin = splitPoint;
294 0 : splitPointZMin.z() = box.min.z();
295 : vector_t splitPointZMax = splitPoint;
296 0 : splitPointZMax.z() = box.max.z();
297 :
298 : boxesOut[0] = box_t(box.min, splitPointZMax);
299 0 : boxesOut[1] = box_t(vector_t(splitPoint.x(), box.min.y(), box.min.z()),
300 0 : vector_t(box.max.x(), splitPoint.y(), box.max.z()));
301 0 : boxesOut[2] = box_t(vector_t(box.min.x(), splitPoint.y(), box.min.z()),
302 0 : vector_t(splitPoint.x(), box.max.y(), box.max.z()));
303 : boxesOut[3] = box_t(splitPointZMin, box.max);
304 0 : }
305 : };
306 :
307 : template <class elem_t>
308 : struct ntree_traits<3, 3, elem_t, NTreeGridData<3> > :
309 : public lg_ntree_traits_base<3, 3, elem_t, NTreeGridData<3> >
310 : {
311 : typedef MathVector<3> vector_t;
312 : typedef AABox<vector_t> box_t;
313 :
314 0 : static void split_box(box_t* boxesOut, const box_t& box, const vector_t& splitPoint)
315 : {
316 : boxesOut[0] = box_t(box.min, splitPoint);
317 0 : boxesOut[1] = box_t(vector_t(splitPoint.x(), box.min.y(), box.min.z()),
318 0 : vector_t(box.max.x(), splitPoint.y(), splitPoint.z()));
319 0 : boxesOut[2] = box_t(vector_t(box.min.x(), splitPoint.y(), box.min.z()),
320 0 : vector_t(splitPoint.x(), box.max.y(), splitPoint.z()));
321 0 : boxesOut[3] = box_t(vector_t(splitPoint.x(), splitPoint.y(), box.min.z()),
322 0 : vector_t(box.max.x(), box.max.y(), splitPoint.z()));
323 :
324 0 : boxesOut[4] = box_t(vector_t(box.min.x(), box.min.y(), splitPoint.z()),
325 0 : vector_t(splitPoint.x(), splitPoint.y(), box.max.z()));
326 0 : boxesOut[5] = box_t(vector_t(splitPoint.x(), box.min.y(), splitPoint.z()),
327 0 : vector_t(box.max.x(), splitPoint.y(), box.max.z()));
328 0 : boxesOut[6] = box_t(vector_t(box.min.x(), splitPoint.y(), splitPoint.z()),
329 0 : vector_t(splitPoint.x(), box.max.y(), box.max.z()));
330 : boxesOut[7] = box_t(splitPoint, box.max);
331 0 : }
332 : };
333 :
334 :
335 :
336 : template <int tree_dim, int world_dim, class grid_elem_t>
337 0 : class lg_ntree : public ntree<tree_dim, world_dim, grid_elem_t*, NTreeGridData<world_dim> >
338 : {
339 : public:
340 : typedef ntree<tree_dim, world_dim, grid_elem_t*, NTreeGridData<world_dim> > base_t;
341 : typedef typename NTreeGridData<world_dim>::position_attachment_t position_attachment_t;
342 :
343 0 : lg_ntree()
344 0 : {}
345 :
346 0 : lg_ntree(Grid& grid, position_attachment_t aPos) :
347 0 : m_gridData(grid, aPos)
348 0 : {}
349 :
350 0 : void set_grid(Grid& grid, position_attachment_t aPos)
351 : {
352 0 : m_gridData = NTreeGridData<world_dim>(grid, aPos);
353 0 : }
354 :
355 : template <class TIterator>
356 0 : void create_tree(TIterator elemsBegin, TIterator elemsEnd)
357 : {
358 : base_t::set_common_data(m_gridData);
359 :
360 0 : base_t::clear();
361 :
362 0 : while(elemsBegin != elemsEnd){
363 : base_t::add_element(*elemsBegin);
364 0 : ++elemsBegin;
365 : }
366 :
367 0 : base_t::rebalance();
368 0 : }
369 :
370 : private:
371 : NTreeGridData<world_dim> m_gridData;
372 : };
373 :
374 :
375 : }// end of namespace
376 :
377 : #endif
|
