xGenerateCoarseSystem.m

function cs = xGenerateCoarseSystem(g,rock,cg,overlap,...
                                    bc,src,weighting,mob,...
                                    coX,coiG,coFaces)

% Constructs coarse system component matrices from fine grid model in
% parallel.
%
% PARAMETERS:
% 
% - G -- Composite variable. Grid structure from MRST.
%
% - rock -- Composite variable. Rock structure from MRST.
% 
% - CG -- composite variable. Coarse grid structure from MRST.
% 
% - overlap -- Composite variable. Number of fine-grid cells in each 
%     physical direction with which to extend the supporting domain of any
%     given basis functions.
%
% - bc -- Composite variable. Boundary condition structure from MRST.
%
% - src -- Composite variable. Explicit source contributions as defined by
%     'addSource'. 
%     Use src = [] as default value.
% 
% - weighting -- Composite variable. Basis function driving source term as
%     supported by function 'evalBasisSource'.
%     use weighting = 'perm' as default value.
%
% - mob -- Composite variable. Total mobility. One scalar value for each
%     cell in the underlying (fine) model.
%
% - coX -- Distributed cell array. Local part on each worker contains all
%   inner products that will be needed on that worker to construct basis
%   functions.
%
% - coiG -- Distributed cell array. One cell pr. worker. that is, for worker 
%   w we have that coiG{w}{i} == k means that coX{w}{i} is inner product 
%   of global cell k.
%
% - coFaces -- Distributed array. Face indices.
%   running: 
%     lpFaces = getLocalPart(coFaces);
%   on a worker, means that lpFaces(:) will be calculated on that worker.
%   (An that the needed ips are available when this function has been
%   run.)
%
% RETURNS:
% - CS -- Composite struct. Value is only stored on first lab. That is you
%     can get the contents of CS on host by running
%       hostCS = CS{1};
%     As in MRST CS contains the following fields:
%        - basis   - Flux basis functions as generated by function
%                   evalBasisFunc.
%        - basisP  - Pressure basis functions as generated by function
%                   evalBasisFunc.
%        - C       - C in coarse hybrid system.
%        - D       - D in coarse hybrid system.
%        - sizeC   - size of C (== size(CS.C))
%        - sizeD   - size of D (== size(CS.D))
%
% NOTE:
% not implemented for utilize evalGlobalBasisFunc: need isempty(global_inf) 
% this is supported in generateCoarseSystem in MRST.
%
% EXAMPLE:
%   See xExample for a complete example of usage. 
%  
% SEE ALSO:
%   xExample, xInitWorkers, xBroadcast, xComputeMimeticIP, 
%   xDistributeIP, xEvalBasisFunc

%{
A part of the xmsmfem module for MRST:
http://www.sintef.no/Projectweb/MRST/
Adapted from the msmfem module with the Parallel Computing Toolbox

Released under the GNU General Public License:
http://www.gnu.org/licenses/gpl.html
 
Written by
Anders Hoff 2012
http://master.andershoff.net
%}


spmd
  weight = evalBasisSource(g, weighting, rock);
end

[coV, coP] = xEvalBasisFunc(g,cg,coX,coiG,coFaces,weight,...
                            mob, src, bc,overlap);
spmd
  V = gather(coV,1);
  P = gather(coP,1);
  
  % included in CS.
  activeFaces = gather(coFaces,1);

  if labindex == 1
    cs.basis  = cell([cg.faces.num,1]);
    cs.basisP = cell([cg.faces.num,1]);
    cs        = assignBasisFuncs(cs, V, P);

    % 1) Compute coarse grid matrices C and D.

    % Compute sizes for matrices 'B', 'C', and 'D'.  Includes all faces, even
    % faces for which there are no associated degrees of freedom.
    sizeB  = size(cg.cells.faces, 1) * [1, 1];
    sizeC  = [sizeB(1), double(cg.cells.num)];
    sizeD  = [sizeB(1), double(cg.faces.num)];

    % Compute topology matrices (C and D).
    cellNo = rldecode(1:cg.cells.num, diff(cg.cells.facePos), 2).';
    cs.C   = topo_mat(cellNo,sizeC(2));
    cs.D   = topo_mat(cg.cells.faces(:,1),sizeD(2));
    
    % 2) Define degrees of freedom and basis function weighting scheme.

    cs.basisWeighting  = weighting;
    
    cs.activeFaces = activeFaces;
    cs.activeCellFaces = find(sum(cs.D(:,activeFaces),2));

    % 3) Assign system matrix sizes for ease of implementation elsewhere.
    % cs.type  = S.type; % hardcode? FIX
    cs.sizeB = sizeB;
    cs.sizeC = sizeC;
    cs.sizeD = sizeD;
  end
end
end

function CS = assignBasisFuncs(CS,V,P) % from MRST
f = cellfun(@(x) x{3}, V);
CS.basis (f) = V; CS.basisP(f) = P;
end

%-----------------------------------------------------------------------
function res = topo_mat(j,n) % from MRST
res = sparse(1:numel(j), double(j), 1, numel(j), n);
end