rbmatlab  1.16.09
 All Classes Namespaces Files Functions Variables Modules Pages
XPartMap.m
1 classdef XPartMap < handle
2  % a geometric tree where every nodes has geometry information attached.
3  %
4  % The geometry relations are as follows:
5  % - Child node''s geometries lie in their parents geometries.
6  % - Sibling''s geometries have zero intersection.
7  % - The union of all child node''s geometries equals a parent''s geometry.
8  % - The geometries are of arbitrary dimension and can be defined by an
9  % arbitrary number of coordinates.
10  %
11  % The main features of this class are:
12  % - Routines to split the geometries and to create new childs are provided.
13  % - The leaf nodes can be accessed over a sorted or an unsorted index.
14  % - We can check in `{\cal O}(\log)`-time in which node''s geometry a certain
15  % coordinate lies.
16  %
17  % @section Implementation details
18  % The nodes store the following information:
19  % -# a number with the ID of the first child of this node or '0'
20  % indicating a leaf node
21  % -# the 'ncoords' coordinates describing the node''s geometry.
22  %
23  % So every node has a size of '1 + ncoords * dimension'.
24  %
25  % We distinguish between three enumerations of nodes:
26  % -# a "node ID" is given to every node in the tree and does not change
27  % during a refinement for nodes that are unchanged.
28  % -# a "leaf element index" is a consecutive enumeration of leaf
29  % indices which does not change for unchanged leaf nodes during a
30  % refinement
31  % -# a "sorted leaf element index" is a consecutive enumeration of leaf
32  % indices enumerating the leaf elements from left-to-right.
33  %
34 
35  properties(Access=public)
36 
37  dimension = 1; % dimension of geometries
38 
39  ncoords = 2; % number of coordinates defining a geometry
40 
41  % actual tree with attached geometries
42  %
43  % This stores a list of pairs '(c, coords)' of the index of the first child
44  % 'c' and the geometry coordinates 'coords'. If 'c' is equal to 0, the node
45  % is a leaf node. The first pair is the root node of the tree.
46  intervals = [0 0 1];
47 
48  % integer specifying how many children are created at refinement step
49  split_size = 2;
50 
51  % an 'unsorted' enumeration of the leaf indices
52  %
53  % If a leaf has not been refined, it keeps its position in this set.
54  leaf_ids = 0;
55 
56  % a 'sorted' enumeration of the leaf indices.
57  %
58  % This set is constructed by a deep traversal of the tree, with
59  % left-to-right enumeration of leafs.
60  leaf_enum = 0;
61 
62  % A vector which stores the refinement level for every node.
63  refine_levels = 0;
64 
65  end
66 
67  properties(Dependent)
68  % the number of the leaf nodes
69  siz;
70 
71  % storage size of a single tree node pair '(c, coords)'
72  skip;
73 
74  % the maximum refinement level of the tree
75  max_level;
76 
77  noof_ids; % the number of node IDs in the tree
78  end
79 
80  methods
81 
82  function skip = get.skip(this)
83  skip = this.dimension*this.ncoords+1;
84  end
85 
86  function siz = get.siz(this)
87  siz = length(this.leaf_ids);
88  end
89 
90  function nids = get.noof_ids(this)
91  nids = length(this.intervals)/this.skip;
92  end
93 
94  function mlevel =get.max_level(this)
95  mlevel = max(this.refine_levels);
96  end
97 
98  function copy = copy(this)
99  % function copy = copy(this)
100  % deep copy for an XPartMap instance
101  %
102  % Return values:
103  % copy: the copied object of type XPartMap
104 
105  copy = XPartMap(this.dimension);
106 
107  fns = setdiff(properties(this), {'siz', 'skip', 'max_level', 'noof_ids'});
108 
109  for i = 1:length(fns)
110  copy.(fns{i}) = this.(fns{i});
111  end
112  end
113 
114  function tpm = XPartMap(dimension, ncoords, init_coords, level1_splits)
115  % function tpm = XPartMap(dimension, ncoords, init_coords)
116  % constructor creating an XPartMap with one single root node.
117  %
118  % Parameters:
119  % dimension: @copybrief #dimension
120  % ncoords: @copybrief #ncoords
121  % init_coords: initial coordinates for the root node''s geometry.
122  if nargin >=1
123  tpm.dimension = dimension;
124  tpm.split_size = 2^dimension;
125  end
126  if nargin >=2
127  tpm.ncoords = ncoords;
128  if tpm.ncoords > 2 && nargin == 2
129  error('if ncoords > 2, init_coords need to be specified');
130  end
131  end
132  if nargin >=3
133  tpm.intervals = [0, init_coords];
134  else
135  icoords = zeros(1, tpm.ncoords * tpm.dimension);
136  icoords(tpm.dimension+1:tpm.dimension*2) = 1;
137  tpm.intervals = [0 icoords];
138  end
139  if nargin >=4
140  if tpm.ncoords ~= 2
141  error('level1 splits are only supported for ncoords == 2');
142  end
143  assert(size(level1_splits, 2) == tpm.dimension);
144  num_splits = size(level1_splits, 1) - 1;
145  tpm.intervals = [tpm.intervals, zeros(1, tpm.skip * num_splits)];
146  tpm.refine_levels = [tpm.refine_levels, ones(1, num_splits)];
147  tpm.leaf_enum = 1:num_splits;
148  tpm.leaf_ids = 1:num_splits;
149  tpm.intervals(1) = 1;
150  left = level1_splits(1, :);
151 
152  for i=2:(num_splits+1)
153  right = level1_splits(i, :);
154  tpm.intervals((i-1) * tpm.skip + 2:i*tpm.skip) = [left, right];
155  left = right;
156  end
157 
158  end
159  end
160 
161  function elems = coords2elem(this, tcoords)
162  % function elems = coords2elem(this, tcoords)
163  % returns element IDs of leaf nodes whose geometries host the given
164  % coordinates.
165  %
166  % Parameters:
167  % tcoords: 'n x dimension' matrix of 'n' coordinate (row) vectors for
168  % which the enclosing entity ID shall be determined.
169  %
170  % Return values:
171  % elems: a 'n' vector holding the element IDs.
172 
173  elems = -1 * ones(1, size(tcoords,1));
174  if this.dimension > 1 || size(tcoords,1) * this.max_level < this.siz
175 
176  for i = 1:size(tcoords,1)
177  tc = tcoords(i,:);
178  j=0;
179  while true
180  if this.is_inside(tc, j)
181  if this.intervals(j*this.skip+1) == 0
182  break
183  else
184  j = this.intervals(j*this.skip+1);
185  end
186  else
187  j = j + 1;
188  if j >= this.noof_ids
189  j = -1;
190  break;
191  end
192  end
193  end
194  elems(i) = j;
195  end
196 
197  else
198  [stcoords, map] = sort(tcoords);
199  i = 1;
200  tc = stcoords(i);
201  for j = this.leaf_enum
202  while this.is_inside(tc, j)
203  elems(map(i)) = j;
204  i = i + 1;
205  if i > length(stcoords)
206  return
207  end
208  tc = stcoords(i);
209  end
210  end
211  end
212  end
213 
214  function new_vec = split(this, iind, split_point)
215  % function new_vec = split(this, iind, split_point)
216  % default implementation for the splitting of one node''s geometry
217  %
218  % This split works as follows:
219  % - If 'split_size' is equal to '2^(dimension)', the geometry is split
220  % in each of the 'dimension' directions.
221  % - If 'split_size' is equal to '2', the geometry is split in the
222  % direction where the geometry has its widest extent.
223  %
224  % Parameters:
225  % iind: the node''s ID
226  % split_point: an optional matrix with coordinate vectors specifying
227  % where the geometry shall be split. (default=barycenter)
228  %
229  % Return values:
230  % new_vec: vector with coordinates of new geometries.
231 
232  new_vec = zeros(this.dimension, this.split_size * this.ncoords);
233  coords = this.intervals(iind*this.skip+2:(iind+1)*this.skip);
234  coords = reshape(coords, this.dimension, this.ncoords)';
235  assert(this.ncoords == 2);
236  dim = this.dimension;
237  if(this.split_size == 2^this.dimension)
238  for i = 1:dim
239  c1 = repmat([coords(1,i),split_point(i)], 1, 2^(i-1));
240  c2 = repmat([split_point(i),coords(2,i)], 1, 2^(i-1));
241  new_vec(i,:) = repmat([c1, c2], 1, this.split_size/2^(i));
242  end
243  elseif this.split_size == 2
244  direction = find(split_point < 0, 1);
245  if isempty(direction)
246  [m, direction]=max(coords(2,:)-coords(1,:));
247  spc = split_point(direction);
248  else
249  spc = -split_point(direction);
250  end
251  c1 = coords(1,:)';
252  c2 = c1;
253  c2(direction) = spc;
254  c4 = coords(2,:)';
255  c3 = c4;
256  c3(direction) = spc;
257  new_vec = [c1, c2, c3, c4];
258  else
259  error('unsupport split_size');
260  end
261  new_vec = [zeros(1,this.split_size);...
262  reshape(new_vec, 2*dim, this.split_size)];
263  new_vec = reshape(new_vec, 1, this.skip*this.split_size);
264  end
265 
266  function [changed_ids, new_ids] = refine(this, iinds, split_points)
267  % function [changed_ids, new_ids] = refine(this, iinds, split_points)
268  % refines nodes given by indices at given points
269  %
270  % Parameters:
271  % iinds: a vector of node IDs for elements that shall be refined.
272  % split_points: a matrix with row vectors of points are given to the
273  % split() method as second argument. This argument is
274  % optional. (default = barycenter)
275  %
276  % Return values:
277  % changed_ids: vector of the IDs of the nodes that were changed during
278  % the refinement.
279  % new_ids: vector of the IDs of the nodes that were newly created
280  % during the refinement.
281  %
282  assert(nargin == 2 || length(iinds) == size(split_points,1));
283 
284  sskip = this.skip * this.split_size;
286  new_vector = zeros(1, sskip * length(iinds));
287  changed_ids = zeros(1, length(iinds));
288  new_ids = zeros(1, length(iinds) * this.split_size);
289 
290  new_levels = zeros(1, length(iinds) * this.split_size);
291  new_l_ids = zeros(1, length(iinds) * (this.split_size-1));
292 
293 
294  offset = length(this.intervals)/this.skip-1; % node ID of first new node
295  for i = 1:length(iinds)
296 
297  j=iinds(i); % node ID of refined element
298  this.intervals(j*this.skip+1) = offset + (i-1)*this.split_size+1;
299 
300  if nargin >= 3
301  split_point = split_points(i,:);
302  else
303  split_point = this.barycenter(j);
304  end
305  new_vector((i-1)*sskip+1:i*sskip) = this.split(j, split_point);
306 
307  changed_ids(i) = j;
308 
309  sprange = ((i-1)*this.split_size+1:i*this.split_size);
310  new_ids(sprange) = offset + sprange;
311  new_levels(sprange) = this.refine_levels(j+1) + 1;
312 
313  % update leaf_ids... do i really need it?
314  this.leaf_ids(find(this.leaf_ids==j, 1)) = offset + sprange(1);
315  new_l_ids((i-1)*(this.split_size-1)+1:i*(this.split_size-1)) = offset + sprange(2:end);
316  end
317 
318  this.intervals = [this.intervals, new_vector];
319  this.leaf_ids = [this.leaf_ids, new_l_ids];
320  this.refine_levels = [this.refine_levels, new_levels];
321 
322  % update sorted leafindexset
323  this.leaf_enum = zeros(1, length(this.leaf_ids));
324  j = -1;
325  for i = 1:length(this.leaf_enum)
326  % traverse to next element
327  while length(j) > 1 ...
328  && (j(end)+1 == length(this.intervals)/this.skip ...
329  || (~this.is_inside(this.barycenter(j(end)+1), j(end-1))))
330  j = j(1:end-1);
331  end
332  j(end) = j(end) + 1;
333 
334  fchild = this.intervals(j(end)*this.skip+1);
335  while(fchild > 0)
336  j = [j, fchild];
337  fchild = this.intervals(j(end)*this.skip+1);
338  end
339 
340  this.leaf_enum(i) = j(end);
341  end
342  end
343 
344  function bary = barycenter(this, iind)
345  % function bary = barycenter(this, iind)
346  % computes the barycenter of a node''s geometry
347  %
348  % Parameters:
349  % iind: node ID of the node whose gemoetry''s barycenter shall be
350  % computed
351  %
352  % Return values:
353  % bary: coordinate vector of the barycenter.
354  coords = this.intervals(iind*this.skip+2:(iind+1)*this.skip);
355  bary = 1/this.ncoords * sum(reshape(coords, this.dimension, this.ncoords));
356  bary = bary';
357  end
358 
359  function yes = intersect(this, coords, iind)
360  % function yes = intersect(this, coords, iind)
361  % checks whether a geometric entity intersects with another one given by
362  % a node''s ID
363  %
364  % Parameters:
365  % coords: '2 x dimension' matrix describing the geometric entity
366  % iind: ID of other node entity
367  %
368  % Return values:
369  % yes: boolean flag indicating whether the intersection exists
370  %
371 
372  if this.ncoords ~= 2
373  error('XPartMap.intersect() is only defined for ncoords = 2');
374  end
375 
376  extreme_coord = this.intervals(this.skip-this.dimension+1:this.skip);
377 
378  yes = false;
379  % check whether one coordinate is inside the node''s geometry
380  for i=1:size(coords,1)
381  if is_inside(this, coords(1,:), iind)
382  yes = true;
383  break;
384  end
385  end
386  % ... if not, they either do not intersect, or all node''s coordinate
387  % lie inside the geometry defined by 'coords', so we only check for
388  % the first one.
389  if ~yes
390  node_coords = get_coords(this, iind);
391  coord = node_coords(1, :);
392  if coords(2,:) == extreme_coord
393  if all(coord >= coords(1,:) & coord <= coords(2,:));
394  yes = true;
395  end
396  else
397  if all(coord >= coords(1,:) & coord < coords(2,:));
398  yes = true;
399  end
400  end
401  end
402  end
403 
404  function ok = is_inside(this, coord, iind)
405  % function ok = is_inside(this, coord, iind)
406  % checks whether a given coord lies in a node''s geometry
407  %
408  % Parameters:
409  % coord: the coordinate vector of length 'dimension'.
410  % iind: ID of node entity to check with
411  %
412  % Return values:
413  % ok: boolean flag
414 
415  if this.ncoords ~= 2
416  error('XPartMap.is_inside() is only defined for ncoords = 2');
417  end
418 
419  extreme_coord = this.intervals(this.skip-this.dimension+1:this.skip);
420  coords = get_coords(this, iind);
421  if coords(2, :) == extreme_coord
422  ok = all(coord >= coords(1,:) & coord <= coords(2,:));
423  else
424  ok = all(coord >= coords(1,:) & coord < coords(2,:));
425  end
426  end
427 
428  function [sxm, old_leaf_enum] = sub_xpart_map(this, coords)
429  % function sxm = sub_xpart_map(this, coords)
430  % creates new object holding a subset of the current XPartMap lying in
431  % the region specified by 'coords'
432  %
433  % Parameters:
434  % coords: a '2xdimension' matrix describing the region in which all the
435  % nodes of the resulting XPartMap shall lie.
436  %
437  % Return values:
438  % sxm: an object of type XPartMap holding a subset of the current map.
439 
440  root = get_coords(this, 0);
441  sxm = XPartMap(this.dimension, this.ncoords, root(:)');
442  sxm.intervals = this.intervals;
443 
444  new_ids = zeros(1:this.noof_ids);
445  skipped_ids = zeros(1:this.noof_ids);
446  skip = 0;
447  leaf_skip = 0;
448 
449  for i = 0:this.noof_ids-1
450  if intersect(this, coords, i)
451  ce = XPartMap.crop_entity(coords, get_coords(this, i));
452  sxm.intervals(i*this.skip+2:(i+1)*this.skip) = ce(:);
453  new_ids(i+1) = i - skip;
454  else
455  if sxm.intervals(i*this.skip+1) == 0
456  leaf_skip = leaf_skip + 1;
457  end
458  sxm.intervals(i*this.skip+1) = -1;
459  new_ids(i+1) = i - skip;
460  skipped_ids(i+1) = 1;
461  skip = skip + 1;
462  end
463  end
464 
465  new_intervals = zeros(1, this.skip * (this.noof_ids - skip));
466  new_refine_levels = zeros(1, this.noof_ids - skip);
467 
468  for i = 0:this.noof_ids-1
469  if sxm.intervals(i*this.skip+1) >= 0
470  ni = new_ids(i+1);
471  child_id = sxm.intervals(i*this.skip+1);
472  new_intervals(ni*this.skip+1) = new_ids(child_id+1);
473  new_intervals(ni*this.skip+2:(ni+1)*this.skip) ...
474  = sxm.intervals(i*this.skip+2:(i+1)*this.skip);
475  new_refine_levels(ni+1) = this.refine_levels(i+1);
476  end
477  end
478 
479  sxm.intervals = new_intervals;
480  sxm.refine_levels = new_refine_levels;
481  sxm.leaf_ids = new_ids(this.leaf_ids (~skipped_ids(this.leaf_ids+1))+1);
482  sxm.leaf_enum = new_ids(this.leaf_enum(~skipped_ids(this.leaf_enum+1))+1);
483 
484  % new_ids is monotone increasing, so...
485  old_leaf_enum = 1:length(this.leaf_enum);
486  old_leaf_enum = old_leaf_enum(~skipped_ids(this.leaf_enum+1));
487  end
488 
489  function coords = get_coords(this, iind)
490  % function coords = get_coords(this, iind)
491  % helper function rearranging the vector of node coordinates into a matrix
492  % with row vectors as coordinates.
493  %
494  % Parameters:
495  % iind: the node''s ID
496  %
497  % Return values:
498  % coords: the rearranged coordinates. ('ncoords x dimension' matrix)
499  coords = this.intervals(iind*this.skip+2:(iind+1)*this.skip);
500  coords = reshape(coords, this.dimension, this.ncoords)';
501  end
502  end
503 
504  methods (Static=true)
505 
506  function ce = crop_entity(coords1, coords2)
507  % function ce = crop_entity(coords1, coords2)
508  % helper function computing the actual intersection of two geometries
509  % given by their coordinates.
510  %
511  % Parameters:
512  % coords1: first '2 x dimension' matrix of cube extreme point
513  % coordinates.
514  % coords2: second '2 x dimension' matrix of cube extreme point
515  % coordinates.
516  %
517  % Return values:
518  % ce: intersection of the two geometries given by the two arguments.
519 
520  assert(all(size(coords1) == size(coords2)));
521  if size(coords1, 1) ~= 2
522  error('XPartMap.crop_entity() is only defined for ncoords = 2');
523  end
524  ce = zeros(size(coords1));
525  ce(1,:) = max(coords1(1,:), coords2(1,:));
526  ce(2,:) = min(coords1(2,:), coords2(2,:));
527 
528  if any(ce(1,:) >= ce(2,:))
529  ce = [];
530  end
531  end
532 
533  end
534 end
dimension
dimension of geometries
Definition: XPartMap.m:61
a geometric tree where every nodes has geometry information attached.
Definition: XPartMap.m:17
Definition: leaf.m:17