function [out,out1,out2,out3]=fe_case(varargin)

%FE_CASE FEM computation case construction.
%
%	Syntax : Case = fe_case(Case,'PropertyName','Description',PropData)
%                fe_case(model,'command' ...)
%
%       FEM computation cases describe boundary conditions, constraints,
%       loads, and parametric design points. FE_CASE is used to initialize
%       (when Case is not provided as first argument) or modify cases 
%       (Case is provided).
%
%       This function scans trough CASES (see doc('sdt/case')) defined 
%       in model.Stack and builds left hand side arguments needed for 
%       response evaluation.
%       
%	In each CASE, information needed to build the load is given in the
%        .Stack, the result is defined using DOFs declared in CASE.DOF
%
%       Stack entries for you give 'PropertyName','Label',PropData are
%        for boundary conditions
%         'KeepDof', 'FixDof', 'MPC'
%        for loads
%         'DOFLoad', 'DOFSet', 'Fvol', 'FSurf'
%        for physical parameters
%         'par'
%
%	Accepted commands are
%         AddToCase (i)  allows specification of the active case for
%            multiple case models
%         Assemble[...]  calls used to assemble the matrices of a model
%            see on-line documentation with doc('fe_case')
%         GetCase        returns the current case fe_case(model,'getcase')
%         GetDataName    returns data for case entry Name
%         GetT           returns a congruent transformation matrix 
%                        which verifies constraints
%         Reset          removes all entries in the current case
%         Remove         model=fe_case(model,'remove','name') removes entry
%
%	See also help fe_mk
%                doc  sdt/case, fe_case

%	Etienne Balmes
%       Copyright (c) 1996-2006 by SDTools, All Rights Reserved.
%       Use under OpenFEM trademark.html license and LGPL.txt library license

if nargin==1 & comstr(varargin{1},'cvs')
 out='$Revision: 1.44 $  $Date: 2006/05/17 09:16:56 $'; return;
end

if nargin==0&nargout==1
 out=struct('Stack',[],'T',[],'DOF',[]); out.Stack={};
 return;
                             % ('Command',model, ... )
elseif nargin>1 & ischar(varargin{1}) & isfield(varargin{2},'Stack')
 [CAM,Cam]=comstr(varargin{1},1);
 Case=varargin{2}; carg=3;
elseif nargin>1 & ischar(varargin{2}) & isfield(varargin{1},'Stack')
 [CAM,Cam]=comstr(varargin{2},1);
 Case=varargin{1}; carg=3;
