function [out,out1,out2]=fe_mknl(varargin)

%FE_MKNL model assembly optimized for non-linear operation
%
% Typically for non-linear operation
%
%   Case=fe_mknl('init',model);                 % performed once
%   mat=fe_mknl('assemble',model,Case,MatType); % performed at every step
%
% For linear matrix assembly. Three calling formats are accepted
%   model=fe_mknl(model)
%   [k,mdof]=fe_mknl(model);
%   [m,k,mdof]=fe_mknl(model);
%
%  The init phase generates Case.GroupInfo which gives
%    GroupInfo = {DofPos Pointers Integ Constit gstate ElMap ...
%       InfoAtNode EltConst}
%    Pointers has one column per element giving %
%     [OutSize1 OutSize2 u3 NdNRule MatDes IntegOffset ConstitOffset ...
%       StateOffset]
%

%	E. Balmes, J. Leclere, C. Delforge
%       Copyright (c) 2001-2006 by INRIA and SDTools,All Rights Reserved.
%       Use under OpenFEM trademark.html license and LGPL.txt library license
%       Use fe_mknl('cvs') for revision information


if nargin==0; help fe_mknl
elseif isstruct(varargin{1}) % Assemble the mass and stiffness

 model=varargin{1}; carg=2;
 if isa(model,'v_handle'); model=model.GetData;end
 if carg<=nargin; Cam=varargin{carg};carg=carg+1;else; Cam='';end
 if isstr(Cam); st=horzcat('assemble',Cam);else; st='assemble';end
 if ~isempty(strfind(comstr(st,-27),'not')); RunOpt.NoT=1;else; RunOpt.NoT=0;end
 if ~isempty(strfind(Cam,'keep')); st1='initKeep';else st1='init';end
 if RunOpt.NoT; st1=[st1 'NoT'];end
 [Case,model.DOF]=fe_mknl(st1,model);
 m=fe_mknl(st,model,Case,2);
 k=fe_mknl(st,model,Case,1);

 model.K={m,k};  if ~isfield(model,'Opt'); model.Opt=1;end
 model.Opt(2,1:2)=[2 1]; model.Opt(2,3:size(model.Opt,2))=0;
 if nargout<2 out=model;
 elseif nargout==2; out=k; 
   if RunOpt.NoT;out1=model.DOF; else;out1=Case.DOF;end
 elseif nargout==3; out=m; out1=k;  
   if RunOpt.NoT;out2=model.DOF; else;out2=Case.DOF;end
 end
 return; 

else
 [CAM,Cam]=comstr(varargin{1},1);carg=2;
 if carg<=nargin; model=varargin{carg};carg=carg+1;
 elseif comstr(Cam,'cvs')
  out='$Revision: 1.92 $  $Date: 2006/05/15 15:10:47 $';
  return;
 end
 if isa(model,'v_handle'); model=model.GetData;end
end

% -----------------------------------------------------------------------
% InitOptimProduct : init and optimize numbering for product evaluations
if comstr(Cam,'initoptimproduct'); [CAM,Cam]=comstr(CAM,9);

 [Case,model.DOF]=fe_mknl('init',model);
 k=fe_mknl('assemble',model,Case,1);
 [i1,i2,st]=fe_time('optimproduct',k);
 model=stack_set(model,'info','FastestProduct',{i2,st});
 Case.DOF=Case.DOF(i1);Case.T=Case.T(:,i1);
 r1=stack_get(model,'case');
 mode=stack_set(model,'case',r1{2},Case);
 %[Case,model.DOF]=fe_mknl('init',model,Case);
 out=model; out1=Case;

% -----------------------------------------------------------------------
% Initialize all the info needed for assembly 
elseif comstr(Cam,'init'); [CAM,Cam]=comstr(CAM,5);


RunOpt=struct('LastError','','Gstate',0,'InitFailed',0);

% Build cell array with 
% DofPos Pointers IntegOpt Constit State ElMap
[EGroup,nGroup]=getegroup(model.Elt);
Case.GroupInfo=cell(nGroup,8);

