2 % A cartesian rectangular grid in two dimensions with axis parallel elements.
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.
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)
20 % constructor of a cartesian rectangular grid in 2 dimensions with
21 % axis parallel elements.
23 % The constructor has sthe following
27 % 2x2 elements with -1 as outer neighbour indices)
28 % -
rectgrid(rgrid) : copy-constructor
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
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.
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.
54 % perhaps later: constructor by duneDGF-file?
55 % perhaps later: contructor-flag: full vs non-full
56 % => only compute redundant information if required.
59 %
'varargin' : variable number of input arguments, see above for description
60 % of possible configurations.
65 % generated fields of grid:
66 % nelements: number of elements
67 % nvertices: number of vertices
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
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
89 % NY :
'NY(i,j) = ' y-coordinate of unit outer normal of edge from el
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
95 % SY : vector of y-coordinate of point `S_j` (for rect: identical to
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) )
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}`
111 % Bernard Haasdonk 9.5.2007
113 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
115 if (nargin==1) && ...
116 isa(varargin{1},
'rectgrid')
119 fnames = fieldnames(varargin{1});
120 for i=1:length(fnames)
121 grid.(fnames{i}) = varargin{1}.(fnames{i});
124 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
125 %
default constructor: unit square
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];
137 params = varargin{1};
140 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
141 % construct from params
143 %
if ~isfield(params,
'verbose')
144 % params.verbose = 0;
147 nx = params.xnumintervals;
148 ny = params.ynumintervals;
151 grid.nvertices =(nx+1)*(ny+1);
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);
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);
164 % set element vertex indices: numbering starting right lower corner
165 % counterclockwise, i.e. edge j connects vertex j and j+1
171 % => grid.VI = [ 2 nx+3 nx+2 1];
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;
183 % midpoint coordinates of grid-cells
184 CX = (0:(nx-1))*(params.xrange(2)-params.xrange(1))/ ...
185 nx+dx/2+params.xrange(1);
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);
192 %disp('stopping after COG computation
');
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))));
205 if isfield(params, 'verbose') && params.verbose>=10
206 disp([xdiff1,xdiff2,xdiff3,xdiff4, ...
207 ydiff1,ydiff2,ydiff3,ydiff4]);
209 if max([xdiff1,xdiff2,xdiff3,xdiff4, ...
210 ydiff1,ydiff2,ydiff3,ydiff4] > eps)
211 error('vertex coordinate and element midpoint consistency!!
');
215 % matrix with edge-lengths
216 grid.EL = repmat([dy, dx, dy, dx ],size(grid.CX,1),1);
218 % matrix with midpoint-distances
219 grid.DC = repmat([dx, dy, dx, dy],size(grid.CX,1),1);
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);
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];
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;
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)];
250 % formerly
default: Dirichlet :
251 % bnd_ind = -1 * ones(1,length(SX));
252 % now: set
default to
"symmetric", i.e. rectangle is a torus
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];
260 if ~isfield(params,
'bnd_rect_index')
261 params.bnd_rect_index = [];
264 if ~isempty(params.bnd_rect_index)
266 if (max(params.bnd_rect_index)>0)
267 error('boundary indices must be negative!');
269 if size(params.bnd_rect_corner1,1) == 1
270 params.bnd_rect_corner1 = params.bnd_rect_corner1';
272 if size(params.bnd_rect_corner2,1) == 1
273 params.bnd_rect_corner2 = params.bnd_rect_corner2';
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);
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
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);
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);
311 % plot_element_data(grid,grid.NBI,params);
312 disp('neighbour indices are not consistent!!');
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]);
325 % make entries of ECX, ECY exactly identical for neighbouring elements!
326 % currently by construction a small eps deviation is possible.
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) );
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)));
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
366 p=plot(grid, params);
368 res = isequal(grid,grid2);
370 grid=set_nbi(grid, nbind, values);
372 function gcopy = copy(grid)
380 glob = local2global(grid, einds, loc, params);
function r = verbose(level, message, messageId)
This function displays messages depending on a message-id and/or a level. Aditionally you can set/res...
A cartesian rectangular grid in two dimensions with axis parallel elements.
Base class for all grid classes.
nelements
number of overall elements (leaf + nonleaf)