rbmatlab  1.16.09
 All Classes Namespaces Files Functions Variables Modules Pages
rectgrid.m
1 classdef rectgrid < gridbase
2  % A cartesian rectangular grid in two dimensions with axis parallel elements.
3  %
4  % General geometrical information is stored including neighbour information.
5  % Therefore also boundary neighbour relations can be specified. The boundary
6  % type is set to be symmetric by default. Additionally, "rectangles" can be
7  % defined: the edges with midpoints within such a rectangle are marked
8  % accordingly. By this, boundary edges can be marked for later special use
9  % Much (partially redundant) information is stored in the grid, which might
10  % be useful in simulations.
11  %
12  % - Sorting of vertices is always counterclockwise SE,NE,NW,SW
13  % - A local edge `j \in \{1,\dots,4\}` is connecting points `j` and `(j+1)
14  % \mod 4`
15 
16  methods
17 
18  function grid=rectgrid(varargin)
19  %function rectgrid(varargin)
20  % constructor of a cartesian rectangular grid in 2 dimensions with
21  % axis parallel elements.
22  %
23  % The constructor has sthe following
24  % Synopsis:
25  %
26  % - rectgrid() : construction of a default rectgrid (2d unit square,
27  % 2x2 elements with -1 as outer neighbour indices)
28  % - rectgrid(rgrid) : copy-constructor
29  % - rectgrid(params) : generate rectgrid with certain options. In this
30  % case either
31  % -# the field grid_initfile is existent in params. Then the file is read
32  % and a pointlist p and a triangle. Procedure is then executing the
33  % following constructor
34  % -# a structured equidistant triangle grid is generated
35  % - 'params.xnumintervals' : number of elements along x directions
36  % - 'params.ynumintervals' : number of elements along y directions
37  % - 'params.xrange,yrange' : intervals covered along the axes
38  % .
39  % with the optional fields
40  % - 'params.bnd_rect_corner1' : coordinates of lower corner of to
41  % be marked boundaries
42  % - 'params.bnd_rect_corner2' : coordinates of upper corner of to
43  % be marked boundaries
44  % - 'params.bnd_rect_index': integer index to be set on the edges
45  % in the above defined rectangle. Should not be positive integer
46  % in the range of the number of elements. use negative indices for
47  % certain later discrimination.
48  % .
49  % For the last three optional boundary settings, also multiple
50  % rectangles can be straightforwardly defined by accepting matrix
51  % of columnwise corners1, corners2 and a vector of indices for the
52  % different rectangles.
53  %
54  % perhaps later: constructor by duneDGF-file?
55  % perhaps later: contructor-flag: full vs non-full
56  % => only compute redundant information if required.
57  %
58  % Parameters:
59  % 'varargin' : variable number of input arguments, see above for description
60  % of possible configurations.
61  %
62  % Return values:
63  % 'grid' : generated rectgrid
64  %
65  % generated fields of grid:
66  % nelements: number of elements
67  % nvertices: number of vertices
68  % nneigh: 4
69  % A : @copybrief gridbaseA
70  % Ainv : vector of inverted element area
71  % X : vector of vertex x-coordinates
72  % Y : vector of vertex y-coordinates
73  % VI : matrix of vertex indices: 'VI(i,j)' is the global index of j-th
74  % vertex of element i
75  % CX : vector of centroid x-values
76  % CY : vector of centroid y-values
77  % NBI : 'NBI(i,j) = ' element index of j-th neighbour of element i
78  % boundary faces are set to -1 or negative values are requested by
79  % 'params.boundary_type'
80  % INB : 'INB(i,j) = ' local edge number in 'NBI(i,j)' leading from
81  % element 'NBI(i,j)' to element 'i', i.e. satisfying
82  % 'NBI(NBI(i,j), INB(i,j)) = i'
83  % EL : 'EL(i,j) = ' length of edge from element 'i' to neighbour 'j'
84  % DC : 'DC(i,j) = ' distance from centroid of element i to NB j
85  % for boundary elements, this is the distance to the reflected
86  % element (for use in boundary treatment)
87  % NX : 'NX(i,j) = ' x-coordinate of unit outer normal of edge from el
88  % 'i' to NB 'j'
89  % NY : 'NY(i,j) = ' y-coordinate of unit outer normal of edge from el
90  % 'i' to NB 'j'
91  % ECX : 'ECX(i,j) = ' x-coordinate of midpoint of edge from el 'i' to NB 'j'
92  % ECY : 'ECY(i,j) = ' y-coordinate of midpoint of edge from el 'i' to NB 'j'
93  % SX : vector of x-coordinates of point `S_i` (for rect: identical to
94  % centroids)
95  % SY : vector of y-coordinate of point `S_j` (for rect: identical to
96  % centroids)
97  % ESX : 'ESX(i,j) = ' x-coordinate of point `S_{ij}` on edge el i to NB j
98  % ESY : 'ESY(i,j) = ' y-coordinate of point `S_{ij}` on edge el i to NB j
99  % DS : 'DS(i,j) = ' distance from `S_i` of element i to `S_j` of NB j
100  % for boundary elements, this is the distance to the reflected
101  % element (for use in boundary treatment)
102  % hmin : minimal element-diameter
103  % alpha: geometry bound (simultaneously satisfying `\alpha \cdot h_i^d
104  % \leq A(T_i)`, `\alpha \cdot \mbox{diam}(T_i) \leq h_i^{d-1}` and
105  % `\alpha \cdot h_i \leq `dist(midpoint `i` to any neigbour) )
106  %
107  % \note for diffusion-discretization with FV-schemes, points `S_i` must exist,
108  % such that `S_i S_j` is perpendicular to edge 'i j', the intersection points
109  % are denoted `S_{ij}`
110  %
111  % Bernard Haasdonk 9.5.2007
112 
113  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
114  % copy constructor
115  if (nargin==1) && ...
116  isa(varargin{1},'rectgrid')
117  grid = rectgrid;
118 
119  fnames = fieldnames(varargin{1});
120  for i=1:length(fnames)
121  grid.(fnames{i}) = varargin{1}.(fnames{i});
122  end
123  else
124  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
125  % default constructor: unit square
126  if (nargin==0)
127  params.xnumintervals = 2;
128  params.ynumintervals = 2;
129  params.xrange = [0,1];
130  params.yrange = [0,1];
131  % mark element in rectangle from [-1,-1] to [+2,+2] with index -1
132  params.bnd_rect_corner1 = [-1,-1]';
133  params.bnd_rect_corner2 = [2 2]';
134  params.bnd_rect_index = [-1];
135  params.verbose = 0;
136  else
137  params = varargin{1};
138  end;
139 
140  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
141  % construct from params
142 
143  % if ~isfield(params,'verbose')
144  % params.verbose = 0;
145  % end;
146 
147  nx = params.xnumintervals;
148  ny = params.ynumintervals;
149 
150  grid.nelements =nx*ny;
151  grid.nvertices =(nx+1)*(ny+1);
152 
153  % get areas of grid-cells
154  dx = (params.xrange(2)-params.xrange(1))/nx;
155  dy = (params.yrange(2)-params.yrange(1))/ny;
156  grid.A = dx*dy*ones(nx*ny,1);
157  grid.Ainv = grid.A.^(-1);
158 
159  % set vertex coordinates
160  vind = (1:(nx+1)*(ny+1))';
161  grid.X = mod(vind-1,nx+1) * dx + params.xrange(1);
162  grid.Y = floor((vind-1)/(nx+1))*dy + params.yrange(1);
163 
164  % set element vertex indices: numbering starting right lower corner
165  % counterclockwise, i.e. edge j connects vertex j and j+1
166  % nx+2--nx+3-- ...
167  % | |
168  % | 1 | 2 ...
169  % | |
170  % 1-----2---- ...
171  % => grid.VI = [ 2 nx+3 nx+2 1];
172 
173  el_ind = 1:length(grid.A);
174  col_ind = transpose(mod((el_ind-1),nx)+1);
175  row_ind = transpose(floor((el_ind-1)/nx)+1);
176  VI = zeros(length(grid.A),4);
177  VI(:,1) = (row_ind-1)*(nx+1)+1 +col_ind;
178  VI(:,2) = (row_ind )*(nx+1)+1 +col_ind;
179  VI(:,3) = (row_ind )*(nx+1) +col_ind;
180  VI(:,4) = (row_ind-1)*(nx+1) +col_ind;
181  grid.VI = VI;
182 
183  % midpoint coordinates of grid-cells
184  CX = (0:(nx-1))*(params.xrange(2)-params.xrange(1))/ ...
185  nx+dx/2+params.xrange(1);
186  CX = CX(:);
187  grid.CX = repmat(CX,ny,1);
188  CY = (0:(ny-1))*(params.yrange(2)-params.yrange(1))/ ...
189  ny+dy/2+ params.yrange(1);
190  CY = repmat(CY,nx,1);
191  grid.CY = CY(:);
192  %disp('stopping after COG computation');
193  %keyboard;
194 
195  % check consistency: grid-midpoints and vertices
196  xdiff1 = max(abs(grid.CX + dx/2 -grid.X(grid.VI(:,1))));
197  xdiff2 = max(abs(grid.CX + dx/2 -grid.X(grid.VI(:,2))));
198  xdiff3 = max(abs(grid.CX - dx/2 -grid.X(grid.VI(:,3))));
199  xdiff4 = max(abs(grid.CX - dx/2 -grid.X(grid.VI(:,4))));
200  ydiff1 = max(abs(grid.CY - dy/2 -grid.Y(grid.VI(:,1))));
201  ydiff2 = max(abs(grid.CY + dy/2 -grid.Y(grid.VI(:,2))));
202  ydiff3 = max(abs(grid.CY + dy/2 -grid.Y(grid.VI(:,3))));
203  ydiff4 = max(abs(grid.CY - dy/2 -grid.Y(grid.VI(:,4))));
204 
205  if isfield(params, 'verbose') && params.verbose>=10
206  disp([xdiff1,xdiff2,xdiff3,xdiff4, ...
207  ydiff1,ydiff2,ydiff3,ydiff4]);
208 
209  if max([xdiff1,xdiff2,xdiff3,xdiff4, ...
210  ydiff1,ydiff2,ydiff3,ydiff4] > eps)
211  error('vertex coordinate and element midpoint consistency!!');
212  end;
213  end;
214 
215  % matrix with edge-lengths
216  grid.EL = repmat([dy, dx, dy, dx ],size(grid.CX,1),1);
217 
218  % matrix with midpoint-distances
219  grid.DC = repmat([dx, dy, dx, dy],size(grid.CX,1),1);
220 
221  % matrix with (unit) normal components
222  grid.NX = repmat([1,0,-1,0],size(grid.CX,1),1);
223  grid.NY = repmat([0,1,0,-1],size(grid.CX,1),1);
224 
225  % matrix with edge-midpoint-coordinates
226  % this computation yields epsilon-differing edge-midpoints on
227  % neighbouring elements. This is adjusted at end of this routine
228  grid.ECX = [grid.CX + dx/2,grid.CX,grid.CX-dx/2,grid.CX];
229  grid.ECY = [grid.CY, grid.CY + dy/2,grid.CY,grid.CY-dy/2];
230 
231  %%% determine indices of neighbouring elements,
232  NBI = repmat((1:nx*ny)',1,4);
233  % first column: +x direction: cyclical shift +1 of indices
234  NBI(:,1) = NBI(:,1)+1;
235  % second column: +y direction
236  NBI(:,2) = NBI(:,2)+nx;
237  % third column: -x direction
238  NBI(:,3) = NBI(:,3)-1;
239  % fourth column: -y direction
240  NBI(:,4) = NBI(:,4)-nx;
241 
242  % correct boundary elements neighbour-indices:
243  bnd_i1 = (1:ny)*nx; % indices of right-column elements
244  bnd_i2 = (1:nx)+ nx*(ny-1); % indices of upper el
245  bnd_i3 = (1:ny)*nx-nx+1; % indices of left-column
246  bnd_i4 = 1:nx; % indices of lower row elements
247  SX = [grid.ECX(bnd_i1,1); grid.ECX(bnd_i2,2); grid.ECX(bnd_i3,3); grid.ECX(bnd_i4, 4)];
248  SY = [grid.ECY(bnd_i1,1); grid.ECY(bnd_i2,2); grid.ECY(bnd_i3,3); grid.ECY(bnd_i4, 4)];
249 
250  % formerly default: Dirichlet :
251  % bnd_ind = -1 * ones(1,length(SX));
252  % now: set default to "symmetric", i.e. rectangle is a torus
253 
254  bnd_ind = [NBI(bnd_i1,1)'-nx, NBI(bnd_i2,2)'-nx*ny, ...
255  NBI(bnd_i3,3)'+nx, NBI(bnd_i4,4)'+nx*ny];
256 
257  % disp('halt 1');
258  % keyboard;
259 
260  if ~isfield(params, 'bnd_rect_index')
261  params.bnd_rect_index = [];
262  end
263 
264  if ~isempty(params.bnd_rect_index)
265  % keyboard;
266  if (max(params.bnd_rect_index)>0)
267  error('boundary indices must be negative!');
268  end;
269  if size(params.bnd_rect_corner1,1) == 1
270  params.bnd_rect_corner1 = params.bnd_rect_corner1';
271  end;
272  if size(params.bnd_rect_corner2,1) == 1
273  params.bnd_rect_corner2 = params.bnd_rect_corner2';
274  end;
275  for i = 1:length(params.bnd_rect_index)
276  indx = (SX > params.bnd_rect_corner1(1,i)) & ...
277  (SX < params.bnd_rect_corner2(1,i)) & ...
278  (SY > params.bnd_rect_corner1(2,i)) & ...
279  (SY < params.bnd_rect_corner2(2,i));
280  bnd_ind(find(indx)) = params.bnd_rect_index(i);
281  end;
282  end;
283  % disp('halt 2');
284 % keyboard;
285 
286  iend1 = length(bnd_i1);
287  iend2 = iend1 + length(bnd_i2);
288  iend3 = iend2 + length(bnd_i3);
289  iend4 = iend3 + length(bnd_i4);
290  NBI(bnd_i1,1) = bnd_ind(1:iend1)'; % set right border-neigbours to boundary
291  NBI(bnd_i2,2) = bnd_ind((iend1+1):iend2)'; % set upper neighbours to boundary
292  NBI(bnd_i3,3) = bnd_ind((iend2+1):iend3)'; % set left neigbours to boundary
293  NBI(bnd_i4,4) = bnd_ind((iend3+1):iend4)'; % set lower neighbours to boundary
294 
295  grid.NBI = NBI;
296 
297  % INB: INB(i,j) = local edge number in NBI(i,j) leading from element
298  % NBI(i,j) to element i, i.e. satisfying
299  % NBI(NBI(i,j), INB(i,j)) = i
300  INB = repmat([3 4 1 2],size(grid.NBI,1),1);
301  grid.INB = INB;
302 
303  % check grid consistency:
304  nonzero = find(NBI>0); % vector with element indices
305  [i,j] = ind2sub(size(NBI), nonzero); % vectors with matrix indices
306  NBIind = NBI(nonzero); % vector with global neighbour indices
307  INBind = INB(nonzero);
308  i2 = sub2ind(size(NBI),NBIind, INBind);
309  i3 = NBI(i2);
310  if ~isequal(i3,i)
311 % plot_element_data(grid,grid.NBI,params);
312  disp('neighbour indices are not consistent!!');
313  keyboard;
314  end;
315 
316  grid.hmin = sqrt(dx^2+dy^2); % minimal diameter of elements
317  % CAUTION: this geometry bound is
318  % adapted to triangles, should be extended
319  % properly to rectangles
320  alpha1 = dx * dy/(dx^2+dy^2); % geometry bound
321  alpha2 = sqrt(dx^2+dy^2)/(2*dx + 2*dy);
322  alpha3 = dx / sqrt(dx^2+dy^2);
323  grid.alpha = min([alpha1,alpha2,alpha3]);
324 
325  % make entries of ECX, ECY exactly identical for neighbouring elements!
326  % currently by construction a small eps deviation is possible.
327 
328  %averaging over all pairs is required
329  nonzero = find(grid.NBI>0); % vector with vector-indices
330  [i,j] = ind2sub(size(grid.NBI), nonzero); % vectors with matrix indices
331  NBIind = NBI(nonzero); % vector with global neighbour indices
332  INBind = INB(nonzero); % vector with local edge indices
333  i2 = sub2ind(size(NBI),NBIind, INBind);
334  % determine maximum difference in edge-midpoints, but exclude
335  % symmetry boundaries by relative error < 0.0001
336  diffx = abs(grid.ECX(nonzero)-grid.ECX(i2));
337  diffy = abs(grid.ECY(nonzero)-grid.ECY(i2));
338  fi = find ( (diffx/(max(grid.X)-min(grid.X)) < 0.0001) & ...
339  (diffy/(max(grid.Y)-min(grid.Y)) < 0.0001) );
340 
341  %disp(max(diffx));
342  %disp(max(diffy));
343  % => 0 ! :-)
344  % keyboard;
345 
346  grid.ECX(nonzero(fi)) = 0.5*(grid.ECX(nonzero(fi))+ grid.ECX(i2(fi)));
347  grid.ECY(nonzero(fi)) = 0.5*(grid.ECY(nonzero(fi))+ grid.ECY(i2(fi)));
348 
349  % for diffusion discretization: Assumption of points with
350  % orthogonal connections to edges. Distances and intersections
351  % determined here. In cartesian case identical to centroids
352  grid.SX = grid.CX;
353  grid.SY = grid.CY;
354  grid.ESX = grid.ECX;
355  grid.ESY = grid.ECY;
356  grid.DS = grid.DC;
357  grid.nneigh= 4;
358 
359  end
360 
361  end
362 
363  % demo
364  demo(dummy);
365 
366  p=plot(grid, params);
367 
368  res = isequal(grid,grid2);
369 
370  grid=set_nbi(grid, nbind, values);
371 
372  function gcopy = copy(grid)
373  % @copybrief gridbasecopy
374  %
375  % @copydoc gridbasecopy
376 
377  gcopy = rectgrid(grid);
378  end
379 
380  glob = local2global(grid, einds, loc, params);
381 
382  end
383 end
384 
function r = verbose(level, message, messageId)
This function displays messages depending on a message-id and/or a level. Aditionally you can set/res...
Definition: verbose.m:17
A cartesian rectangular grid in two dimensions with axis parallel elements.
Definition: rectgrid.m:17
Base class for all grid classes.
Definition: gridbase.m:17
nelements
number of overall elements (leaf + nonleaf)
Definition: gridbase.m:30