i1=strfind(Cam,'-gstate');if ~isempty(i1);
 RunOpt.Gstate=1; CAM(i1+[0:6])=''; [CAM,Cam]=comstr(CAM,1);
end
% Phase 1  : generation of the full list of DOFS - - - - - - - - - - - -

if ~isempty(strfind(Cam,'keep'))
 i1=strfind(Cam,'keep');CAM(i1+[0:3])='';[CAM,Cam]=comstr(CAM,1);
 [mdof,r1,nw,model.Elt,eltid]=feutil('getdof',model); 
else
 [model.DOF,r1,nw,model.Elt,eltid]=feutil('getdof',model);
 mdof=model.DOF;  i1=find(mdof>0);
end
if carg>nargin; Case=fe_case(model,'gett');
else; Case=varargin{carg};carg=carg+1;
  if size(Case.T,1)~=length(model.DOF)|size(Case.T,2)~=length(Case.DOF)
   Case=fe_case(model,'gett',Case);
  end
end


mdof=model.DOF;  i1=find(mdof>0);
if isempty(i1) nd=[];
else
 nd=sparse(round(rem(mdof(i1),1)*100),fix(mdof(i1)),1:length(i1),100, ...
    max([fix(max(mdof(i1)))+1 max(model.Node(:,1))]));
 N=length(model.DOF);if max(max(nd))>N error('Something wrong with MDOF'); end
end

i1=find(model.DOF<0);if isempty(i1) nde=[];
else % element DOFs
   i2=-round(model.DOF(i1)*1000);
   i2=[rem(i2,1000) fix(i2/1e3)];
   i3=max(i2(:,2)+1);if i3>fix(2^31/1000); error('EltId too large');end
   nde=sparse(i2(:,1),i2(:,2),i1,1000,max(i2(:,2)+1));
end

% Phase 2 get PL,IL and global nodes

pl=fe_mat('getpl',model); il=fe_mat('getil',model);
[node,bas,NNode,Case.cGL]=feutil('getnodebas',model);
Case.Node=node; if ~isempty(bas); Case.bas=bas;end
conn=int32(zeros(size(model.Elt,2),size(model.Elt,1)));

for jGroup=1:nGroup  % Loop on groups

% Phase 3 Build DofPos connectivity table 

  [ElemF,opt]= feutil('getelemf',model.Elt(EGroup(jGroup),:),jGroup);
  EGID=opt(1); cEGI = EGroup(jGroup)+1:EGroup(jGroup+1)-1;
  [eCall,opt]= fe_super('call',ElemF,0);
  try; i1=feval(ElemF,'dofcall');catch;i1=[];end
  if ischar(i1);  eval(i1); % Catch for in element building of DofPos
  else; % standard methodology with fixed number of DOF per elt
   i1=fe_super('dof',ElemF);
   if opt(1)==1 % unique superelement
    i4=zeros(length(i1),1);
    i3=find(i1>0); i4(i3)=full(nd(round(i1(i3)*100)-100))-1;
    i3=find(i1<0); if ~isempty(i3)
      i4(i3)=full(nde(round(-i1(i3)*1000-1000)))-1;
    end
    i2=int32(i4); cEGI=[];
   elseif ~any(i1<0)
    i2=round(rem(i1(:),1)*100-100);
    i2=model.Elt(cEGI,fix(i1))'*100+i2(:,ones(length(cEGI),1));
    i2(find(i2<1))=99;  % this allows for variable node numbers (zeros)
    i2=int32(full(nd(i2))-1);
   elseif opt(1)==2 % generic superelement case with element DOFs
    i4=zeros(length(i1),length(cEGI));
    i3=find(i1(:)>0); i2=round(rem(i1(i3),1)*100-100);
    i2=model.Elt(cEGI,fix(i1(i3)))'*100+i2(:,ones(length(cEGI),1));
    i2(find(i2<1))=99;  % this allows for variable node numbers (zeros)
    i4(i3,:)=full(nd(i2))-1;
    i3=find(i1(:)<0);  i2=round(-i1(i3)*1000-2000);
    i2=ones(length(i3),1)*(eltid(cEGI)*1000)'+i2(:,ones(length(cEGI),1));
    i4(i3,:)=full(nde(i2))-1;
    i2=int32(i4);
   end
  end % local implementation or external approach
  Case.GroupInfo{jGroup,1}=int32(i2); % DofPos
  Case.DofPerElt(jGroup)=size(i2,1);% size of DofPos
  if ~isempty(i2)&~isempty(cEGI)
    conn(1:size(i2,1),cEGI)=i2;
  end

