amgmcoarsen.m

五点差分型多重网格方法：各种插值算子的比较）

function AMGmCoarsen(level)

% Ryan McKenzie
% Department of Computational Sciences
% University of Kentucky

%Uses subroutines written by Gundolf Haase

amg_globals;

%if the edge and edge connectivity matrices are empty for this level
if isempty( A(level).edges ) %determine them...
    %derive node numbering for edge matrices from element connectivity
    [Ee,ea] = GetEdgeElements_Schur( A(level).elements , A(level).elconn , PAIRS);
    %Ee is edge matrix having dimensions [2,2,n_edges]
    %ea is edge connectivity matrix having dimensions [2,n_edges]
else %bring them in for use in this algorithm
    Ee = A(level).edges;
    ea = A(level).edconn;
end

% Determine neighborhood for nodes which are connected via an edge.
% neigh(i).nn is the array of neigbors for node i (incl. node i).
neigh = NodeNeighborsEdges(ea);

%----------------------- INITIAL COARSENING PHASE -----------------------%

%initialize coarse, working, and fine sets to be empty
C_set = [];
F_set = [];

nnodes = size(neigh); %get the number of nodes in the fine grid

%initialize the working set to all fine grid nodes
U_set=[];
for m=1:nnodes(2) %for each node in the fine grid
    U_set = union(U_set, m); %add this node to the working set
end
    
for m=1:nnodes(2) %for each node in the fine grid
    SC(m).set = strong(m, U_set, level); %get the set of strongly connected nodes to node m
    SC(m).phi = size( SC(m).set ); %get the cardinality of the strongly connected set of node m
end

n=0; %counts the number of nodes which have not been marked coarse or fine
loopcount = 0;
while n < nnodes(2) %while there are nodes that have not been marked
    loopcount = loopcount+1;
    maxPHI = 0; i=0; %initialize node with maximal strongly connected set
    usize = size( U_set ); %get the size of the working set
    
    %find the node i in the working set with maximal strongly connected set
    for m=1:usize(2) %for each node m in the working set
        if SC( U_set(m) ).phi > maxPHI %if the strogly connected set of m is the maximum encountered
            i=U_set(m); %mark the node having the max value
            maxPHI = SC( U_set(m) ).phi; %remember the max value
        end
    end
    
    C_set = union(C_set, i); %add the node with maximal strongly connected set to the coarse grid
    U_set = setdiff(U_set, i); %remove that node from the working set
    n = n+1; %count that one more node has been marked
    
    SiU = intersect( SC(i).set , U_set ); %construct the set of strongly connected nodes to i which are still in the working set
    ssize = size( SiU ); %get the size of that set
    for j=1:ssize(2) %for each node in that set
        F_set = union(F_set, SiU( j ) ); %tentatively mark that node as fine
        U_set = setdiff(U_set, SiU( j ) ); %remove that node from the working set
        n = n+1; %count that one more node has been marked
        
        SjU = intersect( SC(j).set , U_set ); %construct the set of strongly connected nodes to j which are still in the working set
        sssize = size( SjU ); %get the size of that set
        for k=1:sssize(2) %for each node in that set
            SC( SjU(k) ).phi = SC( SjU(k) ).phi + 1; %incriment the number of strong connection to that node
        end
    end 
end
    
%------------------------ FINAL COARSENING PHASE ------------------------%

%skipped in this version...

%------------------- GENERATE INTERPOLATION MOLECULES -------------------%

fnodes = size(F_set); %get the number of selected fine nodes
cnodes = size(C_set); %get the number of selected coarse nodes

%set up storage for global interpolation operator
% zero matrix with dimensions (all nodes on fine grid) X (coarse selected nodes)
W(level).weight = zeros( nnodes(2), cnodes(2) );

