function [o1,o2,o3]=fe_stres(varargin);

%FE_STRES compute stress or energies
%
%       Syntax: result = fe_stres('Command',MODEL,DEF)
%       
%	MODEL can be specified as
%         - four arguments NODE, ELT, PL, IL (see : doc node ... and FE_MK)
%         - a structure with fields .Node, .Elt, .pl, .il 
%         - a database wrapper with those fields
%       DEF (the deformations) can be specified using
%         - two arguments MODE (deformations) and MDOF (DOF definition vector)
%         - a structure with fields .def, .DOF
%
%	Accepted commands are linked to the following calls
%	
%	- RESULT = fe_stres('ener (element selection)',MODEL,DEF)
%	  where the strain and kinetic energies in the elements selected with
%         (element selection) (see the FEMESH FINDELT command for details)
%         are computed for the deformations in DEF and the result is returned
%         in the structure RESULT with fields .StrainE, .KinE and .IndInElt
%         which specifies which elements were selected.
%
%         For backward compatibility, if no element selection is given
%         FE_STRES returns [StrainE,KinE] as two arguments.
%         To select all elements, use the 'ener groupall' command.
%
%	- RESULT=fe_stres('stress (CritFcn) (rest)',MODEL,DEF,EltSel) 
%         Supported criteria are (CritFcn)
%          sI, sII, sIII principal stresses from max to min
%         Supported restitution are
%          1 (sum at nodes wheighted by number of connected elements)
%        hexa8('test eig stress') gives an example
%
%       - def=fe_stres('expand',model,Case) extrapolates stresses from
%        integration point to nodes
%
%	See also help feplot, fecom, 
%                d_ubeam (demo)

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


% check of inputs

% CAM,node,elt,pl,il,tr,mdof
CAM=varargin{1}; [CAM,Cam]=comstr(CAM,1); st1=''; carg=2;


% check the model input
if carg>nargin & strcmp(Cam,'cvs')
 o1='$Revision: 1.57 $  $Date: 2006/05/12 14:33:02 $';return;
else; model=varargin{carg}; carg=carg+1;
end
if isa(model,'double') % all arguments given
  if size(model,2)~=7 error('improper node definition'); end
  model=struct('Node',model,'Elt',varargin{carg}, ...
  'pl',varargin{carg+1},'il',varargin{carg+2});
  carg=carg+3;

elseif isstruct(model)
  if ~isfield(model,'Node')|~isfield(model,'Elt')
    error('Not a proper model definition');
  end
  model.pl=fe_mat('getpl',model); model.il=fe_mat('getil',model);
elseif isa(model,'sdth')
  try 
    model=struct('Node',model.Node,'Elt',model.Elt,'pl',model.pl, ...
     'il',model.il);
  catch
    error('Not a proper model definition'); 
  end
end
NNode(model.Node(:,1))=1:size(model.Node,1);

% check the deformation input
if carg<=nargin 
 def=varargin{carg}; carg=carg+1;
 if isa(def,'double') % all arguments given
   def=struct('def',def,'DOF',varargin{carg}); carg=carg+1;
 elseif isstruct(def)
  % do some checks
 elseif isa(tr,'sdth'); error('Not a supported case'); end
 if isfield(def,'GroupInfo') % really an expand call
  Case=def;
 else
  if isempty(def.DOF) error('No DOF definition vector');end
  if isempty(def.def) error('No deformation given');end
  if size(def.def,1)~=size(def.DOF,1)
   error('DEF.DEF and DEF.DOF are not compatible');
  end
 end
end

% ----------------------------------------------------------------------------
if comstr(Cam,'en') % 'ener'

[CAM,Cam]=comstr(CAM,'energy','%c');

if isfield(model,'Opt')&model.Opt(1)==3  % Catch SDT type 3 superelement
 if carg<=nargin&ischar(varargin{carg});
   eval(sprintf( ...
   'o1=upcom(model,''ener%s'',''%s'',def);',CAM,varargin{carg}),'o1=[];');
 else eval(sprintf('o1=upcom(model,''ener%s'',def);',CAM),'o1=[];');end
 return;
end

if ~isempty(strfind(Cam,'dens'))
  RunOpt.dens=1; CAM(strfind(Cam,'dens')+[0:3])='';[CAM,Cam]=comstr(CAM,1);