% Phase 4 Get MatId and ProId and call elem('integinfo',[MatId ProId],pl,il)
% to generate integ and constit

  if EGID>0
  % Init IntegOpt and Constit, one COLUMM per property
  i1=fe_super('prop',ElemF);
  if all(i1(1:2)) 
   i2=find(i1(1:2)<=size(model.Elt,2));
   if isempty(i2); i1=[];
   else; i1=model.Elt(cEGI,i1(i2)); [i1,i2,i3]=unique(i1,'rows');
   end
  else i3=1; % xxx what to do when some prop don't exist
  end

  integ=i1'; constit=zeros(1,size(integ,2));  gstate=[];elmap=[];
  for j1=1:size(integ,2)
   SE=[];RunOpt.LastError='';
   try
    [r1,i1,elmap]=feval(ElemF,'integinfo',integ(1:2,j1),pl,il);
   catch
    RunOpt.LastError=lasterr;
    try; 
     SE=stack_get(model,'SE'); if ~isempty(SE) % SE defined in stack
      for j2=1:size(SE,1); 
        eval(sprintf('global SE%s;SE%s=SE{j2,3};',SE{j2,2},SE{j2,2}));
      end
     end
     SE=whos('global',sprintf('SE%s',ElemF));
     if ~isempty(SE);eval(sprintf('global SE%s;SE=SE%s;',ElemF,ElemF));end
     if isfield(SE,'Opt')&isfield(SE,'K'); r1=[];i1=[];elmap=[];
     elseif ~exist(ElemF)
      RunOpt.LastError=sprintf( ...
       'Group %i (%s) not a valid superelement',jGroup,ElemF);
      error(RunOpt.LastError);
     else; error(RunOpt.LastError); end
    catch; 
     if sdtdef('diag'); 
      if isempty(SE)
        [r1,i1,elmap]=feval(ElemF,'integinfo',integ(1:2,j1),pl,il);
      else; error(RunOpt.LastError);
      end
     else; 
      RunOpt.InitFailed=1;
      sdtw('_nb', ...
      '%s(''IntegInfo'') (build integ,constit) failed (group %i)',ElemF,jGroup);
      sdtw('_nb',RunOpt.LastError); r1=[];i1=[];elmap=[];
     end
    end
   end
   if ~isempty(elmap)&size(elmap,1)~=size(Case.GroupInfo{jGroup,1})
     error(sprintf('ElMap and DofPos are inconsistent for group %i',jGroup))
   end
   if any(r1)|~any(integ(1:2,j1))|any(strcmp(ElemF,{'celas','mass2'}))
   elseif isfield(SE,'Opt')&any(SE.Opt(1)==[1 2])
   else
    sdtw('_nb','Null constit for [MatId %i ProId %i]',integ(1:2,j1));
   end
   integ(1:size(i1,1),j1)=i1;
   constit(1:size(r1,1),j1)=r1;
   r2(1:5,j1)=[0 0 j1 j1 0]';
  end
  r2(3,:)=(r2(3,:)-1)*size(integ,1);
  r2(4,:)=(r2(4,:)-1)*size(constit,1);
  pointers([1 2 6 7 10],1:length(i3))=r2(:,i3);