elseif ~isstruct(varargin{1}) % ('Command' ...
 Case=fe_case;Case.Stack={}; carg=1;CAM=''; Cam='';
else                          % (model,'Command')
  Case=varargin{1};carg=2; 
  if carg<=nargin [CAM,Cam]=comstr(varargin{carg},1);carg=carg+1;
  else CAM=''; Cam='';end
end

if isfield(Case,'Node')|isfield(Case,'Elt')|~isfield(Case,'T')
     model=Case; if carg<=nargin r1=varargin{carg}; else r1=[];end
     if comstr(Cam,'getcase')
      [out,out1,CAM,Cam,model]=get_case(CAM,Cam,model);
      return;
     elseif comstr(Cam,'getdata')
      [Case,CaseName,st,st,model]=get_case('getcase','getcase',model);
      [CAM,Cam]=comstr(CAM,8);
      if isempty(CAM) CAM=varargin{carg};carg=carg+1;end
      out=stack_get(Case,'',CAM);
      if isempty(out) %warning('Data not found')
      elseif size(out,1)>1 warning('Multiple matches')
      else out = out{3};end
      return;
     elseif comstr(Cam,'setdata')
      [Case,CaseName,st,st,model]=get_case('getcase','getcase',model);
      [CAM,Cam]=comstr(CAM,8);
      Case=stack_set(Case,'',CAM,varargin{carg});carg=carg+1;
      out=stack_set(model,'case',CaseName,Case);
      return;
     elseif isfield(r1,'Stack')&isfield(r1,'T')&isfield(r1,'DOF')
      Case=r1; carg=carg+1; CaseName='';
     else
      [Case,CaseName,st,st,model]=get_case('','',model);
     end

else 
  model=[]; CaseName=''; 
  if comstr(Cam,'getcase')   out=Case; out1='Case 1'; return; end
end
%if carg==0 carg=1;end
if isempty(Case) Case=fe_case; end % default empty case
if isempty(Cam)&carg==1&ischar(varargin{1}) 
 [CAM,Cam]=comstr(varargin{1},1);carg=2;
end

epsl=sdtdef('epsl');

while ~isempty(Cam) % loop on arguments
 

% Add - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
if comstr(Cam,'add') [CAM,Cam]=comstr(CAM,4);

  [Case,CaseName,CAM,Cam,model]=get_case(CAM,Cam,model);

% assemble [m,k,load] - - - - - - - - - - - - - - - - - - - - - - - - - - - -
elseif comstr(Cam,'assemble') [CAM,Cam]=comstr(CAM,9);

if isempty(model)&carg<=nargin
 model=varargin{carg}; carg=carg+1; 
 [Case,CaseName,CAM,Cam,model]=get_case(CAM,Cam,model);
end

if sp_util('issdt')
  eval(sprintf('[out,out1,out2,out3]=fe_caseg(''assemble%s'',model,Case,CaseName,nargout);',CAM));
else
  error('Not implemented yet');
end % OpenFEM

return;

% Case (return a particular case)  - - - - - - - - - - - - - - - -
elseif comstr(Cam,'case')

if length(Cam)>4
     [Case,CaseName,CAM,Cam,model]=get_case(CAM,Cam,model);
end
out=Case; return


% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
elseif comstr(Cam,'get') [CAM,Cam]=comstr(CAM,4);


% 'GetCase' (get case given input) - - - - - - - - - - - - - - - - - -
if comstr(Cam,'c')

 error('You should never get here');

% 'GetT' (get congruent transformation matrix) - - - - - - - - - - - -
elseif comstr(Cam,'t') [CAM,Cam]=comstr(CAM,2);

  if carg<=nargin Case=varargin{carg};carg=carg+1; end

  if ~isfield(model,'DOF')|isempty(model.DOF)
    if isempty(model)&isfield(Case,'DOF')&~isempty(Case.DOF)
      Case.T=speye(length(Case.DOF)); out=Case; return;
    elseif ~isfield(model,'Elt')
     error('Element definitions are needed for case  T matrix building');
    end
    model.DOF=feutil('getdof',model);
  end

  % early return if transformation is defined
  if isfield(Case,'T')&length(model.DOF)==size(Case.T,1)& ...
     isfield(Case,'DOF')&length(Case.DOF)==size(Case.T,2)& ...
     ~comstr(Cam,'new')
       out=Case;
       if comstr(Cam,'dof') out=Case.DOF;end
       return;
  end

  % deal with possible MPC/rigid connections

  DOF=model.DOF; T=speye(length(DOF));
  if carg<=nargin % dof to keep given as argument
   i1=varargin{carg};carg=carg+1;
   if ~isempty(i1) DOF=DOF(i1); T=T(:,i1); end
  end

  RunOpt.fixdof=[];
  for j1=1:size(Case.Stack,1) % loop on things to account for
    switch comstr(Case.Stack{j1,1},-27)
    case 'fixdof'
         r1=Case.Stack{j1,3}; if isstruct(r1) r1=r1.data;end
         if isnumeric(r1) 
         elseif isfield(r1,'data') r1=r1.data;
         elseif ischar(r1); %allow data to contain the selection string
         else error('Not supported FixDof format');   
         end
         if ischar(r1) 
          i1=strfind(comstr(r1,-27),'-dof');
          if ~isempty(i1)
           i2=comstr(r1(i1+4:length(r1)),[-1 0])/100;r1=r1(1:i1-1);
          else i2=0; end
          r1=feutil(['findnode' r1(:)'],model);
          i2=i2(:)';r1=r1(:,ones(size(i2)))+i2(ones(size(r1,1),1),:);
          r1=r1(:);
         end
         RunOpt.fixdof=[RunOpt.fixdof;r1(:)];
     end
  end

  if ~isfield(model,'Node')
  else
   [Case,pdof]=fe_mpc(model,Case,DOF,RunOpt);
   T=Case.T;DOF=Case.DOF;
  end

  for j1=1:size(Case.Stack,1) % loop on things to account for
    switch comstr(Case.Stack{j1,1},-27)
    case 'fixdof'
         i1=fe_c(DOF,RunOpt.fixdof,'ind',2); DOF=DOF(i1); T=T(:,i1);
    case 'dofset'; 
         r1=Case.Stack{j1,3};
         i1=fe_c(DOF,r1.DOF,'ind');
         if isfield(Case,'TIn')&~isempty(Case.TIn);
          error('Multiple DofSet entries are not supported');
         end
         if ~isempty(i1);Case.TIn=T(:,i1); Case.DofIn=DOF(i1);end
         if isfield(r1,'def');Case.TIn=Case.TIn*r1.def;end
         i1=fe_c(DOF,r1.DOF,'ind');if ~isempty(i1);DOF(i1)=[];T(:,i1)=[];end
         if isempty(fe_c(model.DOF,r1.DOF,'ind'))
           sdtw('DofSet %s entry affects no free DOF',Case.Stack{j1,2})
         end  
    case {'keepdof','par','rigid','dofload','fsurf','fvol','sensdof', ...
          'sensstrain','cyclic','rbe3','mpc','nastran','info'}
    otherwise sdtw('_nb','%s not supported by fe_case GetT',Case.Stack{j1,1});
    end % of supported case
  end
  Case.T=T; Case.DOF=DOF; out=Case;

  if comstr(Cam,'dof') out=Case.DOF;end
  return;

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
else error('Not a known get command'); end

% Cyclic  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
elseif comstr(Cam,'cyclic') [CAM,Cam]=comstr(CAM,7);

   if comstr(Cam,'default')
    name='new cyclic';adof=[];
   else
    if carg>nargin error('You must specify some data'); end
    name=varargin{carg}; carg=carg+1;
    if ~ischar(name) error('You must specify a name of  FixDof case entry');end
    data=varargin{carg}; carg=carg+1;
   end
   Case=stack_set(Case,'cyclic',name,data); Case.T=[];

% UN=0  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
elseif comstr(Cam,'un=0') [CAM,Cam]=comstr(CAM,7);

    name=varargin{carg}; carg=carg+1;
    data=varargin{carg}; carg=carg+1;
    
    r1=data.normal; r1=r1./(max(abs(r1),[],2)*[1 1 1]);
    ind=find(any(abs(r1+1)<sqrt(eps),2));r1(ind,:)=-r1(ind,:);
   r1
   DOF=[data.ID(:)+.01;data.ID(:)+.02;data.ID(:)+.03];
   data=struct('DOF',DOF,'c',[diag(r1(:,1)) diag(r1(:,2)) diag(r1(:,3))]);
   Case=stack_set(Case,'mpc',name,data); Case.T=[];


% FixDOF  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
elseif comstr(Cam,'fixdof') [CAM,Cam]=comstr(CAM,7);

   if comstr(Cam,'default')
    name='new FixDof';adof=[];
   else
    if carg>nargin error('You must specify some data'); end
    name=varargin{carg}; carg=carg+1;
    if ~ischar(name) error('You must specify a name of  FixDof case entry');end
   adof=varargin{carg}; carg=carg+1;
   end
   if isnumeric(adof); adof=adof(:);end
   if ~isempty(Case.Stack) 
    i1=strmatch('FixDof',Case.Stack(:,1),'exact');
    if ~isempty(i1) Case.Stack(i1,1)={'FixDof'};end
   end
   if ~isempty(Case.Stack) 
     i1=strmatch('fixdof',comstr(Case.Stack(:,1),-27),'exact');
     if ~isempty(i1) Case.Stack(i1,1)={'FixDof'};end
   end
   if ischar(adof);
         i1=strfind(comstr(adof,-27),'-id'); if ~isempty(i1);
             [ID,i2,i3,i4]=sscanf(adof(i1+3:end),'%i');
             adof(i1:i1+1+i4)=''; adof=comstr(adof,1);
             adof=struct('data',adof,'ID',ID);
         end
   end
   Case=stack_set(Case,'FixDof',name,adof); Case.T=[];

% AutoSPC - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
elseif comstr(Cam,'autospc') 

 eval('Case=fe_caseg(CAM,model,Case);'); % SDT extension

% KeepDof - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
elseif comstr(Cam,'keepdof') [CAM,Cam]=comstr(CAM,8);

   if comstr(Cam,'default')
    name='new KeepDof';adof=[];
   else
    if carg>nargin error('You must specify some data'); end
    name=varargin{carg}; carg=carg+1;
    if ~ischar(name) error('You must specify a name of a KeepDof case entry');end
    adof=varargin{carg}; carg=carg+1;
   end
   if ~isempty(Case.Stack) 
    i1=strmatch('keepdof',comstr(Case.Stack(:,1),-27),'exact');
    if ~isempty(i1) Case.Stack(i1,1)={'KeepDof'};end
   end

   Case=stack_set(Case,'KeepDof',name,adof(:));Case.T=[];

% sens SDT extension - - - - - - - - - - - - - - - - - - - - - - -
elseif comstr(Cam,'sens')

 eval('[Case,carg]=fe_caseg(CAM,model,Case,varargin(carg:end),carg);')
 if ischar(carg);eval(carg);return; end

% MPC  - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
elseif comstr(Cam,'mpc'); [CAM,Cam]=comstr(CAM,11);
 
   if comstr(Cam,'default')
    name='new MPC';adof=[];
   else
    if carg>nargin error('You must specify some data'); end
    name=varargin{carg}; carg=carg+1;
   end
    
   if ~ischar(name) error('You must specify a name of MPC case entry');end
   adof=varargin{carg}; carg=carg+1;
   if isa(adof,'cell')&size(adof,2)==3
    adof=struct('SID',adof{1},'DOF',adof{2},'c',adof{3});
   elseif isfield(adof,'c')&isfield(adof,'DOF')
   else error('Not a proper MPC format');
   end
   Case=stack_set(Case,'mpc',name,adof);Case.T=[];
 
% Rigid- - - - - - - - - - - - - - - - - - - - - - - - - - - - -
elseif comstr(Cam,'rigid'); [CAM,Cam]=comstr(CAM,6);

 if comstr(Cam,'append')
  error(1)
 else
   if comstr(Cam,'default')
    name='Rigid';adof=[];
   else
    if carg>nargin error('You must specify some data'); end
    name=varargin{carg}; carg=carg+1;
    if ~ischar(name) error('You must specify a name of a Rigid case entry');end
    adof=varargin{carg}; carg=carg+1;
   end

   Case=stack_set(Case,'rigid',name,adof);Case.T=[];

 end

% DOFLoad  - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
elseif comstr(Cam,'dofload'); [CAM,Cam]=comstr(CAM,8);
 
   if comstr(Cam,'default')
    name='new dofload';adof=struct('name',name,'DOF',[],'def',[]);
   else
    if carg>nargin error('You must specify some data'); end
    name=varargin{carg}; carg=carg+1;
    if ~ischar(name) error('You must specify a name of DOFLoad case entry');end
    adof=varargin{carg}; carg=carg+1;
   end

   if isfield(adof,'DOF') & isfield(adof,'def')
    adof.name=name;
    if size(adof.DOF,1)~=size(adof.def,1)
       error('.DOF and .def fields not consistent with DOF load case');
    end
    Case=stack_set(Case,'DOFLoad',name,adof);
   elseif isnumeric(adof)&size(adof,2)==1
    sdtw('_nb','Assuming unit DofLoad entries');
    r1=struct('DOF',adof,'def',eye(size(adof,1)), ...
     'name',name,'lab',{fe_curve('doflab load',adof)});
    Case=stack_set(Case,'DOFLoad',name,r1);

   else
     error('For DOFLoad data.DOF data.def fields must be defined')
   end

% DOFSet  - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
elseif comstr(Cam,'dofset') [CAM,Cam]=comstr(CAM,7);
 
   if comstr(Cam,'default')
    name='new DOFSet';adof=struct('name',name,'DOF',[],'def',[]);
   else
    if carg>nargin error('You must specify some data'); end
    name=varargin{carg}; carg=carg+1;
    if ~ischar(name) error('You must specify a name of DOFSet case entry');end
    adof=varargin{carg}; carg=carg+1;
   end

   if isfield(adof,'DOF') & isfield(adof,'def')
    adof.name=name;
    if size(adof.DOF,1)~=size(adof.def,1)
       error('.DOF and .def fields not consistent with DOF load case');
    end
    Case=stack_set(Case,'DOFSet',name,adof);
   elseif isnumeric(adof)&size(adof,2)==1
    sdtw('_nb','Assuming unit DofSet entries');
    if any(rem(adof,1)==0); % eliminate wild cards
      if isfield(model,'DOF')&~isempty(model.DOF);DOF=model.DOF;
      else; DOF=feutil('getdof',model);
      end
      adof=fe_c(DOF,adof,'dof');
    end
    r1=struct('DOF',adof,'def',eye(size(adof,1)),'name',name);
    Case=stack_set(Case,'DOFSet',name,r1);
   else
     error('For DOFSet data.DOF data.def fields must be defined')
   end

% Fvol  - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
elseif comstr(Cam,'fvol') [CAM,Cam]=comstr(CAM,5);
 
   if comstr(Cam,'default')
    name='new FVol';
    r1=struct('name',name,'sel','GroupAll','dir',{{'0','0','0'}});
   else
    if carg>nargin error('You must specify some data'); end
    name=varargin{carg}; carg=carg+1;
    if ~ischar(name) error('You must specify a name of Fvol case entry');end
    r1=varargin{carg}; carg=carg+1;
    r1.name=name;
   end

   if ~isfield(r1,'sel')
    error('You must specify .sel (element selection) for volume force');
   end
   % either a def/DOF or dir array
   if isfield(r1,'def')&isfield(r1,'DOF');
   elseif ~isfield(r1,'dir'); r1.dir={{'1','0','0'}};
   end

   Case=stack_set(Case,'FVol',name,r1);

% FSurf  - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
elseif comstr(Cam,'fsurf') [CAM,Cam]=comstr(CAM,6);

   if comstr(Cam,'default')
    name='new FSurf';
    r1=struct('name',name,'presel','','sel','','def',[1;1;1],'DOF',[1:3]/100);
   else
    if carg>nargin error('You must specify some data'); end
    name=varargin{carg}; carg=carg+1;
    if ~ischar(name) error('You must specify a name of FSurf case entry');end
    r1=varargin{carg}; carg=carg+1;
    r1.name=name;
   end

   if ~isfield(r1,'sel')&~isfield(r1,'set')
    error('You must specify an face selection for surface force definition');
   end
   if ~isfield(r1,'DOF') r1.DOF=[];end
   Case=stack_set(Case,'FSurf',name,r1);

% Par (SDT parameters) - - - - - - - - - - - - - - - - - - - - - - - - - - -
elseif comstr(Cam,'par'); [CAM,Cam]=comstr(CAM,4);

   if ~isempty(strfind(Cam,'default'))
    if comstr(Cam,'m'); name='new mass parameter';
         r1=struct('sel','groupall','coef',[2  1 .5 2 1]);
    else; name='new stiffness parameter';
         r1=struct('sel','groupall','coef',[1  1 .5 2 1]);
    end
   else
     name=varargin{carg};carg=carg+1;
     [CAM,Cam]=comstr(varargin{carg},1);carg=carg+1;
     r1=struct('sel','groupall','coef',[1 1 .5 2 1],'name',name);
     if comstr(Cam,'-k'); r1.coef(1)=1; [CAM,Cam]=comstr(CAM,3);end
     if comstr(Cam,'-m'); r1.coef(1)=2; [CAM,Cam]=comstr(CAM,3);end
     if ~isempty(CAM); r1.sel=CAM;end
   end
   Case=stack_set(Case,'par',name,r1);

% Remove - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
elseif comstr(Cam,'remove')

  st=varargin{carg};carg=carg+1;
  Case=stack_rm(Case,'',st);

% Info - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
elseif comstr(Cam,'info')

  fprintf('\nCase : %s\n\n',CaseName)
  disp(Case)
  disp(Case.Stack)
  if nargout==0 clear out; return; 
  else out=Case; end

% SetCase - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
elseif comstr(Cam,'set')

  [Case,CaseName,CAM,Cam,model]=get_case(CAM,Cam,model);
  Case=fe_case;

% Reset - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
elseif comstr(Cam,'reset')

  Case=fe_case;

% Connection - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
elseif comstr(Cam,'connection')|comstr(Cam,'rbe3') % advanced SDT connections
  name=varargin{carg}; carg=carg+1;
  eval('model=fe_caseg(sprintf(''%s -name "%s"'',CAM,name),model,varargin{carg:end});');

else error('Not a supported CASE entry'); end

if carg<=nargin [CAM,Cam]=comstr(varargin{carg},1);carg=carg+1;
else Cam='';end

end % loop on input arguments - - - - - - - - - - - - - - - - - - -

if ~isempty(model)&~isempty(CaseName)
  model=stack_set(model,'case',CaseName,Case); out=model;
else out=Case;
end


% ------------------------------------------------------------------------
% ------------------------------------------------------------------------
% ------------------------------------------------------------------------

function  [Case,CaseName,CAM,Cam,model]=get_case(CAM,Cam,model);

if length(Cam)<4 i1=[];else i1=strfind(Cam,'case');end
i2=1; Case=stack_get(model,'case'); CaseName='';

if ~isempty(i1) % the keywork case is present  
    [i2,st1,st2]=comstr(Cam(i1:end),'case','%i');
    if isempty(i2) & comstr(st2,'case') % Deal with GetCaseCase(i)
     [i2,st1,st2]=comstr(st2,'case','%i');
    end
    if isempty(i2)&length(Cam)>i1+3
      [CaseName,CAM,Cam]=comstr(CAM(i1:end),'case','%s');
      if ~isempty(Case)&any(strcmp(CaseName,Case(:,2)))
        i2 =find(strcmp(CaseName,Case(:,2)));
      else i2=size(Case,1)+1; end
    else
      [CAM,Cam]=comstr([CAM(1:i1-1) st1],1);
      if isempty(i2) i2=1;end
    end
else i2=1;
end

if i2>size(Case,1)
  if i2>1&isempty(CaseName)
    warning(sprintf('Case %i undefined, using Case 1',i2));
  end
  Case=fe_case;if isempty(CaseName) CaseName='Case 1';end
else CaseName=Case{i2,2}; Case=Case{i2,3};end

if isempty(CaseName) CaseName=sprintf('Case %i',i2); end
if isempty(Case) Case=fe_case; end

