function [data,fname]=mphmatrix(m,solname,varargin)
% MPHMATRIX Get model matrices
%  
%   str = MPHMATRIX(MODEL,SOLNAME,'Out',{'A'}) returns a MATLAB structure
%   MA containing the matrix A assembled using the solver node SOLNAME and 
%   accessible as str.A. A being taken from the Out property list.
% 
%   str = mphmatrix(model,solname,fname,'Out',{'A','B',...}) returns a 
%   MATLAB structure str containing the matrices A, B, ... assembled using 
%   the solver node SOLNAME and accessible as str.A and str.B. A and B 
%   being taken from the Out property list.
% 
%   Property    Expression    Description
%   -----------------------------------------------------------------------
%   Out         'K'          Stiffness matrix
%               'L'          Load vector
%               'M'          Constraint vector
%               'N'          Constraint Jacobian
%               'D'          Damping matrix
%               'E'          Mass matrix
%               'NF'         Constraint force Jacobian
%               'NP'         Optimization constraint Jacobian (*)
%               'MP'         Optimization constraint vector (*)
%               'MLB'        Lower bound constraint vector (*)
%               'MUB'        Upper bound constraint vector (*)
%               'Kc'         Eliminated stiffness matrix
%               'Lc'         Eliminated load vector
%               'Dc'         Eliminated damping matrix
%               'Ec'         Eliminated mass matrix
%               'Null'       Constraint null-space basis
%               'Nullf'      Constraint force null-space matrix
%               'ud'         Particular solution ud
%               'uscale'     Scale vector
%   Initmethod  {''} | 'init' | 'sol'   Use linearisation point
%   Initsol     {solname} | 'zero'      Solution to use for linearisation
%   Solnum      {'auto'}     Solution number
% 
%   (*) Require the Optimization module
% 
%   Note that the assembly of the eliminated matrices uses the current 
%   solution vector as scaling method. To get the unscaled eliminated
%   system matrices, it is required to set the scaling method to 'none' in
%   the variables step of the solver configuration node.
%   The load vector is assembled using the current solution available as 
%   linearization point unless the initmethod property is provided. In case
%   of the presence of a solver step node in the solver sequence, the load
%   vector correspond then to the residual of the problem. 
% 
%   Example 1:
%     % Load heat_convection_2d.mph from the Model Library
%     model = mphload('heat_convection_2d.mph');
%     % Deactivate the stationary solver node:
%     model.sol('sol1').feature('s1').active(false);
%     % Extract the stiffness matrix and the load vector:
%     str = mphmatrix(model,'sol1','out',{'K','L'});
%     % Plot the sparsity of the matrix:
%     spy(str.K)
%     % Extract the eliminated system:
%     model.sol('sol1').feature('v1').set('scalemethod', 'none');
%     str = mphmatrix(model,'sol1','out',{'Kc','Lc'});
%     % Compare the sparsity of both system matrices (non-eliminated (b) 
%     % and eliminated one (r)):
%     hold on
%     spy(str.Kc,'r')
% 
%   Example 2:
%     % Load heat_radiation_1d.mph from the Model Library:
%     model = mphload('heat_radiation_1d.mph');
%     % Extract the eliminated residual:
%     str = mphmatrix(model,'sol1','out',{'Lc'});
% 
%   Example 3:
%     % Load heat_convection_2d.mph from the Model Library:
%     model = mphload('heat_convection_2d.mph');
%     % Extract the stiffness matrix and the load vector evaluated for the
%     % initial solution:
%     str = mphmatrix(model,'sol1','out',{'Kc' 'Lc'},'initmethod','init');

% Copyright 2010-2011 COMSOL