%generate interpolation molecules for each node selected as fine (row in molecule matrix)
for ni=1:fnodes(2) %for each marked fine node ni
    NNsize = size( neigh( F_set(ni) ).nn, 2 ); %get the size of the node neighborhood around the i-th fine point
    NEsize = size( neigh( F_set(ni) ).ne, 2 ); %get the size of the edge neighborhood around the i-th fine point

    %make sets of C and F nodes in the neighborhood
    fineneighbors = intersect( neigh( F_set(ni) ).nn, F_set );
    fncount=size(fineneighbors);
    coarseneighbors = intersect( neigh( F_set(ni) ).nn, C_set );
    cncount=size(coarseneighbors);
    
    %generate a local node numbering for the neighborhood of i
    %fine points should be sorted to the front of the array
    loc_glob(1:fncount(2))=fineneighbors;
    loc_glob(fncount(2)+1:NNsize)=coarseneighbors;
    %find the local numbering of node ni
    [is,i] = ismember( F_set(ni), loc_glob );
    
    %initialize molecule components for this fine node to the appropriate sizes
    M( ni ).ff = zeros(fncount(2), fncount(2));
    M( ni ).fc = zeros(fncount(2), cncount(2));
    %M( ni ).cf = zeros(cncount(2), fncount(2));
    %M( ni ).cc = zeros(cncount(2), cncount(2));

    for ecount=1:NEsize %for each edge neighboring node i
        %get the global index in the connectivity matrix for this edge
        e = neigh( F_set(ni) ).ne(ecount);
        
        i1 = ea(1,e); %find the first endpoint of this edge in global numbering
        i2 = ea(2,e); %find the second endpoint of this edge in global numbering
        
        i1C=0;i1F=0;i2C=0;i2F=0; %initialize flags
        %attempt to find i1's local numbering from the coarse node neighbors of i
        [i1C, i1_loc] = ismember(i1, coarseneighbors);
        if ~i1C %if i1 is not a local course point
            %find its local numbering from the fine node neighbors of i
            [i1F, i1_loc] = ismember(i1, fineneighbors);
        end
        
        %attempt to find i2's local numbering from the coarse node neighbors of i
        [i2C, i2_loc] = ismember(i2, coarseneighbors);
        if ~i2C %if i1 is not a local course point
            %find its local numbering from the fine node neighbors of i
            [i2F, i2_loc] = ismember(i2, fineneighbors);
        end
        
        %place the edge matrix components of this edge into the appropriate molecule component
        if [i1F i2F] %if both endpoints are finepoints
            %add the edge matrix contributions into the Mcc molecule component
            M( ni ).ff(i1_loc,i1_loc) = M( ni ).ff(i1_loc,i1_loc) + Ee(1,1,e);
            M( ni ).ff(i1_loc,i2_loc) = M( ni ).ff(i1_loc,i2_loc) + Ee(1,2,e);
            M( ni ).ff(i2_loc,i1_loc) = M( ni ).ff(i2_loc,i1_loc) + Ee(2,1,e);
            M( ni ).ff(i2_loc,i2_loc) = M( ni ).ff(i2_loc,i2_loc) + Ee(2,2,e);
        %elseif [i1C i2C] %if both endpoints are coarse points
        %    %add the edge matrix contributions into the Mcc molecule component
        %    M( ni ).cc(i1_loc, i1_loc) = M( ni ).cc(i1_loc, i1_loc) + Ee(1,1,e);
        %    M( ni ).cc(i1_loc, i2_loc) = M( ni ).cc(i1_loc, i2_loc) + Ee(1,2,e);
        %    M( ni ).cc(i2_loc, i1_loc) = M( ni ).cc(i2_loc, i1_loc) + Ee(2,1,e);
        %    M( ni ).cc(i2_loc, i2_loc) = M( ni ).cc(i2_loc, i2_loc) + Ee(2,2,e);
        elseif i1C %endpoint one is coarse and two is fine
            %add the edge matrix contributions into the molecule
            %M( ni ).cc(i1_loc, i1_loc) = M( ni ).cc(i1_loc, i1_loc) + Ee(1,1,e);
            %M( ni ).cf(i1_loc, i2_loc) = M( ni ).cf(i1_loc, i2_loc) + Ee(1,2,e);
            M( ni ).fc(i2_loc, i1_loc) = M( ni ).fc(i2_loc, i1_loc) + Ee(2,1,e);
            M( ni ).ff(i2_loc, i2_loc) = M( ni ).ff(i2_loc, i2_loc) + Ee(2,2,e);
        else %endpoint one is fine and two is coarse
            %add the edge matrix contributions into the molecule
            M( ni ).ff(i1_loc, i1_loc) = M( ni ).ff(i1_loc, i1_loc) + Ee(1,1,e);
            M( ni ).fc(i1_loc, i2_loc) = M( ni ).fc(i1_loc, i2_loc) + Ee(1,2,e);
            %M( ni ).cf(i2_loc, i1_loc) = M( ni ).cf(i2_loc, i1_loc) + Ee(2,1,e);
            %M( ni ).cc(i2_loc, i2_loc) = M( ni ).cc(i2_loc, i2_loc) + Ee(2,2,e);
        end
        
        %see if there are any contributions from coarse edges not connected directly to i
        
        %generate the set of fine node neigbors to i, excluding i
        FNI = setdiff( fineneighbors, i );
        FNIsize = size( FNI, 2 );%get the size of that set
        for j=1:FNIsize %for each non-i fine direct neighbor of i, called j
            [is,j_loc] = ismember( FNI( j ), fineneighbors ); %local numbering of j
            for k=1:cncount(2) %and for each coarse neighbor of i, called k
                k_loc=k; %local numbering of k
                %find the edges which connect the fine direct neighbor 
                %of i, called j, to the coarse direct neighbor of i, called k