amgmcoarsen.m

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

                %NOTE: such an edge is not guaranteed to exist, there can be at most one
                IND = intersect( neigh( FNI(j) ).ne , neigh( coarseneighbors(k) ).ne );
                
                if ~isempty(IND) %if such an edge was found between j and k
                    ne = IND(1); %get the global indedx for the connectivity matrix of the edge
                    
                    %set up the indices for the edge matrices
                    if ea(1,ne) == loc_glob(k_loc) %if the first endpoint is our coarse point k
                    %NOTE: in this case the second endpoint must be our fine node j
                        ee1=1; ee2=2;
                    else %endpoint one is fine and two is coarse
                        ee1=2; ee2=1;
                    end
                    %add the edge matrix contributions into the molecule
                    %M( ni ).cc(k_loc, k_loc) = M( ni ).cc(k_loc, k_loc) + Ee(ee1,ee1,ne);
                    %M( ni ).cf(k_loc, j_loc) = M( ni ).cf(k_loc, j_loc) + Ee(ee1,ee2,ne);
                    M( ni ).fc(j_loc, k_loc) = M( ni ).fc(j_loc, k_loc) + Ee(ee2,ee1,ne);
                    M( ni ).ff(j_loc, j_loc) = M( ni ).ff(j_loc, j_loc) + Ee(ee2,ee2,ne);
                end
            end
        end 
    end
 
    %store the local to global numbering utilized in this ni-th molecule for later use
    %M( ni ).numbering = loc_glob;

    %generate the fncount by cncount matrix Pfc which will be used as the
    %interpolation weight matrix for this fine point
    M( ni ).weight = (-(inv( M( ni ).ff ))) * M( ni ).fc;
    
    %extract the i_th row of Pfc to go into the global interpolation operator
    for cc=1:cncount(2)
        globalcoarse = loc_glob(fncount(2)+cc);
        [is,ci] = ismember( globalcoarse, C_set );
        W(level).weight( F_set(ni), ci ) = M(ni).weight(i,cc);
    end
end

%-------------------------- COARSE LEVEL SETUP --------------------------%

%generate coarse edges
currentCedge = 0;

%get the number of coarsely selected points
cnodes = size(C_set);