% Phase 5 let other initializations be performed for various element types

  eCall='';if isempty(SE);try; 
   eCall=feval(ElemF,'GroupInit',integ,constit);
   eval(eCall);
  catch; 
   if sdtdef('diag')>11; eval(eCall);
   elseif sdtdef('diag')
    st=sprintf( ...
     'Error %s\n during GroupInit use sdtdef(''diag'',12) for debug', ...
     lasterr);error(st)
   elseif ~isempty(eCall); 
    sdtw('_nb','%s(''GroupInit'') failed (group %i)',ElemF,jGroup);
   end
  end;end

  %% initialize gstate for heart model: epsilon_c k tau on each jW %% 
  if size(Case.GroupInfo,2)==8 
  if (isfield(Case.GroupInfo{jGroup,8},'material')...
      & strfind(Case.GroupInfo{jGroup,8}.material,'heart'))
    gstate=zeros(7+Case.DofPerElt(jGroup),Case.GroupInfo{jGroup,8}.Nw*size(Case.GroupInfo{jGroup,1},2));
    
  end
  end

  Case.GroupInfo{jGroup,6}=elmap;
  Case.GroupInfo{jGroup,2}=pointers;
  Case.GroupInfo{jGroup,3}=int32(integ);
  Case.GroupInfo{jGroup,4}=constit;
  Case.GroupInfo{jGroup,5}=gstate;
 end % EGID>0
  i1=fe_super('node',ElemF); 

end % jGroup
if size(Case.GroupInfo,2)<8 Case.GroupInfo{1,8}=[];end
if comstr(Cam,'nocon');
 if RunOpt.InitFailed;Case.InitFailed=1;end
elseif RunOpt.InitFailed==0
 Case.MatGraph=of_mk('meshGraph', ...
	      size(mdof,1),nGroup,size(conn,1), ...
	      int32(EGroup),int32(diff(EGroup)-1),int32(Case.DofPerElt), ...
	      conn );
end

if RunOpt.Gstate % initialize Group state information
 model.DOF=mdof;[Case,dc]=elem0('initgstate',model,Case);
 out=Case; out1=mdof;out2=dc;
else % Third argument is for use in fe_mk
 out=Case; out1=mdof; out2={nd,nde,nw,eltid,node,bas};
end % Group state information


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

% model,MatDes  or model,Case,MatDes or model,Case,def,Matdes
% or model,Case,MatDes,def
if carg>nargin; Case=fe_case(model,'getcase');
else; Case=varargin{carg};carg=carg+1;
end
if ~isstruct(Case); MatDes=Case;Case=fe_case(model,'getcase'); 
elseif carg>nargin; error('You must specify a desired matrix type');
else; MatDes=varargin{carg};carg=carg+1;
end
if isstruct(MatDes); def=MatDes;MatDes=varargin{carg};carg=carg+1;
elseif carg<=nargin; def=varargin{carg};carg=carg+1;
else; def=struct('def',[],'DOF',model.DOF);
end

[EGroup,nGroup]=getegroup(model.Elt);
if ~isfield(Case,'GroupInfo'); Case=fe_mknl('init',model,Case); end
NNode=sparse(Case.Node(:,1),1,1:size(Case.Node,1));

jMat=1; % Provision for when multiple outputs will be supported
% energy computations - - - - - - - - - - - - - - - - - - - - - - -
if any(MatDes(jMat)==[250:299]) 
 
 k=zeros(size(model.Elt,1),size(def.def,2));
 OutType='ener'; MatDes(jMat)=MatDes(jMat)-250; 
 RunOpt.Ener=zeros(size(model.Elt,1),size(def.def,2));

% 300-400 modify internal state in Case
elseif MatDes(jMat)>=300&MatDes(jMat)<400

 OutType='state'; 
 if isfield(def,'def'); FieldAtDof=def.def; else; FieldAtDof=def;end 
 
% RHS right hand side inits - - - - - - - - - - - - - - -
elseif MatDes(jMat)>=100&MatDes(jMat)<200  

 OutType='rhs'; if carg<=nargin 
    NNode=varargin{carg}; carg=carg+1;
 end
 o1=zeros(size(model.DOF,1),1);

% Standard matrix inits - - - - - - - - - - - - - - - - - - - - - - - -
else

 if ~isfield(Case,'MatGraph')|isempty(Case.MatGraph)
   error('Case.MatGraph is not filled, the fe_mknl(''init'') phase failed');
 end
 k=Case.MatGraph; of_mk( 'fillvalue', k, 0 );
 OutType='mat';

