We describe the coarsening algorithm implemented in coarsen.m
for two dimensional bisection grids. We assume when bisecting a triangle, the left child is stored in a prori of the right child. More details can be found in the paper,
Every bisection grids can be obtained from an initial grid with a sequence of compatible bisections. In two dimensions, there are two types of compatible bisections.
node = [1,0; 0,1; -1,0; 0,-1];
elem = [2 3 1; 4 1 3];
figure(1);
subplot(1,2,1); showmesh(node,elem);
[node,elem] = bisect(node,elem,[1 2]);
subplot(1,2,2); showmesh(node,elem);
findnode(node,5,'noindex');
The nodes added by compatible bisections are called good-to-carsen nodes or simpley good nodes.
elem = [2 3 1];
subplot(1,2,1); showmesh(node,elem);
[node,elem] = bisect(node,elem,1);
subplot(1,2,2); showmesh(node,elem);
findnode(node,5,'noindex');
Our coarsening algorithm consists of finding good nodes and removing them. The coarsening procedure is illustrated by the following figures.
load Lshapemesh
figure
for k = 1:3
subplot(1,3,k); showmesh(node,elem);
[node,elem] = coarsen(node,elem,'all');
end
We use the following characterization to find a good node p
load Lshapemesh
N = size(node,1);
NT = size(elem,1);
valence = accumarray(elem(:),ones(3*NT,1),[N 1]);
valenceNew = accumarray(elem(:,1),ones(NT,1), [N 1]); % for newest vertex only
allNodes = (1:N)';
intGoodNode = allNodes((valence==valenceNew) & (valence==4));
bdGoodNode = allNodes((valence==valenceNew) & (valence==2));
figure(2); showmesh(node,elem);
findnode(node,intGoodNode,'noindex');
findnode(node,bdGoodNode,'noindex');
The star of good nodes are found using the incides matrix between triangles and vertices.
t2v = sparse([1:NT,1:NT,1:NT], elem(1:NT,:), 1, NT, N);
[ii,jj] = find(t2v(:,intGoodNode));
nodeStar = reshape(ii,4,length(intGoodNode));
For interior good nodes, there are 3 configurations of 4 triangles in the
star of a good node showing the following figures. Other combinations are
ruled out by our ordering assumption: the left child is stored in a prori
of the right child. Therefe nodeStar(:,1)
is always in L type and
nodeStar(:,4)
is in R type. Here the oritentation is respect the
direction from the new added vertex to the opposite vertex in each
triangle.
figure(1);
node1 = [1,0; 0,1; -1,0; 0,-1; 0,0];
elem1 = [5 3 2; 5 1 2; 5 4 1; 5 3 4];
subplot(1,3,1); showmesh(node1,elem1); findelem(node1,elem1);
findnode(node1,5,'noindex');
elem1 = [5 3 2; 5 4 1; 5 1 2; 5 3 4];
subplot(1,3,2); showmesh(node1,elem1); findelem(node1,elem1);
findnode(node1,5,'noindex');
elem1 = [5 3 2; 5 4 1; 5 3 4; 5 1 2];
subplot(1,3,3); showmesh(node1,elem1); findelem(node1,elem1);
findnode(node1,5,'noindex')
% We switch the local indices of nodal star to the first case.
ix = (elem(nodeStar(1,:),3) == elem(nodeStar(4,:),2));
iy = (elem(nodeStar(2,:),2) == elem(nodeStar(3,:),3));
% nodeStar(:,ix & iy) = nodeStar([1 2 3 4],ix & iy);
nodeStar(:,ix & ~iy) = nodeStar([1 3 2 4],ix & ~iy);
nodeStar(:,~ix & iy) = nodeStar([1 4 2 3],~ix & iy);
% nodeStar(:,~ix & ~iy)= nodeStar([1 4 3 2],~ix & ~iy);
%%
% We label the local indices of nodes in the star in the following figure.
figure(2)
node1 = [0,0; -1,0; 1,0; 0,1; 0,-1];
elem1 = [1 4 2; 1 3 4; 1 2 5; 1 5 3];
showmesh(node1,elem1); findelem(node1,elem1); findnode(node1);
Then t1 and t2 are merged to a big triangle and stored in t1 and similarly t3 and t4 are merged together and stored in t3. Elements t2 and t4 will be removed later on.
t1 = nodeStar(1,:);
t2 = nodeStar(2,:);
t3 = nodeStar(3,:);
t4 = nodeStar(4,:);
% p1 = elem(t1,1);
p2 = elem(t1,3);
p3 = elem(t2,2);
p4 = elem(t1,2);
p5 = elem(t3,2);
elem(t1,:) = [p4 p2 p3];
elem(t2,1) = 0;
elem(t3,:) = [p5 p3 p2];
elem(t4,1) = 0;
For boundary good nodes, there is only one case.
figure(2)
node1 = [0,0; -1,0; 1,0; 0,1];
elem1 = [1 4 2; 1 3 4];
showmesh(node1,elem1); findelem(node1,elem1); findnode(node1);
[ii,jj] = find(t2v(:,bdGoodNode));
nodeStar = reshape(ii,2,length(bdGoodNode));
t1 = nodeStar(1,:);
t2 = nodeStar(2,:);
p1 = elem(t1,1);
p2 = elem(t1,3);
p3 = elem(t2,2);
p4 = elem(t1,2);
elem(t1,:) = [p4 p2 p3];
elem(t2,1) = 0;
The empty rows in node and elem matrices should be relased for the efficent usage of meomory.
elem((elem(:,1) == 0),:) = [];
isGoodNode = false(N,1);
isGoodNode(intGoodNode) = true;
isGoodNode(bdGoodNode) = true;
node(isGoodNode,:) = [];
Note that the indices of nodes in the coarsened mesh is changed while
elem still use the indices in the fine mesh. We build the indices map
between the fine grid to the coarse grid. For example, indexMap(10) = 6
means the 10-th node in the fine grid is now the 6-th node in the coarse
one.
indexMap = zeros(N,1);
indexMap(~isGoodNode)= 1:size(node,1);
elem = indexMap(elem);
In the application to adaptive finite element methods, we will remove good-to-coraen nodes whose star are marked for coarsening.
In the application to multigrid methods, we will record all neighboring
nodes in the nodal star of good nodes; see the functoin: uniformcoarsen
.
The 3-D coarsening algorithm is slightly complicated than the 2-D case.
bdFlag
stores information on boundary conditions. It is updated asbdFlag(t1,:) = [bdFlag(t1,2) bdFlag(t2,1) bdFlag(t1,1)];
bdFlag(t3,:) = [bdFlag(t3,2) bdFlag(t4,1) bdFlag(t3,1)];
bdFlag(t5,:) = [bdFlag(t5,2) bdFlag(t6,1) bdFlag(t5,1)];
bdFlag((elem(:,1) == 0),:) = [];
tree(:,1:3)
stores the binary tree of the coarsening. tree(:,1)
is the index of parent element in the coarsened mesh and tree(:,2:3)
are two children indices in the original mesh.tree = zeros(NTdead,3,'uint32');
tree(1:NTdead,1) = [t1'; t3'; t5(idx)'];
tree(1:NTdead,2) = [t1'; t3'; t5(idx)'];
tree(1:NTdead,3) = [t2'; t4'; t6(idx)'];