2 % A triangular conforming grid in two dimensions.
4 % General geometrical information is stored including neighbour
5 % information. Therefore also boundary neighbour relations can be
6 % specified. Boundary edges can be marked by
"painting rectangles": the
7 % edges with midpoints within such a rectangle are marked accordingly.
8 % By
this boundary edges can be marked
for later special use. Much
9 % (partially redundant) information is stored in the grid, which might be
10 % useful in simulations.
13 % number of interior edges
16 % number of boundary edges
24 % constructor of a triangular conform grid in 2 dimensions with following
27 % -
triagrid() : construction of a default triagrid (2d unit square,
28 % loaded from file, -1 as outer neighbour indices)
29 % - @link triagrid() triagrid(tgrid) @endlink : copy-constructor
30 % - @link triagrid() triagrid(params) @endlink : in this 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
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.
53 % - @link triagrid() triagrid(p,t,params) @endlink : generate triagrid from
54 % triangle-data with certain options.
55 % - p is assumed to be a 2 x npoints matrix with coordinates
56 % - t is assumed to be a XX x ntriangles matrix, but only first three
57 % rows are used == vertex indices, in clockwise order as default, all
58 % nondefined edges are set to -1, then the following "rectangles" are
59 % set as specified in params
61 % Using this class, grids from PDETOOLS can be used:
64 % %%% => generate your grid and export 'p' and 't' to MATLAB workspace
66 % save('mygrid','p','t')
67 % grid = triagrid(struct('grid_initfile','mygrid'));
70 % perhaps later: constructor by duneDGF-file?
71 % perhaps later: contructor-flag: full vs non-full
72 % => only compute redundant information if required.
75 % 'varargin' : variable number of input arguments, see above for description
76 % of possible configurations.
79 % 'grid' : generated triagrid
81 % generated fields of grid:
82 % nelements: number of elements
83 % nvertices: number of vertices
85 % A : vector of element area
86 % Ainv : vector of inverted element area
87 % X : vector of vertex x-coordinates
88 % Y : vector of vertex y-coordinates
89 % VI : matrix of vertex indices: 'VI(i,j)' is the global index of j-th
91 % CX : vector of centroid x-values
92 % CY : vector of centroid y-values
93 % NBI : 'NBI(i,j) = ' element index of j-th neighbour of element i
94 % boundary faces are set to -1 or negative values are requested by
95 % 'params.boundary_type'
96 % INB : 'INB(i,j) = ' local edge number in 'NBI(i,j)' leading from
97 % element 'NBI(i,j)' to element 'i', i.e. satisfying
98 % 'NBI(NBI(i,j), INB(i,j)) = i'
99 % EL : 'EL(i,j) = ' length of edge from element 'i' to neighbour 'j'
100 % DC : 'DC(i,j) = ' distance from centroid of element i to NB j
101 % for boundary elements, this is the distance to the reflected
102 % element (for use in boundary treatment)
103 % NX : 'NX(i,j) = ' x-coordinate of unit outer normal of edge from el
105 % NY : 'NY(i,j) = ' y-coordinate of unit outer normal of edge from el
107 % ECX : 'ECX(i,j) = ' x-coordinate of midpoint of edge from el 'i' to NB 'j'
108 % ECY : 'ECY(i,j) = ' y-coordinate of midpoint of edge from el 'i' to NB 'j'
109 % SX : vector of x-coordinates of point `S_i` (for rect: identical to
111 % SY : vector of y-coordinate of point `S_j` (for rect: identical to
113 % ESX : 'ESX(i,j) = ' x-coordinate of point `S_{ij}` on edge el i to NB j
114 % ESY :
'ESY(i,j) = ' y-coordinate of point `S_{ij}` on edge el i to NB j
115 % DS :
'DS(i,j) = ' distance from `S_i` of element i to `S_j` of NB j
116 %
for boundary elements,
this is the distance to the reflected
117 % element (
for use in boundary treatment)
118 % hmin : minimal element-diameter
119 % alpha: geometry bound (simultaneously satisfying `\alpha \cdot h_i^d
120 % \leq A(T_i)`, `\alpha \cdot \mbox{diam}(T_i) \leq h_i^{d-1}` and
121 % `\alpha \cdot h_i \leq `dist(midpoint `i` to any neigbour) )
122 % JIT : Jacobian inverse transposed 'JIT(i,:,:)' is a 2x2-matrix of the Jac.
123 % Inv. Tr. on element 'i'
125 % \note for diffusion-discretization with FV-schemes, points `S_i` must exist,
126 % such that `S_i S_j` is perpendicular to edge 'i j', the intersection points
127 % are denoted `S_{ij}`
130 % Bernard Haasdonk 9.5.2007
132 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
134 if (nargin==1) & isa(varargin{1},
'triagrid')
137 fnames = fieldnames(varargin{1});
138 for i=1:length(fnames)
139 grid.(fnames{i}) = varargin{1}.(fnames{i});
144 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
145 %
default constructor: unit square
147 % p = [0.0 0.5 1.0 0.0 0.5 1.0 0.0 0.5 1.0; ...
148 % 0.0 0.0 0.0 0.5 0.5 0.5 1.0 1.0 1.0 ];
149 % t = [1 2 2 3 4 5 5 6; ...
150 % 2 5 3 6 5 8 6 9; ...
155 % mark boundary in rectangle from [-1,-1] to [+2,+2] with index -1
156 params.bnd_rect_corner1 = [-1,-1]';
157 params.bnd_rect_corner2 = [2 2]
';
158 params.bnd_rect_index = [-1];
160 params = varargin{1};
161 if isfield(params,'grid_initfile
')
162 tmp = load(params.grid_initfile);
166 nx = params.xnumintervals;
167 ny = params.ynumintervals;
168 nvertices =(nx+1)*(ny+1);
170 dx = (params.xrange(2)-params.xrange(1))/nx;
171 dy = (params.yrange(2)-params.yrange(1))/ny;
173 % set vertex coordinates
174 vind = (1:(nx+1)*(ny+1));
175 X = mod(vind-1,nx+1) * dx + params.xrange(1);
176 Y = floor((vind-1)/(nx+1))*dy + params.yrange(1);
179 % identify triangles:
185 all_ind = 1:(nx+1)*(ny+1);
186 % index vector of lower left corners, i.e. not last col/row
187 lc_ind = find((mod(all_ind,nx+1)~=0) & (all_ind < (nx+1)*ny));
188 t1 = [lc_ind;lc_ind+1;lc_ind+2+nx]; % lower triangles
189 t2 = [lc_ind;lc_ind+nx+2;lc_ind+nx+1]; % upper triangles
192 else % read p-e-t information from varargin
198 params = varargin{3};
201 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
202 % construct from p,t and params
208 % nx = params.xnumintervals;
209 % ny = params.ynumintervals;
211 grid.nelements = size(t,2);
212 grid.nvertices = size(p,2);
215 % set vertex coordinates
219 % set element vertex indices: numbering counterclockwise
220 % counterclockwise, edge j connects vertex j and j+1
224 %
get areas of grid-cells after VI and X,Y are available
226 % A(a b c) = 0.5 |Det ((b-a) (c-a)) |
227 grid.A = abs(area_triangle(...
228 [grid.X(grid.VI(:,1)), ...
229 grid.Y(grid.VI(:,1))], ...
230 [grid.X(grid.VI(:,2)), ...
231 grid.Y(grid.VI(:,2))], ...
232 [grid.X(grid.VI(:,3)), ...
233 grid.Y(grid.VI(:,3))]));
235 % a11 = grid.X(grid.VI(:,2))-grid.X(grid.VI(:,1));
236 % a21 = grid.Y(grid.VI(:,2))-grid.Y(grid.VI(:,1));
237 % a12 = grid.X(grid.VI(:,3))-grid.X(grid.VI(:,1));
238 % a22 = grid.Y(grid.VI(:,3))-grid.Y(grid.VI(:,1));
239 % grid.A = 0.5 * (a11.*a22 - a21.*a12);
240 grid.Ainv = grid.A.^(-1);
242 % midpoint coordinates of grid-cells
243 % equals cog of corners
245 % coordinates of vertices: XX(elnum, vnum)
246 XX = reshape(grid.X(grid.VI(:)),size(grid.VI));
247 YY = reshape(grid.Y(grid.VI(:)),size(grid.VI));
248 grid.CX = mean(XX,2);
249 grid.CY = mean(YY,2);
251 % edge-vectors: [DX(el,edge) DX (el,edge)]
252 DX = [XX(:,2:3), XX(:,1)] - XX;
253 DY = [YY(:,2:3), YY(:,1)] - YY;
255 % matrix with edge-lengths
256 grid.EL = sqrt(DX.^2 + DY.^2);
258 % matrix with (unit) normal components
259 grid.NX = grid.EL.^(-1) .* DY;
260 grid.NY = - grid.EL.^(-1) .* DX;
262 % matrix with edge-midpoint-coordinates
263 % this computation yields epsilon-differing edge-midpoints on
264 % neighbouring elements, perhaps solve later.
265 grid.ECX = XX + DX * 0.5;
266 grid.ECY = YY + DY * 0.5;
268 %%% determine indices of neighbouring elements, default: -1
269 NBI = -1 * ones(grid.nelements, grid.nneigh);
270 INB = zeros(grid.nelements, grid.nneigh);
271 % and now it gets complicated...
273 % setup edge list e(1,i), e(2,i) are the point indices of the i-th edge
274 % eid(i) is the element number, which has the edge
275 % led(i) : local edge number, which is the edge
277 ee = zeros(2,grid.nelements*3);
278 % li = [1 2 2 3 3 1 4 5 5 6 6 4 ... ]
279 li = 1:grid.nelements*3;
281 li = reshape(li,6,grid.nelements);
282 li = [li(2:end,:); li(1,:)];
283 li = reshape(li,2,grid.nelements*3);
285 ee = sort(ee); % now smallest vindex above
286 elid = repmat(1:grid.nelements,3,1);
288 led = repmat((1:3)',1,grid.nelements);
291 % generate sort key from both indices
292 key = ee(1,:) + grid.nvertices * ee(2,:);
293 [keysort, ind] = sort(key);
295 elid_sort = elid(ind);
298 % find indices of duplicates by simple forward difference
299 % => these are inner edges!
301 diff = keysort(1:end-1) - [keysort(2:end)];
302 inner_edges = find(diff==0);
304 % prepare entries in NBI: for all i in inner_edges
305 % NBI(elid_sort(i),led_sort(i)) = elid_sort(i+1);
306 % INB(elid_sort(i),led_sort(i)) = led_sort(i+1);
309 li = sub2ind(size(NBI),...
310 elid_sort(inner_edges),...
311 led_sort(inner_edges));
312 NBI(li) = elid_sort(inner_edges+1);
313 INB(li) = led_sort(inner_edges+1);
314 li = sub2ind(size(NBI),...
315 elid_sort(inner_edges+1),...
316 led_sort(inner_edges+1));
317 NBI(li) = elid_sort(inner_edges);
318 INB(li) = led_sort(inner_edges);
320 % % find all element-pairs, which have a common edge
321 % % i.e. neighs(i,1) and neighs(i,2) are two elements, which
322 % % have the same edge i, where this edge is the
323 % % eind(i,1)-th edge (1..3) in the first element
324 % % eind(i,2)-th edge (1..3) in the first element
326 % nedges = max(t(:));
327 % ehist = hist(t(:),1:nedges);
329 % error('edge occuring in more than 2 triangles!!');
331 % inner_edges = find(ehist==2);
332 % neighs = zeros(length(inner_edges),2);
333 % eind = zeros(length(inner_edges),2);
335 % % somehow map edge-numbers => triangle numbers
336 % % need two of them, as occasionally global and local edge numbers
338 % % edgemap1(i,locedgenum) = elementno if edge i is local edge in elementno
339 % % edgemap2(i,locedgenum) = elementno if edge i is local edge in elementno
341 % edgemap1 = nan * ones(nedges,grid.nneigh);%
343 % [locedgenum,triangnum] = ind2sub(size(t),1:length(t(:)));
344 % li = sub2ind(size(edgemap1),t(:),locedgenum(:));
345 % edgemap1(li) = triangnum; % here due to duplicates, only the last
348 % % remove these entries from t and repeat for 2nd edgemap
350 % [edgenum,locedgenum] = find(~isnan(edgemap1));
351 % li = sub2ind(size(edgemap1),edgenum(:),locedgenum(:));
353 % li2 = sub2ind(size(t2),locedgenum(:),edgemap1(li));
357 % edgemap2 = nan * ones(nedges,grid.nneigh);
358 % nonzeros = find(t2~=0);
360 % [locedgenum,triangnum] = ind2sub(size(t2),1:length(nonzeros));
361 % li = sub2ind(size(edgemap2),t2(nonzeros),locedgenum(:));
362 % edgemap2(li) = triangnum;
363 % % now together edgemap1 and edgemap 2 represent the edge information.
365 % % extract inner edges
366 % edgemap = [edgemap1(inner_edges,:), edgemap2(inner_edges,:)];
367 % % check, that exactly two entries per row are non-nan
368 % entries = sum(~isnan(edgemap),2);
369 % if (min(entries)~=2) | max(entries)~=2
370 % disp('error in neighbour detection, please check');
374 % [mi, imi] = min(edgemap,[],2);
375 % [ma, ima] = max(edgemap,[],2);
379 % % then fill NBI and INB in a vectorized way, i.e.
381 % % NBI(neighs(i,1), eind(i,1)) = neighs(i,2);
382 % % NBI(neighs(i,2), eind(i,2)) = neighs(i,1);
383 % % INB(neighs(i,1), eind(i,1)) = eind(i,2);
384 % % INB(neighs(i,2), eind(i,2)) = eind(i,1);
386 % l1 = sub2ind( size(NBI), neighs(:,1) , eind(:,1));
387 % l2 = sub2ind( size(NBI), neighs(:,2) , eind(:,2));
388 % NBI(l1) = neighs(:,2);
389 % NBI(l2) = neighs(:,1);
390 % INB(l1) = eind(:,2);
391 % INB(l2) = eind(:,1);
393 % search non-inner boundaries and their midpoints
396 SX = grid.ECX(li(:));
397 SY = grid.ECY(li(:));
399 % default-boundary = -1
400 bnd_ind = -1 * ones(length(li),1);
402 if isfield(params,'bnd_rect_index')
403 if isfield(params,'boundary_type')
404 error(['Do not specify both bnd_rect_index and boundary_type', ...
407 if (max(params.bnd_rect_index)>0)
408 error('boundary indices must be negative!');
410 if size(params.bnd_rect_corner1,1) == 1
411 params.bnd_rect_corner1 = params.bnd_rect_corner1';
413 if size(params.bnd_rect_corner2,1) == 1
414 params.bnd_rect_corner2 = params.bnd_rect_corner2';
416 for i = 1:length(params.bnd_rect_index)
417 indx = (SX > params.bnd_rect_corner1(1,i)) & ...
418 (SX < params.bnd_rect_corner2(1,i)) & ...
419 (SY > params.bnd_rect_corner1(2,i)) & ...
420 (SY < params.bnd_rect_corner2(2,i));
421 bnd_ind(indx) = params.bnd_rect_index(i);
424 if isfield(params,'boundary_type')
425 bnd_ind = params.boundary_type([SX(:),SY(:)],params);
429 NBI(li) = bnd_ind'; % set neighbours to boundary
434 % check grid consistency:
435 nonzero = find(NBI>0); % vector with vector-indices
436 [i,j] = ind2sub(size(NBI), nonzero); % vectors with matrix indices
437 NBIind = NBI(nonzero); % vector with global neighbour indices
438 INBind = INB(nonzero);
439 i2 = sub2ind(size(NBI),NBIind, INBind);
442 % plot_element_data(grid,grid.NBI,params);
443 disp('neighbour indices are not consistent!!');
447 % matrix with centroid-distances
448 grid.DC = nan * ones(size(grid.CX,1),grid.nneigh);
449 nonzero = find(grid.NBI>0);
450 [elind, nbind] = ind2sub(size(grid.NBI),nonzero);
452 CXX = repmat(grid.CX(:),1,grid.nneigh);
453 CYY = repmat(grid.CY(:),1,grid.nneigh);
456 CXN(nonzero) = grid.CX(grid.NBI(nonzero));
459 CYN(nonzero) = grid.CY(grid.NBI(nonzero));
461 DC = sqrt((CXX-CXN).^2 + (CYY - CYN).^2);
463 grid.DC(nonzero) = DC(nonzero);
465 % computation of boundary distances: twice the distance to the border
467 nondef = find(isnan(grid.DC));
469 [elind, nind] = ind2sub(size(grid.DC),nondef);
470 nind_plus_one = nind + 1;
471 i = find(nind_plus_one > grid.nneigh);
472 nind_plus_one(i) = 1;
473 nondef_plus_one = sub2ind(size(grid.DC),elind,nind_plus_one);
474 nondef_plus_one = nondef_plus_one(:)';
476 d = zeros(length(nondef),1);
478 q = [grid.CX(elind)'; grid.CY(elind)'];
479 p1 = [grid.X(grid.VI(nondef)), ...
480 grid.Y(grid.VI(nondef))];
481 p2 = [grid.X(grid.VI(nondef_plus_one)), ...
482 grid.Y(grid.VI(nondef_plus_one))];
484 grid.DC(nondef) = 2 * dist_point_line(q,p1,p2);
486 % make entries of ECX, ECY exactly identical for neighbouring elements!
487 % currently by construction a small eps deviation is possible.
488 %averaging over all pairs is required
489 nonzero = find(grid.NBI>0); % vector with vector-indices
490 [i,j] = ind2sub(size(grid.NBI), nonzero); % vectors with matrix indices
491 NBIind = NBI(nonzero); % vector with global neighbour indices
492 INBind = INB(nonzero); % vector with local edge indices
493 i2 = sub2ind(size(NBI),NBIind, INBind);
494 % determine maximum difference in edge-midpoints, but exclude
495 % symmetry boundaries by relative error < 0.0001
496 diffx = abs(grid.ECX(nonzero)-grid.ECX(i2));
497 diffy = abs(grid.ECY(nonzero)-grid.ECY(i2));
498 fi = find ( (diffx/(max(grid.X)-min(grid.X)) < 0.0001) & ...
499 (diffy/(max(grid.Y)-min(grid.Y)) < 0.0001) );
506 grid.ECX(nonzero(fi)) = 0.5*(grid.ECX(nonzero(fi))+ grid.ECX(i2(fi)));
507 grid.ECY(nonzero(fi)) = 0.5*(grid.ECY(nonzero(fi))+ grid.ECY(i2(fi)));
509 % for diffusion discretization: Assumption of points with
510 % orthogonal connections to edges. Distances and intersections
511 % determined here. For triangles, this can simply be the
512 % circumcenters and the corresponding intersections
514 q = [grid.X(grid.VI(:,1)), ...
515 grid.Y(grid.VI(:,1))];
516 p1 = [grid.X(grid.VI(:,2)), ...
517 grid.Y(grid.VI(:,2))];
518 p2 = [grid.X(grid.VI(:,3)), ...
519 grid.Y(grid.VI(:,3))];
520 S = circumcenter_triangle(q,p1,p2);
524 % compute DS for inner edges (identical as DC computation!)
525 grid.DS = nan * ones(size(grid.CX,1),grid.nneigh);
526 nonzero = find(grid.NBI>0);
527 [elind, nbind] = ind2sub(size(grid.NBI),nonzero);
529 SXX = repmat(grid.SX(:),1,grid.nneigh);
530 SYY = repmat(grid.SY(:),1,grid.nneigh);
533 SXN(nonzero) = grid.SX(grid.NBI(nonzero));
536 SYN(nonzero) = grid.SY(grid.NBI(nonzero));
538 DS = sqrt((SXX-SXN).^2 + (SYY - SYN).^2);
540 grid.DS(nonzero) = DS(nonzero);
542 % computation of boundary distances: twice the distance to the border
544 nondef = find(isnan(grid.DS));
546 [elind, nind] = ind2sub(size(grid.DS),nondef);
547 nind_plus_one = nind + 1;
548 i = find(nind_plus_one > grid.nneigh);
549 nind_plus_one(i) = 1;
550 nondef_plus_one = sub2ind(size(grid.DS),elind,nind_plus_one);
551 nondef_plus_one = nondef_plus_one(:)';
553 d = zeros(length(nondef),1);
555 q = [grid.SX(elind)'; grid.SY(elind)'];
556 p1 = [grid.X(grid.VI(nondef)), ...
557 grid.Y(grid.VI(nondef))];
558 p2 = [grid.X(grid.VI(nondef_plus_one)), ...
559 grid.Y(grid.VI(nondef_plus_one))];
563 grid.DS(nondef) = 2 * dist_point_line(q,p1,p2);
565 if ~isempty(find(isnan(grid.DS)))
566 error('nans produced in grid generation!');
569 % intersections of S_i are the centroids of edges
573 % compute geometry constants for CFL-condition
574 h = max(grid.EL,[],2); % elementwise maximum edgelength
576 % alpha1 * h_j^2 <= |T_j|
577 alpha1 = min(grid.A .* h.^(-2));
578 % alpha2 * |boundary T_j| <= h_j
580 alpha2 = min(h .* b.^(-1));
581 % alpha3 h_j <= d_jl, d_jl distance of circumcenters
582 alpha3 = min(min(grid.DS,[],2) .* h.^(-1));
583 grid.alpha = min([alpha1,alpha2,alpha3]);
586 % disp(['Caution: Grid-geometry constant alpha less/equal zero', ...
587 % ' problematic ' ...
588 % ' in automatic CFL-conditions!']);
591 % DF = [(p2-p1) (p3-p1)] a 2x2 matrix
592 DF = zeros(grid.nelements,2,2); % jacobian
593 DF(:,1,1) = XX(:,2)-XX(:,1);
594 DF(:,2,1) = YY(:,2)-YY(:,1);
595 DF(:,1,2) = XX(:,3)-XX(:,1);
596 DF(:,2,2) = YY(:,3)-YY(:,1);
597 DetDF = DF(:,1,1).*DF(:,2,2)- DF(:,2,1).*DF(:,1,2); % determinant
598 DetDFinv = DetDF.^(-1);
599 JIT = zeros(grid.nelements,2,2);
600 % inversetransposed of A = (a,b; c,d) is A^-1 = 1/det A (d,-c; -b,a);
601 JIT(:,1,1)=DF(:,2,2).*DetDFinv;
602 JIT(:,1,2)=-DF(:,2,1).*DetDFinv;
603 JIT(:,2,1)=-DF(:,1,2).*DetDFinv;
604 JIT(:,2,2)=DF(:,1,1).*DetDFinv;
607 if (size(JIT,1)>=5) &(max(max(abs(...
608 transpose(reshape(JIT(5,:,:),2,2))*...
609 reshape(DF(5,:,:),2,2)-eye(2)...
611 error('error in JIT computation!');
615 grid.JIT = JIT; % perhaps later also store area, etc.
617 grid.nedges_boundary = length(find(grid.NBI<0));
619 grid.nedges_interior = 0.5*(3*grid.nelements - grid.nedges_boundary);
620 if ~isequal(ceil(grid.nedges_interior), ...
621 grid.nedges_interior);
622 error('number of inner edges odd!!!');
632 gridp = gridpart(grid,eind);
634 lcoord = llocal2local(grid, faceinds, llcoord);
636 glob = local2global(grid, einds, loc, params);
638 p = plot(gird, params);
640 grid = set_nbi(grid, nbind, values);
642 function gcopy = copy(grid)
646 % gcopy:
object of type
triagrid. This is a deep copy of the current instance.
656 [C, G] = aff_trafo_glob2loc(x0, y0);
658 [C, G] = aff_trafo_loc2glob(x0, y0);
660 [C, G] = aff_trafo_orig2ref(x0, y0, varargin);
662 loc = global2local(grid, elementid, glob);
664 micro2macro = micro2macro_map(microgrid, macrogrid);