%for each selected coarse point i
for icount=1:cnodes(2)
    i = C_set( icount ); %get the global point ID for i
    
    %for each selected coarse point j
    for jcount=icount+1:cnodes(2)
        j = C_set( jcount ); %get the global point ID for j
        
        %examine the possibility of making a coarse edge between i and j
        if i ~= j %do not consider cases where it is a self connection...
            
            CE(icount,jcount).cc = zeros(2,2);
            CE(icount,jcount).cf = zeros(2,2);
            CE(icount,jcount).fc = zeros(2,2);
            
            %if i and j are directly connected on the fine grid
            if ismember( j, neigh(i).nn )
                %create a coarse edge between i and j
                currentCedge = currentCedge + 1;
                %place the edge in the connectivity matrix
                A(level+1).edconn(1,currentCedge) = icount; %coarse grid numbering
                A(level+1).edconn(2,currentCedge) = jcount; %coarse grid numbering
                
                %find the edge connecting i and j on the fine grid
                end_i=-1; end_j=-1;
                for( edgenum=1:size(ea,2) ) %check all fine endges
                    [is loc] = ismember([i;j], ea(:,edgenum), 'rows');
                    if is %if current fine edge is ij or ji
                        if loc(1)==1 %if current fine edge is ij
                            end_i = 1; %ea(1,edgenum);
                            end_j = 2; %ea(2,edgenum);
                        else %current fine edge is ji
                            end_i = 2; %ea(2,edgenum);
                            end_j = 1; %ea(1,edgenum);
                        end
                        break; %stop looking for the edge
                    end
                end
                
                %add the contribution of that edge to the coarse edge molecule
                %sprintf('end_i: %d end_j: %d',end_i,end_j)
                CE(icount,jcount).cc(1,1) = Ee(end_i, end_i, edgenum);
                CE(icount,jcount).cc(1,2) = Ee(end_i, end_j, edgenum);
                CE(icount,jcount).cc(2,1) = Ee(end_j, end_i, edgenum);
                CE(icount,jcount).cc(2,2) = Ee(end_j, end_j, edgenum);
            else %i and j are not directly connected on the fine grid            
                %for each fine selected point k
                fnodes = size(F_set);
                for kcount=1:fnodes(2)
                    k = F_set( kcount ); %get the global point ID for k
                
                    %construct the set of coarse points which are strongly connected to fine point k
                    Sck = intersect( C_set, SC(k).set );
                
                    %if k is connected to both coarse points i and j
                    if ismember( [i j] , Sck )
                        %create a coarse edge between i and j
                        currentCedge = currentCedge + 1;
                        %place the edge in the connectivity matrix
                        A(level+1).edconn(1,currentCedge) = icount; %coarse grid numbering
                        A(level+1).edconn(2,currentCedge) = jcount; %coarse grid numbering
                        
                        %find the edge connecting i and k on the fine grid
                        end_i=-1; end_k=-1;
                        for( edgenum=1:size(ea,2) ) %check all fine endges
                            [is loc] = ismember([i;k], ea(:,edgenum), 'rows');
                            if is %if current fine edge is ik or ki
                                if loc(1)==1 %if current fine edge is ik
                                    end_i = 1; %ea(1,edgenum);
                                    end_k = 2; %ea(2,edgenum);
                                else %current fine edge is ki
                                    end_i = 2; %ea(2,edgenum);
                                    end_k = 1; %ea(1,edgenum);
                                end
                                break; %stop looking for the edge
                            end
                        end                        
                        
                        %add the contribution of that edge to the coarse edge molecule
                        CE(icount,jcount).cc(1,1) = Ee(end_i, end_i, edgenum);
                        CE(icount,jcount).cf(1,2) = Ee(end_i, end_k, edgenum);
                        CE(icount,jcount).fc(2,1) = Ee(end_k, end_i, edgenum);
                        CE(icount,jcount).ff(2,2) = Ee(end_k, end_k, edgenum);
                        
                        %find the edge connecting k and j on the fine grid
                        end_k=-1; end_j=-1;
                        for( edgenum=1:size(ea,2) ) %check all fine endges
                            [is loc] = ismember([k;j], ea(:,edgenum), 'rows');
                            if is %if current fine edge is kj or jk
                                if loc(1)==1 %if current fine edge is kj
                                    end_k = 1; %ea(1,edgenum);
                                    end_j = 2; %ea(2,edgenum);
                                else %current fine edge is jk
                                    end_k = 2; %ea(2,edgenum);
                                    end_j = 1; %ea(1,edgenum);
                                end
                                break; %stop looking for the edge
                            end
                        end                     
                        
                        %add the contribution of that edge to the coarse edge molecule
                        CE(icount,jcount).ff(1,1) = Ee(end_k, end_k, edgenum);
                        CE(icount,jcount).fc(1,2) = Ee(end_k, end_j, edgenum);
                        CE(icount,jcount).cf(2,1) = Ee(end_j, end_k, edgenum);
                        CE(icount,jcount).cc(2,2) = Ee(end_j, end_j, edgenum);                        
                    end
                end
            end                 
        end
    end
end

%generate the coarse edge matrices
for( CEcount=1:size(A(level+1).edconn,2) )
    i=A(level+1).edconn(1,CEcount);
    j=A(level+1).edconn(2,CEcount);
    A(level+1).edges(:,:,CEcount) = CE(i,j).cc - ( CE(i,j).cf * CE(i,j).fc );
end

%generate coarse grid matrix via the Edge Accumulation Method
A(level+1).matrix = AccumulateMatrix( A(level+1).edges , A(level+1).edconn );
tester2 = A(level+1).matrix;

%generate coarse grid matrix via the Galerkin approach
A(level+1).matrix = W(level).weight' * A(level).matrix * W(level).weight;
tester1 = A(level+1).matrix;

%difference = abs( tester1 - tester2 )