else; RunOpt.dens=0; 
end

if isempty(Cam); RunOpt.MatDes=1;
elseif Cam(1)=='k'; RunOpt.MatDes=1; [CAM,Cam]=comstr(CAM,2);
elseif Cam(1)=='m'; RunOpt.MatDes=2; [CAM,Cam]=comstr(CAM,2);
else; RunOpt.MatDes=1;
end

if isempty(Cam)&carg<=nargin&ischar(varargin{carg}); % selection in other arg
 CAM=varargin{carg};carg=carg+1;[CAM,Cam]=comstr(CAM,1);
end
if ~isempty(Cam)
 el0=model.Elt;[ind,model.Elt,ind2] = feutil(['findelt' CAM],model);
end
[eltid,model.Elt]=feutil('eltidfix',model.Elt);
[EGroup,nGroup]=getegroup(model.Elt);
if nGroup == 0		error('no element group specified in ELT'); end
opt=0;

% Verify that DOF numbering is consistent with mdof
[Case,model.DOF]=fe_mknl('initnocon NoT',model);
if isequal(def.DOF,model.DOF)
elseif sp_util('issdt'); def=feutil('placeindof',model.DOF,def);
else; def.def=fe_c(model.DOF,def.DOF,def.def')';def.DOF=model.DOF;
end
if RunOpt.dens; % energy density
  o1=fe_mknl('assemble',model,Case,250+RunOpt.MatDes,def);
  def.def=zeros(size(def.def,1),1);
  def.def(fe_c(def.DOF,.01,'ind'),1)=1; % RB mode for density
  o2=fe_mknl('assemble',model,Case,252,def);
  r2=o2(:,end);r2=r2/sum(r2); ind=find(r2);
  o1(ind,:)=o1(ind,:)./r2(ind,ones(size(o1,2),1));
else
  o1=fe_mknl('assemble',model,Case,250+RunOpt.MatDes,def);
end
if isreal(def);% modified 12/05/06  the 1/2 factor was in doc but not software
     o1=o1/2;
else; o1=o1/4; % if complex integration over a cycle
end  
if any(any(o1<0)); sdtw('_nb','Some energies are negative.');end

ind2=find(eltid);if max(ind2)>size(o1,1);ind2=ind2(find(ind2<size(o1,1)));end
o1=struct('data',o1(ind2,:),'EltId',eltid(ind2));
switch RunOpt.MatDes
case 1; o1.name='Ener k'; case 2;o1.name='Ener m'; 
otherwise o1.name=sprintf('Ener %i',RunOpt.MatDes');
end
if RunOpt.dens; o1.name=[o1.name '_dens'];end

% ----------------------------------------------------------------------------
% Extrapolation stress field at integration points to stress at nodes
% def=fe_stres('expand',model,Case)
elseif comstr(Cam,'expand')

[EGroup,nGroup]=getegroup(model.Elt);
mdof=feutil('findnode groupall',model);
NNode=sparse(mdof+1,1,1:length(mdof));
conn=int32(zeros(size(model.Elt,2),size(model.Elt,1)));

% Matrix Init
for jGroup=1:nGroup;
 [ElemF,opt,ElemP]= feutil('getelemf',model.Elt(EGroup(jGroup),:),jGroup);
 if any(strcmp(ElemF,{'celas','mass1','mass2'})); continue;end
 cEGI = EGroup(jGroup)+1:EGroup(jGroup+1)-1;
 i1=fe_super('node',ElemF);
 Case.DofPerElt(jGroup)=length(i1);
 i2=model.Elt(cEGI,i1);
 i2=int32(reshape(full(NNode(i2+1)-1),size(i2,1),size(i2,2))');
 Case.GroupInfo{jGroup,1}=i2;  % DofPos
 pointers=Case.GroupInfo{jGroup,2};  % Pointers
 pointers(1,:)=length(i1)^2;  % OutSize 1
 pointers(5,:)=1;             % desired matrix
 pointers(6:7,:)=0;           % integ constit offset
 Case.GroupInfo{jGroup,2}=int32(pointers);  % Pointers
  
 Case.GroupInfo{jGroup,3}=int32([1;1;length(i1)]);  % integ
 Case.GroupInfo{jGroup,4}=1;  % Constit
 Case.GroupInfo{jGroup,6}=reshape(1:length(i1)^2,length(i1),length(i1));
    % ElMap
 
 conn(1:size(i2,1),cEGI)=i2;

 opt=Case.GroupInfo{jGroup,8};  %EltConst
 opt.StrainDefinition={[1 1 1 1 length(i1)]};
 opt.StrainLabels={{'s'}}; opt.ConstitTopology={1};
 opt=integrules('matrixrule',opt);
 Case.GroupInfo{jGroup,8}=opt;
 
end

 Case.MatGraph=of_mk('meshGraph', ...
	      size(mdof,1),nGroup,size(conn,1), ...
	      int32(EGroup),int32(diff(EGroup)-1),int32(Case.DofPerElt), ...
	      conn );
 k=Case.MatGraph; of_mk( 'fillvalue', k, 0 );

% Loop to assemble pseudo-mass and associated RHS
elt=model.Elt;NNode=sparse(model.Node(:,1),1,1:size(model.Node,1));
b=zeros(size(k,1),size(Case.GroupInfo{jGroup,5},3));
o1=struct('def',[],'DOF',mdof+.99,'label','Strain extended at nodes', ...
  'lab',{{}});

for jGroup=1:nGroup;

  [ElemF,i1]= feutil('getelemf',model.Elt(EGroup(jGroup),:),jGroup);
  if any(strcmp(ElemF,{'celas','mass1','mass2'})); continue;end
  EGID=i1(1); 
  inode=fe_super('node',ElemF);
  cEGI = EGroup(jGroup)+1:EGroup(jGroup+1)-1;

  % 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};

  NodePos=reshape(NNode(elt(cEGI,inode)'),length(inode),length(cEGI));
  NodePos=int32(full(NodePos));
  cEGI = EGroup(jGroup)+1:EGroup(jGroup+1)-1;
  of_mk('matrixintegration',DofPos,NodePos,Case.Node, ...
       pointers,integ,constit,gstate, ...
       elmap,InfoAtNode,EltConst,[],k,int32([Case.DofPerElt(jGroup);0;0]));
  Nw=EltConst.Nw;
  for jElt=1:length(cEGI)
      % Two ways to store stresses : (jW,jElt,jType) or (jType*jW,jElt)
      nodeE=model.Node(NodePos(:,jElt),5:7);
      of_mk('buildndn',pointers(4,jElt),EltConst,nodeE);
      if size(gstate,3)>1
       be=(squeeze(gstate(:,jElt,:))'*diag(EltConst.jdet)*EltConst.N)';
      else
       be=(reshape(gstate(:,jElt),size(gstate,1)/Nw,Nw) ...
         *diag(EltConst.jdet)*EltConst.N)';
      end
      if isempty(o1.lab);o1.lab={}; for j2=1:size(be,2);
        o1.lab{j2,1}=sprintf('Field %i',j2);
      end;end
      i2=double(DofPos(:,jElt))+1;
      if size(b,2)<size(be,2);b(1,size(be,2))=0;end
      b(i2,:)=b(i2,:)+be;
  end

end
kd=ofact(k); o1.def=kd\b; ofact('clear',kd);


% ----------------------------------------------------------------------------
elseif comstr(Cam,'thermal') 
% Computation of thermal induced stresses in volume elements of the *b family


[Case,model.DOF]=fe_mknl('init',model);

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

% build a defT defined at nodes
defT=zeros(size(Case.Node,1),1);
ind=fe_c(def.DOF,.20,'ind');
defT(NNode(fix(def.DOF(ind))),1)=def.def(ind,1);
elt=model.Elt;node=Case.Node;
ci_ts_eg=[1 5 9 8 7 4]';

for jGroup=1:nGroup  % Loop on groups

  % See documentation elem0#call
  DofPos=Case.GroupInfo{jGroup,1};
  pointers=Case.GroupInfo{jGroup,2};
  integ=Case.GroupInfo{jGroup,3};
  constit=Case.GroupInfo{jGroup,4};
  elmap=int32(Case.GroupInfo{jGroup,6});
  InfoAtNode=Case.GroupInfo{jGroup,7};
  EltConst=Case.GroupInfo{jGroup,8};

  if ~isfield(EltConst,'Nw'); continue;end
  gstate=zeros(6*EltConst.Nw,size(Case.GroupInfo{1,1},2));
  cEGI = EGroup(jGroup)+1:EGroup(jGroup+1)-1;
  NodePos=reshape(NNode(elt(cEGI,1:EltConst.Nnode)'), ...
   EltConst.Nnode,length(cEGI));

  if size(constit,1)<48; error('Thermal properties not defined');end
  idd=find(EltConst.ConstitTopology{1});

  for jElt=1:size(DofPos,2) %
   point=pointers(:,jElt);
   DD=double(EltConst.ConstitTopology{1});
   DD(idd)=constit(DD(idd)+point(7));
   at=reshape(constit([39:47]+point(7)),3,3);

   nodeE=node(NodePos(:,jElt),5:7);

   if (point(4)~=2 & point(4)~=3 & point(4)~=23) 
    error('Not a supported NDN strategy'); 
   end
   of_mk('buildndn',point(4),EltConst,nodeE);

   Nstrain=size(DD,1);

   for jW=0:EltConst.Nw-1
    % Temperature at integration point
    r1=(defT(NodePos(:,jElt))-constit(48+point(7)))'*EltConst.NDN(:,jW+1);
    % Stress at integration point
    r1=at*r1;r1=DD*r1(ci_ts_eg);
    gstate(jW*Nstrain+[1:Nstrain],jElt)=r1;
   end % loop on integration points
  end % loop on elements
  Case.GroupInfo{jGroup,5}=gstate;

end

o1=Case;


% ----------------------------------------------------------------------------
elseif comstr(Cam,'st') 
[CAM,Cam]=comstr(CAM,'stress','%c');[CAM,Cam]=comstr(CAM,1);

CritFcn=stack_get(model,'info','StressCritFcn','getdata');
if isempty(CritFcn);
 CritFcn='r1=Principal(1,r1,[],TensorTopology);'; DIRS=[];
end

% Attempt to deal with where the stress is evaluated.
RunOpt=struct('typ',1,'header','At node'); % sum on nodes
i1=strfind(Cam,'atinteg'); if ~isempty(i1);
 RunOpt.typ=3;CAM(i1+[0:6])='';RunOpt.header='At integ';
end
i1=strfind(Cam,'atcenter'); if ~isempty(i1);
 RunOpt.typ=2;CAM(i1+[0:7])='';RunOpt.header='At center';
end
i1=strfind(Cam,'atnode'); if ~isempty(i1)RunOpt.typ=1;CAM(i1+[0:5])='';end
i1=strfind(Cam,'gstate'); 
if ~isempty(i1) RunOpt.typ=3;CAM(i1+[0:5])='';CritFcn='';end
[CAM,Cam]=comstr(CAM,1);

if RunOpt.typ==3
elseif comstr(Cam,'siii');    CritFcn='r1=Principal(3,r1,[],TensorTopology);'; 
elseif comstr(Cam,'sii');    CritFcn='r1=Principal(2,r1,[],TensorTopology);'; 
elseif comstr(Cam,'sidir');
  CritFcn='[r1,dir]=Principal(21,r1,bas,TensorTopology);DIRS(cEGI(jElt),1:length(dir),jDef)=dir;'; bas=eye(3);
elseif comstr(Cam,'si');    CritFcn='r1=Principal(1,r1,[],TensorTopology);'; 
elseif comstr(Cam,'mises'); CritFcn='r1=Principal(4,r1,[],TensorTopology);'; 
elseif comstr(Cam,'sxx');   CritFcn='r1=Principal(11,r1,[],TensorTopology);'; 
elseif comstr(Cam,'syy');   CritFcn='r1=Principal(12,r1,[],TensorTopology);'; 
elseif comstr(Cam,'szz');   CritFcn='r1=Principal(13,r1,[],TensorTopology);'; 
elseif comstr(Cam,'none');   CritFcn='';
else; sdtw('_nb','Using max principal stress');
end
if ischar(CritFcn)&~isempty(CritFcn);
  RunOpt.header(2,1:length(CritFcn))=CritFcn;
end

if carg<=nargin & ischar(varargin{carg})
 [eltid,model.Elt]=feutil('eltidfix',model.Elt);
 [ind,elt] = feutil('findelt',model,varargin{carg});carg=carg+1;
else elt=model.Elt; end

%mods=model;mods.Elt=elt;Case=fe_mknl('initnocon',mods);

[EGroup,nGroup]=getegroup(elt);
if nGroup == 0;	error('no element group specified in ELT'); end

NNode=sparse(model.Node(:,1),1,1:size(model.Node,1));
mode=zeros(size(model.Node,1),3);


model.Elt=elt;
% Verify that DOF numbering is consistent with mdof
[Case,mdof]=fe_mknl('initnocon',model);node=Case.Node;
if isequal(def.DOF,mdof)
elseif sp_util('issdt'); def=feutil('placeindof',mdof,def);
else; def.def=fe_c(mdof,def.DOF,def.def')';def.DOF=mdof;
end

switch RunOpt.typ(1)
case 1 % nodal stress field
 o1=zeros(size(model.Node,1),size(def.def,2)); o2=zeros(size(model.Node,1),1);
case 2 % element stress value
 o1=zeros(size(elt,1),size(def.def,2)); o2=[];
case 3 % Stress at integration points
otherwise error('Not a valid stress option');
end

% loop over element groups - - - - - - - - - - - - - - - - - - - - - -
for jGroup = 1:nGroup

[ElemF,gOpt,ElemP]= feutil('getelemf',elt(EGroup(jGroup),:),jGroup);

EltConst=Case.GroupInfo{jGroup,8};
if size(Case.GroupInfo,2)>7 & ...
  isfield(EltConst,'MatrixIntegrationRule')& ...
  isfield(EltConst,'StrainDefinition')
  try;
   if RunOpt.typ==1;   ;EltConst=integrules('stressrule',EltConst,'node');
   elseif RunOpt.typ==2;EltConst=integrules('stressrule',EltConst,'center');
   else; EltConst=integrules('stressrule',EltConst);
   end
   ElemP='mat_og';
  catch;  
   if sdtdef('diag')>=10; error(lasterr);
   else;disp(lasterr);sdtw('_nb','Stress rule problem');
   end
  end
end

cEGI = EGroup(jGroup)+1:EGroup(jGroup+1)-1;
inode=fe_super('node',ElemF);
    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
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};

if ~isempty(integ)&~isa(integ,'int32') error('Improper class');end

prop=fe_super('prop',ElemF);i1=fe_super('node',ElemF);

try;
if gOpt(1)<0 % not a display element
elseif strcmp(ElemP,'mat_og')

   [stress,Nw]=elem0('stress_og',NodePos,jGroup,Case,def, ...
    pointers,integ,constit,EltConst); 
   Nelt=size(NodePos,2); Ns=size(stress,1)/Nw;
   if isfield(EltConst,'TensorTopology');TensorTopology=EltConst.TensorTopology;
   elseif Ns==6; TensorTopology=[1 6 5;6 2 4;5 4 3];
   else;  TensorTopology=[1 3;3 2];
   end

   if RunOpt.typ==3&isempty(CritFcn); % gstate
     Case.GroupInfo{jGroup,5}=stress;
   else % do some post-processing on elements
     for jElt=1:Nelt; for jDef=1:size(def.def,2);
      r1=reshape(stress(:,jElt+Nelt*(jDef-1)),Ns,Nw);
      eval(CritFcn); 
      if RunOpt.typ~=3;
      [o1,o2]=post_expand(o1,o2,NodePos(:,jElt),jDef,cEGI(jElt),r1,RunOpt.typ);
      else; Case.GroupInfo{jGroup,5}(:,jElt,jDef)=r1;
      end
     end;end
   end % CritFcn or Not
else

switch ElemF

% 3-D elements - - - - - - - - - - - - - - - - - - - - - - -
case {'hexa8','penta6','hexa20','tetra4','tetra10','penta15'} % 3-D elements

 % hexa8 : iopt=[ndof SizeOfOut 200 Unused telt tcar noref]
 %         for car format is fe_mat('of_mk')
 iopt=[length(fe_super('dof',ElemF)) length(i1)*6 200 0 0 2];
 TensorTopology=[1 4 6;4 2 5;6 5 3];
 for jElt=1:length(cEGI)
    nodeE=Case.Node(NodePos(:,jElt),[5:7 1]); i2=NodePos(:,jElt);
    %i2=NNode(elt(cEGI(jElt),i1));nodeE=node(i2,[5:7 1]);
    
    point=pointers(:,jElt); point(5)=200;
    point(1)=double(integ(double(point(6))+4))*6;
    state=zeros(3*size(NodePos,1),1);
    i3=double(Case.GroupInfo{jGroup,1}(:,jElt))+1;i4=find(i3);i3=i3(i4);
    for jDef=1:size(def.def,2)
     state(i4)=def.def(i3,jDef);
     bi=of_mk(ElemF,int32(point),integ,constit,nodeE,[],state);
     r1=reshape(bi,6,length(i1));
     eval(CritFcn);
     [o1,o2]=post_expand(o1,o2,i2,jDef,cEGI(jElt),r1,RunOpt.typ);
    end
 end

% 2-D elements - - - - - - - - - - - - - - - - - - - - - - -
case {'t3p','q4p','t6p','q5p','q8p'}

 % t3p : iopt=[ndof SizeOfOut 200 Unused telt tcar noref]
 %         for car format is fe_mat('of_mk') 
 % 3 for sigma xx yy xy
 
 for jElt=1:length(cEGI)
    nodeE=Case.Node(NodePos(:,jElt),[5:7 1]); i2=NodePos(:,jElt);
    point=pointers(:,jElt); point(5)=200;
    % axi : S_rr, S_rz, S_zz, S_tt
    % 2d  : S_xx, S_yy, S_xy 
    % r1=reshape(bi,3,1);
    if   integ(double(point(6))+7)==3
     point(1)=4;TensorTopology=[1 2 5;2 3 5;5 5 4]; %axi S_rr, S_rz, S_zz, S_tt
    else point(1)=3; TensorTopology=[1 3;3 2];% plane Sxx yy xy
    end 

    state=zeros(2*size(NodePos,1),1);
    i3=double(Case.GroupInfo{jGroup,1}(:,jElt))+1;i4=find(i3);i3=i3(i4);
    for jDef=1:size(def.def,2)
     state(i4)=def.def(i4,jDef);
     bi=of_mk(ElemF,int32(point),integ,constit,nodeE, ...
        [0 0 0 ...           % alpha
        zeros(1,10)],state, ... % deformation at element
        zeros(10));  % Ti
     r1=bi(:); eval(CritFcn);
     [o1,o2]=post_expand(o1,o2,i2,jDef,cEGI(jElt),r1,RunOpt.typ);
    end
 end
% 2-D plate elements - - - - - - - - - - - - - - - - - - - - - - -
case {'dktp'}


for jDef=1:size(def.def,2) 
  r1=feutil('dof2mode',full(def.def(:,jDef)),def.DOF);
  mode(NNode(r1(:,1)),1:size(r1,2)-1)=r1(:,2:end);
  for jElt=1:length(cEGI)
    i2=NNode(elt(cEGI(jElt),i1));nodeE=node(i2,[5:7 1]);
    point=pointers(:,jElt); point(1)=3; point(5)=200;

    bi=of_mk(ElemF,int32(point),integ,constit,nodeE, ...
        [0],mode(i2,[3 4 5 1])',[0],[0]); % deformation at element
    r1=bi(:);
    eval(CritFcn);
    [o1,o2]=post_expand(o1,o2,i2,jDef,cEGI(jElt),r1,RunOpt.typ);
  end
end


case {'quad4'}  % this IS WRONG AND FOR TESTING ONLY

TensorTopology=[1 3;3 2];

for jDef=1:size(def.def,2) 
  r1=feutil('dof2mode',full(def.def(:,jDef)),def.DOF);
  mode(NNode(r1(:,1)),1:size(r1,2)-1)=r1(:,2:end);
  pointers(5,:)=210;
  for jElt=1:length(cEGI)
    nodeE=Case.Node(NodePos(:,jElt),[5:7 1]); i2=NodePos(:,jElt);
    [bi,bas]=quad4(nodeE,elt(jElt,:),pointers(:,jElt),integ,constit, ...
        [],mode(i2,1:3)');

    % 2d  : S_xx, S_yy, S_xy 
    % r1=reshape(bi,3,1);
    r1=bi(:);
    eval(CritFcn);
    [o1,o2]=post_expand(o1,o2,i2,jDef,cEGI(jElt),r1,RunOpt.typ);
  end
end

otherwise
 try;
   eCall=feval(ElemF,'call');
 catch
  sdtw('_nb','stress computation not supported %s',ElemF);break;
 end
end; end % display or not
catch; % elements
  sdtw('_nb','stress computation failed for group %i(%s)',jGroup,ElemF);
end

end % loop on element groups

switch RunOpt.typ(1)  % Post process
case 1         % at node
 i1=find(o2); o1(i1,:)=o1(i1,:)./o2(i1,ones(size(o1,2),1))*size(o1,2);
 o1=struct('DOF',model.Node(:,1)+.99,'data',o1,'header',RunOpt.header);
case 2         % at center
 i1=feutil('eltidfix',elt);
 i2=find(i1);

 o1=struct('EltId',i1(i2),'data',o1(i2,:),'header',RunOpt.header);
 if ~isempty(DIRS) o1.DIRS=DIRS(i2,:,:); end
case 3; o1=Case;
end

% ---------------------------------------------------------------------------
else error([CAM ' : not a available command']);
end

% -------------------------------------------------------------------------
function [out,dir] = Principal(ind,r1,bas,TensorTopology);

r2=zeros(size(TensorTopology));
if size(r1,1)==6 % 3-D stress

  % Sxx Syy Szz Sxy Syz Szx
  for j1=1:size(r1,2)
    r2(1:9)=r1(TensorTopology,j1);s=sort(eig(r2)); 
    switch ind
    case 1; out(j1)=s(1);
    case 2; out(j1)=s(2);
    case 3; out(j1)=s(3);
    case 4; out(j1)=sqrt(((s(1)-s(2))^2+(s(2)-s(3))^2+(s(3)-s(1))^2)/2);
    case 11; out(j1)=r1(1,j1);
    case 12; out(j1)=r1(2,j1);
    case 13; out(j1)=r1(3,j1);
    end
  end

elseif size(r1,1)==4 % 2-D axi stress

  r1(5)=0;% S_rr, S_rz, S_zz, S_tt
  for j1=1:size(r1,2)
    r2(1:9)=r1(TensorTopology,j1);s=sort(eig(r2)); 
    switch ind 
    case 4;   out(j1)=sqrt(((s(1)-s(2))^2+(s(2)-s(3))^2+(s(3)-s(1))^2)/2);
    case 11; out(j1)=r1(1,j1);
    case 12; out(j1)=r1(2,j1);
    case 13; out(j1)=r1(3,j1);
    case 21;  
     [dir,s]=eig(r2);[s,i2]=max(abs(diag(s)));
     dir=dir(:,i2);out(j1)=s(end);
     dir=bas*[dir;0];
    otherwise; out(j1)=s(ind); 
    end
  end


elseif size(r1,1)==3 % 2-D stress

  if any(ind==3) error('Can''t compute sIII for 2-D stress');end
  % Sxx yy xy
  for j1=1:size(r1,2)
    r2(1:4)=r1(TensorTopology,j1);s=sort(eig(r2)); 
    switch ind 
    case 4; 
     s(3)=0;out(j1)=sqrt(((s(1)-s(2))^2+(s(2)-s(3))^2+(s(3)-s(1))^2)/2);
    case 11; out(j1)=r1(1,j1);
    case 12; out(j1)=r1(2,j1);
    case 13; out(j1)=r1(3,j1);
    case 21;  
     [dir,s]=eig(r2);[s,i2]=max(abs(diag(s)));
     dir=dir(:,i2(end));out(j1)=s(end);
     dir=bas*[dir;0];
    otherwise; out(j1)=s(ind); 
    end
  end
end

% -------------------------------------------------------------------------
function   [o1,o2]=post_expand(o1,o2,i2,jDef,jElt,r1,opt);

    switch opt(1)
    case 1     % reexpand to nodes
      o1(i2,jDef)=o1(i2,jDef)+r1(:);
      o2(i2, 1)=o2(i2, 1)+1;
    case 2     % Give at element center
      o1(jElt,jDef)=o1(jElt,jDef)+mean(r1);
    end

% -------------------------------------------------------------------------
% Computation of stress state associated with a deformation mode
function [out,nw]=rule_in_mfile(NodePos,jGroup,Case,def,pointers,integ,constit,EltConst)
%   ke=rule_in_mfile(NodePos,jGroup,Case,def,pointers,integ,constit,EltConst);