end

% these are standard variable names as stated in elem0
pl=fe_mat('getpl',model); il=fe_mat('getil',model);
elt=model.Elt; node=Case.Node;
t=cputime; st1='';if isfield(Case,'bas'); bas=Case.bas;end


for jGroup=1:nGroup  % Loop on groups

  [ElemF,i1]= feutil('getelemf',model.Elt(EGroup(jGroup),:),jGroup);
  EGID=i1(1); 
  inode=fe_super('node',ElemF);

  % See documentation elem0#call
  DofPos=Case.GroupInfo{jGroup,1};
  pointers=Case.GroupInfo{jGroup,2};
  integ=Case.GroupInfo{jGroup,3};
  constit=Case.GroupInfo{jGroup,4};
  gstate=Case.GroupInfo{jGroup,5};
  elmap=int32(Case.GroupInfo{jGroup,6});
  InfoAtNode=Case.GroupInfo{jGroup,7};
  EltConst=Case.GroupInfo{jGroup,8};
  switch OutType
  case {'mat','mat_of','mat_og'}
    try
      [fHandle,SymFlag]=feval(ElemF,'matcall',integ,constit);
    catch
      SymFlag=1; fHandle=feval(ElemF,'call',integ,constit);
    end
    pointers(5,:)=MatDes(jMat);
    if isempty(fHandle) fHandle=feval(ElemF,'call'); end
    if ischar(fHandle)&strcmp(fHandle,'mat_of')  % standard of_mk call
       OutType='mat_of';SymFlag=1;
       if length(constit)<2; error('Initialization failed');end
    elseif ischar(fHandle)&strcmp(fHandle,'mat_og')  % of_mk matrix integration
       OutType='mat_og';SymFlag=0;
       if length(constit)<2; error('Initialization failed');end
    else  % m file elements returning full element matrix
      fullmap=find(ones(Case.DofPerElt(jGroup),Case.DofPerElt(jGroup)));
      fullmap=int32(reshape(fullmap, ...
         Case.DofPerElt(jGroup),Case.DofPerElt(jGroup)));
      OutType='mat';
    end
    opt=MatDes; % obsolete kept for compatibility with old elements
    if ~isequal(model.DOF,Case.DOF)&~isempty(def.def)&isequal(def.DOF,Case.DOF)
     def.def=Case.T*def.def; def.DOF=model.DOF;
     if MatDes(1)==5; 
      sdtw('_nb','def resizing is not acceptable for NL rhs computations')
     end
    end

  case 'ener'

    try
      [fHandle,SymFlag]=feval(ElemF,'matcall',integ,constit);
      pointers(5,:)=MatDes(jMat);
      if ~strcmp(fHandle,'mat_og'); error('Use call');end
    catch; [fHandle,SEopt]=fe_super('call',ElemF,MatDes);
    end
    if isempty(fHandle) cEGI=[];
     disp([ElemF ' : element function/superelement not found and ignored' 7]);
    elseif comstr(fHandle,'mat_og')
    elseif comstr(fHandle,ElemF)
      disp(fHandle);cEGI=[];
    elseif SEopt(1,1)==1; cEGI=EGroup(jGroup); % single superelement
    end
    if i1(1)<=0 cEGI=[];end % do not compute energies
    opt=0; % obsolete kept for compatibility with old elements

  case 'state'

    try;
     fHandle=feval(ElemF,'state');
    catch
     fHandle=''; 
    end
    pointers(5,:)=int32(300);gstate=Case.GroupInfo{jGroup,5};

  case {'rhs','rhs_of','rhs_og'}

    pointers(5,:)=int32(MatDes(jMat));% desired output
    if isfield(EltConst,'RhsDefinition'); fHandle=elem0('rhsmat_og'); 
    else; [fHandle,SEopt]=fe_super('rhscall',ElemF,MatDes);
    end
    if strcmp(fHandle,'rhs_of')
       OutType='rhs_of'; pointers(30,:)=0; % space for face definitions 
       if length(constit)<2; error('Initialization failed');end
    elseif strcmp(fHandle,'rhs_og')
       OutType='rhs_og'; 
    else; OutType='rhs';
    end

  otherwise;  error('Not a supported output');
  end
  if isempty(fHandle) cEGI=[];  % in case the element is not assembled
  else cEGI = EGroup(jGroup)+1:EGroup(jGroup+1)-1;
    NodePos=elt(cEGI,inode)'; 
    if min(min(NodePos))==0; 
     i1=find(NodePos);NodePos(i1)=NNode(NodePos(i1));
    else;NodePos=reshape(NNode(NodePos),length(inode),length(cEGI));
    end
    NodePos=int32(full(NodePos));
  end
  if EGID>0

   switch OutType
   case 'mat' % .m file element matrix assembly

     for jElt=1:length(cEGI)
      i1=NodePos(:,jElt); nodeE=Case.Node(i1(find(i1)),[5:7 1]);
      eval([fHandle ';']);
      if isequal(k1,Inf);break;end
      i1=[Case.DofPerElt(jGroup);0;0;DofPos(:,jElt)];
      of_mk('asmsparse',k,int32(i1),k1(:),fullmap);
     end

   case 'mat_of' % Standard of_mk matrix assembly
     for jElt=1:length(cEGI)
      nodeE=Case.Node(NodePos(:,jElt),[5:7 1]);
      point=pointers(:,jElt);
      ke=of_mk(ElemF,int32(point),integ,constit,nodeE);
      i1=[Case.DofPerElt(jGroup);SymFlag;0;DofPos(:,jElt)];
      of_mk('asmsparse',k,int32(i1),ke,elmap);
     end
   case 'mat_og' % of_mk MatrixIntegration assembly
     pointers=int32(pointers);
     if isfield(EltConst,'VectMap') 
      EltConst.VectMap=int32(EltConst.VectMap);
     end
