2 % a geometric tree where every nodes has geometry information attached.
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.
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
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.
23 % So every node has a size of
'1 + ncoords * dimension'.
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
31 % -# a
"sorted leaf element index" is a consecutive enumeration of
leaf
32 % indices enumerating the
leaf elements from left-to-right.
35 properties(Access=
public)
37 dimension = 1; % dimension of geometries
39 ncoords = 2; % number of coordinates defining a geometry
41 % actual tree with attached geometries
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.
48 % integer specifying how many children are created at refinement step
51 % an
'unsorted' enumeration of the
leaf indices
53 % If a
leaf has not been refined, it keeps its position in
this set.
56 % a
'sorted' enumeration of the
leaf indices.
58 % This set is constructed by a deep traversal of the tree, with
59 % left-to-right enumeration of leafs.
62 % A vector which stores the refinement level
for every node.
68 % the number of the
leaf nodes
71 % storage size of a single tree node pair
'(c, coords)'
74 % the maximum refinement level of the tree
77 noof_ids; % the number of node IDs in the tree
82 function skip =
get.skip(
this)
83 skip = this.dimension*this.ncoords+1;
86 function siz =
get.siz(
this)
87 siz = length(this.leaf_ids);
90 function nids =
get.noof_ids(
this)
91 nids = length(this.intervals)/this.skip;
94 function mlevel =
get.max_level(
this)
95 mlevel = max(this.refine_levels);
98 function copy = copy(
this)
99 %
function copy = copy(
this)
100 % deep copy
for an
XPartMap instance
103 % copy: the copied
object of type
XPartMap
107 fns = setdiff(properties(
this), {
'siz',
'skip',
'max_level',
'noof_ids'});
109 for i = 1:length(fns)
110 copy.(fns{i}) = this.(fns{i});
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.
119 % dimension: @copybrief #dimension
120 % ncoords: @copybrief #ncoords
121 % init_coords: initial coordinates
for the root node
''s geometry.
124 tpm.split_size = 2^dimension;
127 tpm.ncoords = ncoords;
128 if tpm.ncoords > 2 && nargin == 2
129 error(
'if ncoords > 2, init_coords need to be specified');
133 tpm.intervals = [0, init_coords];
135 icoords = zeros(1, tpm.ncoords * tpm.dimension);
136 icoords(tpm.dimension+1:tpm.dimension*2) = 1;
137 tpm.intervals = [0 icoords];
141 error(
'level1 splits are only supported for ncoords == 2');
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, :);
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];
161 function elems = coords2elem(
this, tcoords)
162 %
function elems = coords2elem(
this, tcoords)
163 % returns element IDs of
leaf nodes whose geometries host the given
167 % tcoords:
'n x dimension' matrix of
'n' coordinate (row) vectors for
168 % which the enclosing entity ID shall be determined.
171 % elems: a 'n' vector holding the element IDs.
173 elems = -1 * ones(1, size(tcoords,1));
174 if this.dimension > 1 || size(tcoords,1) * this.max_level < this.siz
176 for i = 1:size(tcoords,1)
180 if this.is_inside(tc, j)
181 if this.intervals(j*this.skip+1) == 0
184 j = this.intervals(j*this.skip+1);
188 if j >= this.noof_ids
198 [stcoords, map] = sort(tcoords);
201 for j = this.leaf_enum
202 while this.is_inside(tc, j)
205 if i > length(stcoords)
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
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.
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)
230 % new_vec: vector with coordinates of new geometries.
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)
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));
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);
249 spc = -split_point(direction);
257 new_vec = [c1, c2, c3, c4];
259 error('unsupport split_size');
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);
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
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)
277 % changed_ids: vector of the IDs of the nodes that were changed during
279 % new_ids: vector of the IDs of the nodes that were newly created
280 % during the refinement.
282 assert(nargin == 2 || length(iinds) == size(split_points,1));
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);
290 new_levels = zeros(1, length(iinds) * this.split_size);
291 new_l_ids = zeros(1, length(iinds) * (this.split_size-1));
294 offset = length(this.intervals)/this.skip-1; % node ID of first new node
295 for i = 1:length(iinds)
297 j=iinds(i); % node ID of refined element
298 this.intervals(j*this.skip+1) = offset + (i-1)*this.split_size+1;
301 split_point = split_points(i,:);
303 split_point = this.barycenter(j);
305 new_vector((i-1)*sskip+1:i*sskip) = this.split(j, split_point);
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;
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);
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];
322 % update sorted leafindexset
323 this.leaf_enum = zeros(1, length(this.leaf_ids));
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))))
334 fchild = this.intervals(j(end)*this.skip+1);
337 fchild = this.intervals(j(end)*this.skip+1);
340 this.leaf_enum(i) = j(end);
344 function bary = barycenter(this, iind)
345 % function bary = barycenter(this, iind)
346 % computes the barycenter of a node''s geometry
349 % iind: node ID of the node whose gemoetry''s barycenter shall be
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));
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
365 % coords: '2 x dimension' matrix describing the geometric entity
366 % iind: ID of other node entity
369 % yes:
boolean flag indicating whether the intersection exists
373 error('
XPartMap.intersect() is only defined for ncoords = 2');
376 extreme_coord = this.intervals(this.skip-this.dimension+1:this.skip);
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)
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
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,:));
397 if all(coord >= coords(1,:) & coord < coords(2,:));
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
409 % coord: the coordinate vector of length 'dimension'.
410 % iind: ID of node entity to check with
416 error('
XPartMap.is_inside() is only defined for ncoords = 2');
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,:));
424 ok = all(coord >= coords(1,:) & coord < coords(2,:));
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'
434 % coords: a '2xdimension' matrix describing the region in which all the
435 % nodes of the resulting
XPartMap shall lie.
438 % sxm: an
object of type
XPartMap holding a subset of the current map.
440 root = get_coords(this, 0);
441 sxm =
XPartMap(this.dimension, this.ncoords, root(:)');
442 sxm.intervals = this.intervals;
444 new_ids = zeros(1:this.noof_ids);
445 skipped_ids = zeros(1:this.noof_ids);
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;
455 if sxm.intervals(i*this.skip+1) == 0
456 leaf_skip = leaf_skip + 1;
458 sxm.intervals(i*this.skip+1) = -1;
459 new_ids(i+1) = i - skip;
460 skipped_ids(i+1) = 1;
465 new_intervals = zeros(1, this.skip * (this.noof_ids - skip));
466 new_refine_levels = zeros(1, this.noof_ids - skip);
468 for i = 0:this.noof_ids-1
469 if sxm.intervals(i*this.skip+1) >= 0
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);
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);
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));
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.
495 % iind: the node''s ID
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)';
504 methods (Static=true)
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.
512 % coords1: first '2 x dimension' matrix of cube extreme point
514 % coords2: second '2 x dimension' matrix of cube extreme point
518 % ce: intersection of the two geometries given by the two arguments.
520 assert(all(size(coords1) == size(coords2)));
521 if size(coords1, 1) ~= 2
522 error('
XPartMap.crop_entity() is only defined for ncoords = 2');
524 ce = zeros(size(coords1));
525 ce(1,:) = max(coords1(1,:), coords2(1,:));
526 ce(2,:) = min(coords1(2,:), coords2(2,:));
528 if any(ce(1,:) >= ce(2,:))
dimension
dimension of geometries
a geometric tree where every nodes has geometry information attached.