try
    if sdtdef('diag')>10|~isfield(EltConst,'MatrixIntegrationRule') 
      error('Trying mat file');
    end
    of_mk('matrixintegration',DofPos,NodePos,Case.Node, ...
       pointers,integ,constit,gstate, ...
       elmap,InfoAtNode,EltConst,def.def, ...
        k,int32([Case.DofPerElt(jGroup);SymFlag;0]));
catch
     try;
      if sdtdef('diag')>11|~isfield(EltConst,'MatrixIntegrationRule') 
       error('Trying mat file');
      end
      for jElt=1:length(cEGI)
       ke=of_mk('matrixintegration',jElt,NodePos,Case.Node, ...
       pointers,integ,constit,gstate, ...
       elmap,InfoAtNode,EltConst,def);
       i1=[Case.DofPerElt(jGroup);SymFlag;0;DofPos(:,jElt)];
       of_mk('asmsparse',k,int32(i1),ke,elmap);
      end
     catch % loop on elements is done in ELEM0 for development purposes
      elem0('matrixintegration',DofPos,NodePos,Case.Node, ...
      pointers,integ,constit,gstate, ...
      elmap,InfoAtNode,EltConst,def.def, ...
      k,int32([Case.DofPerElt(jGroup);SymFlag;0]));
     end
end

   case 'ener'
     pointers=int32(pointers);
     if isequal(fHandle,'mat_og'); 
      RunOpt.cEGI=int32(cEGI);
      EltConst.VectMap=int32(EltConst.VectMap);
      try;
       of_mk('matrixintegration',DofPos,NodePos,Case.Node, ...
       pointers,integ,constit,gstate, ...
       elmap,InfoAtNode,EltConst,def.def, ...
       RunOpt,int32([Case.DofPerElt(jGroup);-1;0;]));
      catch;
       elem0('matrixintegration',DofPos,NodePos,Case.Node, ...
       pointers,integ,constit,gstate, ...
       elmap,InfoAtNode,EltConst,def.def, ...
       RunOpt,int32([Case.DofPerElt(jGroup);-1;0;]));
      end
     else % standard element
     for jElt=1:length(cEGI)
      i1=NodePos(:,jElt); nodeE=Case.Node(i1(find(i1)),[5:7 1]);
      eval(fHandle);
      in1=double(Case.GroupInfo{jGroup,1}(:,jElt))+1;
      in2=find(in1); r1=zeros(size(in1));in1=in1(in2);
      k1=k1(in2,in2); 
      for j2=1:size(def.def,2)
         r1(in2)=def.def(in1,j2);
         RunOpt.Ener(cEGI(jElt),j2) = real(r1'*k1*r1);
      end
      end % mat_og or other element
      if cputime-t>15 % display every 15 seconds
         st1 = comstr(st1,-7,sprintf('Done %i elements',cEGI(jElt)));
         t=cputime;
      end
     end

   case 'rhs'
     if isempty(fHandle) % ignore these elements
     else;
      for jElt=1:length(cEGI)
       i1=NodePos(:,jElt); nodeE=Case.Node(i1(find(i1)),[5:7 1]);
       defe=def(:,i1);  point=pointers(:,jElt);
       eval(fHandle)
       in1=double(Case.GroupInfo{jGroup,1}(:,jElt))+1;
       in2=find(in1);in1=in1(in2);
       o1(in1,1)=o1(in1,1)+be(in2);
      end
     end
   case 'rhs_og'
      o1=elem0('rhs_og',DofPos,NodePos,Case.Node,pointers,integ, ...
       constit,gstate,elmap,InfoAtNode,EltConst,def,o1);

   case 'rhs_of'  % standardized of_mk call

     if isstruct(def)
      for jElt=1:length(cEGI)
       i1=NodePos(:,jElt);  nodeE=Case.Node(i1,[5:7 1]);
       point=pointers(:,jElt);
       be=of_mk(ElemF,int32(point),integ,constit,nodeE,gstate, ...
         full(def.def(DofPos(:,jElt)+1,:)),EltConst,InfoAtNode);

       in1=double(Case.GroupInfo{jGroup,1}(:,jElt))+1;
       in2=find(in1);in1=in1(in2);
       o1(in1,1)=o1(in1,1)+be(in2);
      end
     else
      for jElt=1:length(cEGI)
       i1=NodePos(:,jElt);  nodeE=Case.Node(i1,[5:7 1]);
       point=pointers(:,jElt);
       be=of_mk(ElemF,int32(point),integ,constit,nodeE,gstate, ...
         def(:,NodePos(:,jElt)),EltConst,InfoAtNode);

       in1=double(Case.GroupInfo{jGroup,1}(:,jElt))+1;
       in2=find(in1);in1=in1(in2);
       o1(in1,1)=o1(in1,1)+be(in2);
      end
     end

   case 'state' % Modification of element internal state

     for jElt=1:length(cEGI)
      i1=NodePos(:,jElt);
      nodeE=Case.Node(i1,[5:7 1]);
      FieldAtEltDof=FieldAtDof(double(Case.GroupInfo{jGroup,1}(:,jElt))+1);
      point=pointers(:,jElt);
      eval(fHandle)
     end

   otherwise;  error('Not a supported output');
   end
  end % if EGID>0

end % loop on groups

switch OutType
case {'mat','mat_of','mat_og'}
  if size(Case.T,1)==size(k,1)&isempty(strfind(Cam,'not'));
     out=tkt_femknl(Case.T,k); 
  elseif isempty(strfind(Cam,'nop')); out=k+spalloc(size(k,1),size(k,2),0);
  else; out=k;
  end
  if nargout==2; out1=Case;end
case 'ener';  out=RunOpt.Ener;
case {'rhs','rhs_of','rhs_og'}
   out=o1;
case 'state';
otherwise error('Not a known OutType output');
end


% -----------------------------------------------------------------------
else sdtw('''%s'' is unknown',CAM);end

% ---------------------------------------------------------------------------
function K=tkt_femknl(T,K)

if sp_util('issdt') K=feutilb('tkt',T,K);
else
 Tt=T';
 K=Tt*K*T;
end
