function [out,out1,out2,out3,out4]=feutil(varargin);

%FEUTIL Finite element model handling utilities
%
%	This function groups a number of mesh handling utilities with all
%	arguments given (unlike feplot, femesh, fecom which use context 
%	information).
%	sdtdef('epsl') gives the mesh precision (smallest distance considered)
%	whose default is 1e-6.
%
%	Accepted calls are (see the manual for details)
%
%       [AllNode,ind]=femesh('AddNode',OldNode,NewNode);
%       Pairs = feutil('closestnode',node1,node2) -> closest node matching
%       Mode      = feutil('dof2mode',def,dof) -> dof to [node disp] form
%       [def,dof] = feutil('mode2dof',mode)    -> [node disp] to dof form
%       [model,def]=feutil('DefLocalToGlocal',model,def) deformation in global
%       elt=feutil('divide group 1 withnode{x>10}',model) ->
%                                                    divide group 1 in 2 groups
%       [node,elt]=feutil('divide HexaTransition',node,elt,nodeID);
%       [EGID,eltnames]=feutil('egid',elt) or 
%            elt=feutil('egid',elt,EGID,eltnames)
%       [EltId]  = feutil('eltid [,fix]',elt)
%       [EltInd] = feutil('eltind',elt)
%       [ind,snode] = feutil('findnode (selection)',node,elt,elt0, ...)
%	  See the femesh findnode command for more details
%       [ind,els,indwithead] = feutil('findelt (selection)',node,elt,el0, ...)
%	  See the femesh findelt command for more details
%
%	[Line,nInf,node_ind]  = feutil('get Line',node,elt,nInf)
%	  node empty => use indices
%	Patch = feutil('get Patch',node,elt) : node empty => use node numbers
%	[ElemF,opt,ElemP]= feutil('getelemf',elt(EGroup(jGroup),:),[jGroup])
%	[mdof,tr,nword]  = feutil('get DOF',elt,adof,tr)
%	[nor,cg]  = feutil('get normal[elt,node]',node,elt)
%	Map  = feutil('get normalMap[elt,node]',node,elt)
%
%	feutil('info elt Name',elt) information on model description matrix
%       elt=feutil('joingroup1:2',model.Elt)
%	[el_match,point,el0_normal]=feutil('match',node,elt,el0) find faces 
%         in elt matching faces in el0.
%	[out]=feutil('[MatId ProId MPId]',elt) mat/pro ID for each element
%	elt=feutil('mpid',elt,mpid) set mat/pro ID for each element
%       Object : object constructors ( see more details in doc('femesh') )
%          HoleInPlate, BeamLine, Mass, Beam, Quad, Hexa, Circle, Cylinder
%       [node,elt]=feutil('optimmodel',model) removes unused nodes from FEnode
%       model=feutil('optimeltcheck',model) : ill conditionning fixes
%       [node,elt]=feutil('optimnodenum',model) optimizes numbering
%       elt=feutil('orient',model) reorient elements if necessary
%       elt=feutil('removeelt (selection)',model) 
%	st=feutil('stringdof',mdof) string labels for DOFs
%       elt=feutil('trace2elt',ldraw) trace line for elt format
%
%	See also help feplot
%                doc  femesh, findnode, findelt, fem

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


try;epsl=sp_util('epsl');catch;epsl=sdtdef('epsl'); end
ModelStack={};
[CAM,Cam]=comstr(varargin{1},1);carg=2;

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

% ('addnode',OldNode,NewNode) - - - - - - - - - - - - - - - - - - -
if comstr(Cam,'node'); [CAM,Cam] = comstr(CAM,5);

 r1=varargin{carg};carg=carg+1;
 r2=varargin{carg};carg=carg+1; 

 if ~isempty(r1) & any(~finite(r1(end,:))) r1 = r1(1:end-1,:); end
 if ~isempty(r2) & any(~finite(r2(end,:))) r2 = r2(1:end-1,:); end
 if size(r2,2)==3 r2=[[1:size(r2,1)]' zeros(size(r2,1),3) r2];end

 out=r1;  out1 = zeros(size(r2,1),1);

 [epsl,CAM,Cam]=test_epsl(epsl,CAM,Cam);
 if comstr(Cam,'from')   [opt,CAM,Cam]=comstr(Cam,'from','%i'); opt=opt-1;
 elseif ~isempty(out) opt=max(out(:,1)); else opt=0; end

 if ~isempty(strfind(Cam,'new'))
  opt=comstr(CAM(strfind(Cam,'new')+4:end),[-1 opt]);
  out1=size(out,1)+[1:size(r2,1)]';
  out(end+[1:size(r2,1)],1:7)=[opt+[1:size(r2,1)]' r2(:,2:7)];

 elseif size(r2,1)>10 % try a new strategy based on sortrows

  % catch for SDT strategy to do the same thing
  if sp_util>5.0007;
   try;eval('[out,out1]=feutilb(''addnode'',r1,r2,epsl,opt);');return;end
  end
  % i3 contains a list corresponding to sort of [r1;r2]
  if ~isempty(r1); r3=[r1(:,5:7);r2(:,5:7)];
  else; r3=[r2(:,5:7)];
  end
  [u,s,v]=svd(r3,0);r3=r3*v; % orient along main mesh directions
  [r4,i3]=sort(r3(:,1));r3=r3(i3,:);ind=i3; 

  inx=find(abs(diff(r3(:,1)))<epsl); % x is matching
  while ~isempty(inx)
   % one segment of matching x
   if length(inx)==1 i2=1; else i2=min(find(diff(inx)>1));end
   if isempty(i2) i2=length(inx);end
   iny=inx(1):inx(i2)+1;
   [r4,i4]=sort(r3(iny,2));r3(iny,:)=r3(iny(i4),:); 
   i3(iny)=i3(iny(i4)); ind(iny)=ind(iny(i4));
   iny=iny(find(abs(diff(r4))<epsl));

   while ~isempty(iny) % y is matching
    if length(iny)==1 i5=1; else i5=min(find(diff(iny)>1));end
    if isempty(i5) i5=length(iny);end
    inz=iny(1):iny(i5)+1;
    [r4,i4]=sort(r3(inz,3)); r3(inz,:)=r3(inz(i4),:); 
    i3(inz)=i3(inz(i4));ind(inz)=ind(inz(i4));
    inz=inz(find(abs(diff(r4))<epsl));

    while ~isempty(inz) % z is matching
     if length(inz)>1 i6=min(find(diff(inz)>1));else i6=1;end
     if isempty(i6) i6=length(inz);end
     i7=inz(1):inz(i6)+1;
     [r4,i4]=min(i3(i7)); % verify coincidence
     r4=sqrt(sum((r3(i7,:)-ones(size(i7(:)))*r3(i7(i4),:)).^2,2))<epsl;
     for j1=find(r4(:)'==0&i3(i7)'>size(r1,1)) % when failed do a direct search
       r5=sqrt(sum((r3(i7,:)-ones(size(i7(:)))*r3(i7(j1),:)).^2,2));
       i5=find(r5<epsl);i3(i7(j1))=min(i3(i7(i5)));
     end % when failed do a direct search
     i7=i7(find(r4));  i3(i7)=min(i3(i7));
     inz=inz(i6+1:end);
    end

    iny=iny(i5+1:end);
   end % y is matching
   inx=inx(i2+1:end);
  end % x is matching

  % check coincident nodes
  i4=find(i3<size(r1,1));  

  i4=find(i3>size(r1,1)); if ~isempty(i4) % 
    i4=unique(i3(i4));
    out=[r1;opt+[1:length(i4)]' r2(i4-size(r1,1),2:end)];
    % new indices of r2 in out
    nind=[[1:size(r1,1)]';i4];nind(nind)=1:length(nind);
    out1(ind)=nind(i3); out1=out1(size(r1,1)+1:end);
  else
    out=r1;
    nind=[1:size(r1,1)];nind(nind)=1:length(nind);
    out1(ind)=nind(i3);
    out1=out1(size(r1,1)+1:end);
  end
  
 else % old strategy that scans through each node
  opt=opt(1)+1;
  for j1=1:size(r2,1)
    if isempty(out) out1(j1)=1;out=[opt r2(j1,2:7)];opt=opt+1;
    else
      i1 = find(abs(out(:,5)-r2(j1,5))<epsl);
      if ~isempty(i1)   i1 = i1(find(abs(out(i1,6)-r2(j1,6))<epsl)); end
      if ~isempty(i1)   i1 = i1(find(abs(out(i1,7)-r2(j1,7))<epsl)); end
      if ~isempty(i1)   out1(j1) = i1(1);
      else
        out = [out;opt(1) r2(j1,2:7)];opt=opt+1; out1(j1,1)=size(out,1);
      end
    end % of ~isempty(out)
  end
 end % 'new or not'

% ('addtest',model,model2) - - - - - - - - - - - - - - - - - - -
elseif comstr(Cam,'test'); [CAM,Cam] = comstr(CAM,5);

model=varargin{carg};carg=carg+1;
el0=varargin{carg};carg=carg+1; 
i1=strfind(Cam,'keeptest'); if ~isempty(i1)
 RunOpt.KeepTest=1;
else;RunOpt.KeepTest=0;
end
RunOpt.EGID=[];i1=strfind(Cam,'-egid');if ~isempty(i1)
 [RunOpt.EGID,i2,i3,i4]=sscanf(CAM(i1+5:end),'%i',1);
 CAM(i1:i1+3+i4)='';[CAM,Cam]=comstr(CAM,1);
end

NNode=sparse(model.Node(:,1),1,1:size(model.Node,1));
el0.OrigNodeId=el0.Node(:,1);
if isfield(el0,'bas'); [el0.Node,el0.bas]=basis(el0.Node,el0.bas);end
FEnode=model.Node; FEelt=model.Elt;
i1=strfind(Cam,'epsl'); if ~isempty(i1);
 [epsl,i2,i3,i4]=sscanf(Cam(i1+4:end),'%g');
 CAM(i1+[0:i4+2])=''; [CAM,Cam] = comstr(CAM,1);
end

if isempty(el0); opt(2)=Inf;  %  - - - - - - - - - - - - - - - -  same nodes 
elseif comstr(Cam,'combine') % combine node sets. Nodes with same number must
% be at same location or error

i2=find(~ismember(el0.Node,model.Node,'rows'));
[r3,i4,i5]=intersect(el0.Node(i2,1),model.Node(:,1));
if ~isempty(r3); 
 for j1=1:length(i4)
  r3(j1)=norm(el0.Node(i2(i4(j1)),:)-model.Node(i5(j1),:));
  if r3(j1)<epsl; i4(j1)=0;i5(j1)=0;r3(j1)=0; end
 end
 i4=i4(find(i4));i5=i5(find(i5));
 if ~isempty(i4)
  sdtw('_clip 60 1','%i ',[model.Node(i5,1)]')
  error('Some nodes with identical numbers are different');
 end
else % append new nodes
 FEnode(end+[1:length(i2)],:)=el0.Node(i2,:);
end
nind=sparse(FEnode(:,1),1,1:size(FEnode,1)); i1=nind(el0.Node(:,1));
opt=[-1 0];

elseif comstr(Cam,'merge')  %  'AddTestMerge' merge nodes - - - - - - - - - -

 [FEnode,i1]=feutil(sprintf('AddNode epsl %.15g',epsl),model.Node,el0.Node);
 r2=FEnode(i1,5:7)-el0.Node(:,5:7);if max(sqrt(sum(r2.^2,2)))>epsl*sqrt(3)
  FEnode=model.Node;
  [FEnode,i1]=feutil(sprintf('AddNode epsl %.15g',epsl),FEnode,FEnode);
  if length(unique(i1))<size(FEnode,1);
   sdtw('_nb','Renumbering non-merged nodes in initial model');
   model=feutil('renumber',model,FEnode(i1,1)); FEnode=model.Node;
   [FEnode,i1]=feutil(sprintf('AddNode epsl %.15g',epsl),FEnode,el0.Node);
  else;
   error('This case is seen as a bug');
  end
 end
 % attempt to use el0 node numbers if NodId is free
 [i2,i3]=setdiff(FEnode(i1,1),model.Node(:,1));
 [i2,i4]=setdiff(el0.Node(i3,1),model.Node(:,1));
 FEnode(i1(i3(i4)),1)=el0.Node(i3(i4),1);
 opt=[-1 0]; 

else                        %  - - - - - - - - - - - - - - - -  node shift 
 r2 = max(el0.Node(:,1)); if length(NNode)<r2 NNode(r2)=0;end
 opt = comstr(CAM,[-1 0]); opt(2)=0; i1=NNode(el0.Node(:,1)); 
 if opt(1)==0
   if all(i1) & norm(FEnode(i1,5:7)-el0.Node(:,5:7))<epsl 
     sdtw('_nb','AddTest : test nodes kept unchanged');opt(2)=Inf;
   elseif any(i1) opt(1)= max(FEnode(:,1)); end
 end
end


if isinf(opt(2))  % same nodes
  i1=NNode(el0.Node(:,1));
elseif opt(1)==-1
elseif opt(1)~=0  % node shift
  el0.Node(:,1)=el0.Node(:,1)+opt(1);
  if length(NNode)<max(el0.Node(:,1)) NNode(max(el0.Node(:,1)))=0;end
  i1=NNode(el0.Node(:,1));
  if all(i1) & norm(FEnode(i1,5:7)-el0.Node(:,5:7))<epsl 
    warning('test nodes kept unchanged');opt(2)=Inf;
  end
end


if ~isempty(FEnode)&~finite(FEnode(end,5)) FEnode=FEnode(1:end-1,:);end
% NNode giving old node positions in i1 (new node numbers)
% i1 gives the associated node numbers
if opt(1)==-1 % merged nodes
  NNode=[]; NNode=sparse(el0.Node(:,1),1,1:size(i1,1));
elseif finite(opt(2))&finite(opt(1))  % new test nodes added
  i1 = size(FEnode,1)+[1:size(el0.Node,1)]';FEnode=[FEnode;el0.Node];
  NNode=[];NNode(el0.Node(:,1)-opt(1))=1:size(el0.Node,1);
elseif ~isempty(el0.Node)&~finite(opt(2))
  NNode=[];NNode(el0.Node(:,1)-opt(1))=1:size(el0.Node,1);
end

elt=el0.Elt;
i4=size(FEelt,1);
if i4&isempty(strfind(Cam,'keeptest'));
 FEelt=feutil('removeelt EGID -1',FEnode,FEelt);
 if size(FEelt,1)~=i4
    sdtw('_nb','Removed test wire-frame element groups');
 end
end


if isempty(elt)
 FEelt(end+1,1:8)=[Inf abs('mass2') 0 -1];
 FEelt(end+[1:size(el0.Node,1)],1)=el0.Node(:,1);
elseif finite(elt(1,1)) % this is a trace line
 i3=find(elt(:,83:end));
 i2=elt(:,83:end);i3=find(i2); i2(i3)=FEnode(i1(NNode(i2(i3))),1);
 elt(:,83:end)=i2; 
 FEel0=feutil('trace2elt',elt); femesh('addsel',';');

else
    [EGroup,nGroup]=getegroup(elt);
    for jGroup=1:nGroup
     ElemF= feutil('getelemf',elt(EGroup(jGroup),:),jGroup);
     cEGI = EGroup(jGroup)+1:EGroup(jGroup+1)-1;
     if ~isempty(RunOpt.EGID); 
      elt(EGroup(jGroup),length(ElemF)+3)=RunOpt.EGID;
     end
     i2=fe_super('node',ElemF); % node indices
     elt(cEGI,i2)=reshape(FEnode(i1(NNode(elt(cEGI,i2))),1), ...
                          length(cEGI),length(i2));
    end
    FEelt(end+[1:size(elt,1)],1:size(elt,2))=elt;
end

out1=i1;  % Return new nodes
% Attempt to merge the model information
try;
 model.pl(end+[1:size(el0.pl,1)],1:size(el0.pl,2))=el0.pl;
 model.pl=unique(model.pl,'rows');
end
try;
 model.il(end+[1:size(el0.il,1)],1:size(el0.il,2))=el0.il;
 model.il=unique(model.il,'rows');
end
try; model.Stack(end+[1:size(el0.Stack,1)],:)=el0.Stack; end

% Refine the OrigNumbering Output
if ~isequal(out1,el0.OrigNodeId)
 i1=stack_get(model,'info','OrigNumbering','getdata');
 i1=[i1;int32([el0.OrigNodeId FEnode(out1,1)])];
 model=stack_set(model,'info','OrigNumbering',i1);

 if isfield(el0,'tdof'); % Shift node numbers in tdof definition
   nind=sparse(el0.OrigNodeId,1,FEnode(out1,1));
   model.tdof=el0.tdof; 
   if size(el0.tdof,2)==1 
     model.tdof=full(nind(fix(el0.tdof)))+rem(el0.tdof,1);
   else
     model.tdof(:,2)=full(nind(fix(el0.tdof(:,2))))
   end
 end

end



out=model; out.Elt=FEelt; out.Node=FEnode;

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

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

r1=varargin{carg};carg=carg+1;
r2=varargin{carg};carg=carg+1;

if size(r1,2)==7 out1=r1(:,1); r1=r1(:,5:7); else; out1=[1:size(r1,1)]';end
if size(r2,2)==7 i2=r2(:,1); r2=r2(:,5:7); else; i2=[1:size(r2,1)]';end

i1=zeros(size(r1,1),1);
for j1=1:size(r1,1)
 r3=(r2-r1(ones(size(r2,1),1)*j1,:)).^2*[1;1;1];
 [r3,i1(j1)]=min(r3);
 out1(j1,2)=i2(i1(j1));
end
i1
out=i1;

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

 elt = varargin{carg};carg=carg+1;
 if isfield(elt,'Elt') elt=elt.Elt; 
 elseif finite(elt(1))&size(elt,2)==7 elt = varargin{carg};carg=carg+1;end
 out=[1:size(elt,1)]';

 [EGroup,nGroup]=getegroup(elt);
 for jGroup=1:nGroup
   [ElemF,opt]= feutil('getelemf',elt(EGroup(jGroup),:),jGroup);
   [i2,i3]=fe_super('prop',ElemF);
   if i3==2 out(EGroup(jGroup))=0; end
 end

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

elt = varargin{carg};carg=carg+1;
if isfield(elt,'Elt') elt=elt.Elt; end
if finite(elt(1))&size(elt,2)==7; elt = varargin{carg};carg=carg+1;end

i1=strfind(Cam,'-elt'); if ~isempty(i1); 
 RunOpt.Back='elt';CAM(i1+[0:3])='';[CAM,Cam]=comstr(CAM,1);
else; RunOpt.Back='eltid';
end
if comstr(Cam,'from')
 i1=comstr(Cam(5:end),[-1]); if isempty(i1); error('Not a valid EltIdFrom');end
 RunOpt.Mode='Renumber'; RunOpt.From=i1;  RunOpt.EltId=[];

elseif carg<=nargin; % return a string list by EltId
  RunOpt.EltId= varargin{carg};carg=carg+1;
  i3=find(RunOpt.EltId);cind(RunOpt.EltId(i3))=i3; st='';
  RunOpt.Mode='list';
else; RunOpt.EltId=[];st='';RunOpt.Mode='nominal';
end

i1=strfind(Cam,'silent'); if ~isempty(i1)&~isempty(Cam)
 RunOpt.Silent=1;Cam(i1+[0:5])='';CAM(i1+[0:5])='';
elseif any(Cam==';');
 RunOpt.Silent=1;Cam(find(Cam==';'))='';CAM(find(CAM==';'))='';
else; RunOpt.Silent=0; end

[EGroup,nGroup]=getegroup(elt); eltrows=finite(elt(:,1));

out=zeros(size(elt,1),1);out1=out;
for jGroup=1:nGroup 
   [ElemF,opt]= feutil('getelemf',elt(EGroup(jGroup),:),jGroup); 
   cEGI=EGroup(jGroup)+1:EGroup(jGroup+1)-1; 
   [i2,i3]=fe_super('prop',ElemF);
   if i3==1 % Unique superelement
    i2=min(find(elt(EGroup(jGroup),:)==0)); 
    if ~isempty(i2)&size(elt,2)>=i2+2 
      out(EGroup(jGroup))=elt(EGroup(jGroup),i2+2);
    elseif isempty(i2); i2=size(elt,2);
    end
    out1(EGroup(jGroup))=i2+2;eltrows(EGroup(jGroup))=1;
   else % Generic element/superelement
    if length(i2)<3|~i2(3)   out1(cEGI)=0;out(cEGI)=0; 
    elseif i2(3)>size(elt,2) out1(cEGI)=i2(3);out(cEGI)=0; 
    else out1(cEGI)=i2(3);out(cEGI)=elt(cEGI,i2(3));  end 
    if ~isempty(RunOpt.EltId) & strcmp(RunOpt.Mode,'list')
      i2=out(cEGI); i3=find(i2); 
      if max(i2(i3))>length(cind); cind(max(i2(i3)))=0;end 
      if any(cind(i2(i3))) 
        st=sprintf('%s\n%s (%i)',st,ElemF,opt(1)); 
        for j1=find(cind(i2(i3)))' 
          st=sprintf('%s\n%s',st,sprintf('%i ',elt(cEGI(i3(j1)),:))); 
        end 
      end 
    end 
   end % Superelement type 
end 

i4=find(sparse(out+1,1,(out~=0)*1.1)>1.1)-1;
if isempty(strfind(Cam,'fix'))&strcmp(RunOpt.Mode,'nominal');

 if ~isempty(i4)
  ind=1:length(i4); 
  if length(ind)>100; ind=ind([1:10 end+[-9:0]]);
   fprintf('Displaying first 10 and last 10 only'); 
  end
  for j1=ind
   fprintf('EltID %i ( %s)\n',i4(j1),sprintf('%i ',find(out==i4(j1))))
  end
  warning(['Repeated EltID ' sprintf(' %i',i4)]);
 end

else % fix

% i2=find(RunOpt.EltId); RunOpt.EltId(i2)=i1+[0:length(i2)-1]';
  if strcmp(RunOpt.Mode,'Renumber');
   i4=find(out); out(i4)=RunOpt.From+[0:length(i4)-1]';
  elseif ~isempty(i4) 
   if ~RunOpt.Silent; fprintf('Fixed repeated eltid\n');end
   [i5,i4]=sort(out); i1=find(diff(i5)==0&i5(1:length(i5)-1)~=0)+1;
   out(i4(i1))=max(out)+[1:length(i1)]';
  end
  i1=find(out==0&eltrows); % fix zero values
  if ~isempty(i1)
    if ~RunOpt.Silent;fprintf('Fixed zero eltid\n');end
    out(i1)=max(out)+[1:length(i1)]';
  end

  i2=find(out1); i3=i2+size(elt,1)*(out1(i2)-1);
  if prod(size(elt))<max(i3) elt(1,ceil(max(i3)/size(elt,1)))=0; end
  elt(i2+size(elt,1)*(out1(i2)-1))=out(i2);
  out1=elt;
end
if nargout==0 & ~isempty(st); disp(st);
elseif strcmp(RunOpt.Back,'elt'); out=elt;
end

% ------------------------------------------------------------------------
elseif comstr(Cam,'egid')

 elt = varargin{carg};carg=carg+1;
 if isfield(elt,'Elt'); elt=elt.Elt; 
 elseif finite(elt(1))&size(elt,2)==7; elt = varargin{carg};carg=carg+1;end

 if carg<=nargin opt=varargin{carg}; carg=carg+1; else opt=[];end
 % Specify name transcription table
 if ~isempty(strfind(Cam,'-vol'));
   RunOpt.NameTable={'quad4','q4p';'tria3','t3p';'quadb','q8p';'tria6','t6p'};
 else;RunOpt.NameTable={};
 end
 if carg<=nargin; out1=varargin{carg}; carg=carg+1; % eltname
 else; out1={};
 end 
 [EGroup,nGroup]=getegroup(elt);
 if ~isempty(opt)&length(opt)~=nGroup
  error('you must provide nGroup EGID values to set');
 end
 if ~isempty(out1)&length(out1)~=nGroup
  error('you must provide nGroup EltName values to set');
 end

 out=zeros(nGroup,1);
 for jGroup=1:nGroup
   [ElemF,i1]= feutil('getelemf',elt(EGroup(jGroup),:),jGroup);
   if ~isempty(RunOpt.NameTable)
    i2=strmatch(ElemF,RunOpt.NameTable(:,1),'exact');
    if ~isempty(i2);
     out1{jGroup}=RunOpt.NameTable{i2,2};
     elt(EGroup(jGroup),1:length(out1{jGroup})+2+length(i1))= ...
       [Inf abs(out1{jGroup}) 0 i1(:)'];
    end
   end
   if isempty(i1); i1=jGroup;end; out(jGroup)=i1(1); 
   if length(out1)<jGroup|isempty(out1{jGroup})
    out1{jGroup}=ElemF;
   end
   if ~isempty(opt);
     elt(EGroup(jGroup),1:length(out1{jGroup})+3)= ...
      [Inf abs(out1{jGroup}) 0 opt(jGroup)];
   end
 end
 if ~isempty(opt)|~isempty(RunOpt.NameTable); out=elt;end


% ------------------------------------------------------------------------
elseif comstr(Cam,'matid')|comstr(Cam,'proid')|comstr(Cam,'mpid')

 elt = varargin{carg};carg=carg+1;
 if isempty(elt)
 elseif isfield(elt,'Elt') elt=elt.Elt; 
 elseif finite(elt(1))&size(elt,2)==7 elt = varargin{carg};carg=carg+1;
 end
 [EGroup,nGroup]=getegroup(elt);
 if carg<=nargin r1=varargin{carg};carg=carg+1;out=elt;
 else r1=[];out=zeros(size(elt,1),2); end

 for jGroup=1:nGroup
   cEGI=EGroup(jGroup)+1:EGroup(jGroup+1)-1;
   [ElemF,opt]= feutil('getelemf',elt(EGroup(jGroup),:),jGroup);
   [i2,i3]=fe_super('prop',ElemF); i2(4)=0;
   if i3~=2 % single superelement
    out(EGroup(jGroup),3)=jGroup;
   elseif ~isempty(r1) % set value
    if i2(1)&any(r1(cEGI,1)); out(cEGI,i2(1))=r1(cEGI,1); end
    if i2(2)&any(r1(cEGI,2)); out(cEGI,i2(2))=r1(cEGI,2); end    
   else
    if i2(1)&size(elt,2)>=i2(1) out(cEGI,1)=elt(cEGI,i2(1)); end
    if i2(2)&size(elt,2)>=i2(2) out(cEGI,2)=elt(cEGI,i2(2)); end
    out(cEGI,3)=jGroup;
   end
 end
 if ~isempty(r1)
 elseif comstr(Cam,'matid')     out=out(:,1);
 elseif comstr(Cam,'proid') out=out(:,2); end

% ------------------------------------------------------------------------
elseif comstr(Cam,'dof2mode')  [CAM,Cam]=comstr(CAM,5);

 def= varargin{carg};carg=carg+1;
 mdof= varargin{carg};carg=carg+1;
 % NodeId ux uz uz rx ry rz
 i1=find(sparse(round(mdof),1,1));
 r2=fe_c(mdof,[i1+.01;i1+.02;i1+.03],'place');
 r3=fe_c(mdof,[i1+.04;i1+.05;i1+.06],'place');

 out=zeros(length(i1),4,size(def,2));
 if nnz(r3) r2=[r2;r3];i2=6; else i2=3; end

 for j2=1:size(def,2)
  out(:,1,j2)=i1; 
  out(:,2:i2+1,j2)=reshape([r2]*full(def(:,j2)),length(i1),i2);
 end

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

mode=varargin{carg}; carg=carg+1;
if ~isa(mode,'cell');mode={mode};end 
if isa(mode,'cell')
 def=[]; dof=[]; 
 r1=mode{1};
 i1=r1(:,1)';
 if size(r1,2)==4 dof=[i1+.01;i1+.02;i1+.03];
 else dof=[1:size(r1,2)-1]'/100;
      dof=dof(:,ones(size(i1)))+i1(ones(size(dof)),:);
 end
 dof=dof(:); i1=i1(:); def=zeros(length(dof),length(mode));
 def(:,1)=reshape(r1(:,2:end)',length(dof),1);
 for j1=2:length(mode)
    r1=mode{j1};
    if norm(r1(:,1)-i1) error('Variable node not supported'); end
    def(:,j1)=reshape(r1(:,2:end)',length(dof),1);
 end
end
if ~isempty(strfind(Cam,'struct'));out=struct('def',def,'DOF',dof);
else; out=def; out1=dof; 
end

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

[carg,FEnode,FEelt,FEel0,ModelStack]=get_nodeelt(varargin,carg,ModelStack);

% 'Divide InGroups'
if comstr(Cam,'in')

[EGroup,nGroup]=getegroup(FEelt);
NNode=sparse(FEnode(:,1),1,1:size(FEnode,1));
elt=[];

for jGroup = 1:nGroup %loop on element groups
   ElemF= feutil('getelemf',FEelt(EGroup(jGroup),:),jGroup);
   cEGI = EGroup(jGroup)+1:EGroup(jGroup+1)-1;
   i1=fe_super('node',ElemF); % node indices
   i2 = [length(i1) length(cEGI) length(i1)*length(cEGI)];
   i3 = 1:i2(2);
   % build matrix of node use for elements

   i2=[2:length(i1) 1];
   i2 = sparse(NNode(FEelt(cEGI,i1)),NNode(FEelt(cEGI,i2)), ...
               ones(length(i1)*length(cEGI),1));
   i2=triu(i2+i2'); i3=etree(i2);
   i4=find(~i3); i3(i4)=-[1:length(i4)];
   i4=find(i3>0);while ~isempty(i4) i3(i4)=i3(i3(i4));i4=find(i3>0);end;
   i3=-i3;i4=i3(NNode(FEelt(cEGI,i1(1))));


   for j1=1:max(i4)
     i5=find(i4==j1); if ~isempty(i5)
     elt=[elt;FEelt(EGroup(jGroup),:);FEelt(cEGI(i5),:)];
   end;end
end % of loop on groups
out=elt;
 
  
elseif comstr(Cam,'group')  % 'DivideGroups'

 [opt,CAM,Cam]=comstr(CAM,'group','%i');
 if length(opt)~=1
 error('DivideGroup can only be applied to a single group');
 end

 [EGroup,nGroup]=getegroup(FEelt);
 if opt>nGroup | opt<1 error('Not a valid group number');end

 if carg<=nargin r1 = varargin{carg};carg=carg+1; else r1=[];end
 [i1,elt]=feutil(['findelt' CAM],FEnode, ...
     FEelt(EGroup(opt):EGroup(opt+1)-1,:),FEel0,varargin{2:end});
 if ~isempty(i1) & ~(length(i1)==EGroup(opt+1)-EGroup(opt))
  i2 = [EGroup(opt)+1:EGroup(opt+1)-1]';i2(i1-1)=0;i2=i2(find(i2));
  i1 = [[1:EGroup(opt)]';i1+EGroup(opt)-1;EGroup(opt);i2; ...
        [EGroup(opt+1):size(FEelt,1)]'];
  FEelt = FEelt(i1,:); 
 else warning('selection lead to no division'); end
 out=FEelt;

elseif comstr(Cam,'hexatransition')  % 'DivideHexaTransition'

[EGroup,nGroup]=getegroup(FEelt);
NNode=sparse(FEnode(:,1),1,1:size(FEnode,1));
elt=[];

i1=varargin{carg};carg=carg+1; 
cEGI=feutil('findelt withnode',FEnode,FEelt,[],i1);
% Reference division
xi = [-1 -1 -1;1 -1 -1;1 1 -1;-1 1 -1;-1 -1 1;1 -1 1;1 1 1;-1 1 1;
        -1/3 -1/3 0; 1/3 -1/3 0; 1/3  1/3 0; -1/3  1/3 0; -1/3 -1/3 1;
        1/3 -1/3 1;   1/3  1/3 1; -1/3  1/3 1  ];
i2=[1 5 13 9 2 6 14 10; 2 6 14 10 3 7 15 11; 3 7 15 11 4 8 16 12;
        4 8 16 12 1 5 13 9; 9:16; 1:4 9:12 ];

r1=ones(size(xi,1)*length(cEGI),3); % new nodes
i3=ones(size(i2,1)*length(cEGI),10); % new elts

for jElt=1:length(cEGI)
    nodeE=FEnode(NNode(FEelt(cEGI(jElt),1:8)),5:7);
    r2=integrules('hexa8',[xi zeros(size(xi,1),1)],nodeE);
    r1(16*jElt+[-15:0],1:3)=r2.N*nodeE;
    i3(6*jElt+[-5:0],1:10)=[16*(jElt-1)+i2 ones(6,1)*FEelt(cEGI(jElt),9:10)];
end % jElt
[FEnode,ind]=feutil('addnode',FEnode,r1);
elt(end+[1:size(i3,1)],1:size(i3,2))=[ind(i3(:,1:8)) i3(:,9:10)];


FEelt(cEGI(2:end,:),:)=[];
FEelt=FEelt([1:cEGI(1)-1 cEGI(1)*ones(1,size(elt,1)) cEGI(1)+1:size(FEelt,1)],:);
FEelt(cEGI(1)+[0:size(elt,1)-1],1:10)=elt;


out=FEnode;out1=FEelt;
 

else error('Unknown feutil divide');end

% ------------------------------------------------------------------------
elseif comstr(Cam,'find')  [CAM,Cam]=comstr(CAM,5);

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
% feutil('findnode',FEnode,FEelt,FEel0,varargin)
if comstr(Cam,'node') [CAM,Cam]=comstr(CAM,5);

%i4 current node set, i5 operator positions
ModelStack={}; SetStack={};
[carg,FEnode,FEelt,FEel0,ModelStack]=get_nodeelt(varargin,carg,ModelStack);
NNode=sparse(FEnode(:,1),1,1:size(FEnode,1));
[epsl,CAM,Cam]=test_epsl(epsl,CAM,Cam);

if comstr(Cam,'stack') [CAM,Cam]=comstr(CAM,6);st2='return';else st2=''; end
if carg<=nargin&isa(varargin{carg},'cell') Stack=varargin{carg};carg=carg+1;
elseif isempty(CAM)
 if carg<=nargin; ind=varargin{carg};carg=carg+1;
  out=full(NNode(ind)); out1=FEnode(out,:); return; 
 else; error('Not a valid findnode call');
 end
 
else  Stack=BuildSelStack(CAM,FEnode,FEelt,FEel0,1,varargin{carg:end});
end % of Stack is given as argument

if comstr(st2,'return') out=Stack;return; end

if ~strcmp(Stack{end,1},'}') Stack{end+1,1}='}';end
out=[]; j1=0;

% loop on stack - - - - - - - - - - - - - - - - - - - - - - - - - - -
while j1 <size(Stack,1)-1

 j1=j1+1; Cam=comstr(Stack{j1,2},-27); Bole=Stack{j1,1};

 % finds nodes betwen the planes defined by two nodes
 if comstr(Cam,'between') [st,Cam]=comstr(Cam,8);

  if length(Stack{j1,4})~=2 error('improper between entry'); end
  r1=FEnode(NNode(Stack{j1,4}),5:7);r1(2,:)=r1(2,:)-r1(1,:);
  r2=norm(r1(2,:)); if r2 r1(2,:)=r1(2,:)/r2;end
  r1=(FEnode(:,5:7)-r1(ones(size(FEnode,1),1),:))*r1(2,:)';
  i4=find(r1+epsl>=0&r1-epsl<=r2);

 % finds nodes in a group  ggggggggggggggggggggggggggggggggggggggggggggg
 elseif comstr(Cam,'group') [st,Cam]=comstr(Cam,6);

 % model description matrix can be given as second argument to femesh
 if ~isempty(Cam)&Cam(1)=='a' elt=FEel0; else elt=FEelt; end
 [EGroup,nGroup]=getegroup(elt);
 if strcmpi(Stack{j1,4},'all') opt=1:nGroup;
 else opt=IEqTest(1:nGroup,Stack{j1,3},Stack{j1,4}); end

 i1=[];for jGroup=opt; if jGroup>0&jGroup<=nGroup
   ElemF= feutil('getelemf',elt(EGroup(jGroup),:),jGroup);
   [i2,SEopt]=fe_super('node',ElemF);
   if SEopt(1,1)==2 i2 = elt(EGroup(jGroup)+1:EGroup(jGroup+1)-1,i2);end
   i1 =[i1;i2(:)];
 end;end
 i4 = find(sparse(i1+1,1,i1))-1;

 % GID Node group identification ------------------------------------
 elseif comstr(Cam,'gid') [st,Cam]=comstr(Cam,4);

   opt=Stack{j1,4}; if ischar(opt) opt=str2num(opt);end
   st=Stack{j1,3}; i1=FEnode(:,4);
   if comstr(st,'~='); i4=find(~ismember(i1,opt));
   elseif comstr(st,'>=') i4=find(i1>=opt(1));
   elseif comstr(st,'<=') i4=find(i1<=opt(1));
   elseif comstr(st,'>')  i4=find(i1>opt(1));
   elseif strcmp(st,'<')  i4=find(i1<opt(1));
   else 
    i4=find(ismember(i1,opt));
   end
   i4=FEnode(i4,1);

 % NodeId eeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeee
 elseif comstr(Cam,'nodeid') [st,Cam]=comstr(Cam,7);

   opt=Stack{j1,4}; if ischar(opt) opt=str2num(opt);end
   st=Stack{j1,3}; i1=FEnode(:,1);
   if comstr(st,'~=') 
     if max(opt)>length(NNode) NNode(max(opt))=0; end
     i4=full(NNode(opt));i4=i4(find(i4)); i1=1:size(FEnode,1);
     if ~isempty(i4)  i1(i4)=0; end
     i4=find(i1);
   elseif comstr(st,'>=') i4=find(i1>=opt(1));
   elseif comstr(st,'<=') i4=find(i1<=opt(1));
   elseif comstr(st,'>')  i4=find(i1>opt(1));
   elseif strcmp(st,'<')  i4=find(i1<opt(1));
   else
     if max(opt)>length(NNode) NNode(max(opt))=0; end
     i4=full(NNode(opt));i4=i4(find(i4));
   end
   i4=FEnode(i4,1);

 % Egid eeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeee
 elseif comstr(Cam,'egid')|comstr(Cam,'eltid')|comstr(Cam,'eltname') ...
        | comstr(Cam,'matid') | comstr(Cam,'proid')

   [r1,elt]=feutil('findelt',FEnode,FEelt,FEel0,Stack(j1,:));
   if isempty(elt) 
     warning(sprintf('No element selected for %s',Cam)); i4=[];
   else 
     i4=feutil('findnode',FEnode,elt,[],{'','Group','==','All'});
   end
 
 % inelt eeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeee
 elseif comstr(Cam,'inel') 

  ind=SubStackIndex(Stack,j1+1);

  if  Cam(5)=='t'
   model=struct('Node',FEnode,'Elt',FEelt,'El0',FEel0,'Stack',[]);
   model.Stack=ModelStack;[r1,elt]=feutil('findelt',model,Stack(ind,:));
  elseif  Cam(5)=='0'
   model=struct('Node',FEnode,'Elt',FEel0,'El0',[],'Stack',[]);
   model.Stack=ModelStack;[r1,elt]=feutil('findelt',model,Stack(ind,:));
  else error('Not a valid InElt command')
  end
  if isempty(elt) error('No element selected for InElt'); end
  i4=feutil('findnode',FEnode,elt,[],{'','Group','==','All'});
  j1=max(ind);

 % notin eeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeee
 elseif comstr(Cam,'notin') 

  ind=SubStackIndex(Stack,j1+1);

  [r1,elt]=feutil('findelt',FEnode,FEelt,FEel0,Stack(ind,:));
  if isempty(elt) error('No element selected for NotIn'); end
  ind=FEnode(:,1);
  i4=feutil('findnode',FEnode,elt,[],{'','Group','==','All'});
  ind(full(NNode(i4)))=0; i4=ind(find(ind));
  j1=max(ind);

 % finds nodes on a plane  uuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuu
 elseif comstr(Cam,'plane') [st,Cam]=comstr(Cam,6);

   st=Stack{j1,3};opt=Stack{j1,4}; if ischar(opt) opt=str2num(opt);end
   i3=(FEnode(:,5:7)-opt(ones(size(FEnode,1),1),1:3))*opt(4:6)';

   if     comstr(st,'>=') i4=find(i3>-epsl);
   elseif comstr(st,'<=') i4=find(i3<epsl);
   elseif comstr(st,'>')  i4=find(i3>0);
   elseif strcmp(st,'<')  i4=find(i3<0);
   else                   i4=find(abs(i3)<epsl);
   end
   if ~isempty(i4) i4=FEnode(i4,1); end

 % finds nodes in a set defined in the model stack
 elseif comstr(Cam,'set'); [st,Cam]=comstr(Cam,4); st='';

  i4=FeutilMatchSet(ModelStack,Stack(j1,:));

 % finds nodes within a radius  uuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuuu
 elseif comstr(Cam,'rad') [st,Cam]=comstr(Cam,4); st='';

   st=Stack{j1,3};opt=Stack{j1,4}; if ischar(opt) opt=str2num(opt);end
   i2=length(opt);

   if i2<2 error('feutil findnode rad: must specify radius and node');
   elseif i2==2
    i2=find(FEnode(:,1)==opt(2));
    if isempty(i2) opt(2:4)=[0 0 0];else opt(2:4)=FEnode(i2(1),5:7);end
   elseif i2<4 opt(4)=0; end
   i3=FEnode(:,5:7)-opt(ones(size(FEnode,1),1),2:4);
   i3=sum((i3').^2)-opt(1)^2;

   if     comstr(st,'==') i4=find(abs(i3)<epsl);
   elseif comstr(st,'>=') i4=find(i3>=-epsl);
   elseif comstr(st,'>')  i4=find(i3>0);
   elseif strcmp(st,'<')  i4=find(i3<0);
   else                   i4=find(i3<=epsl);
   end
   if ~isempty(i4) i4=FEnode(i4,1); end

 % finds nodes by coordinates u u u u u u u u u u u u u u u
 elseif ~isempty(Cam) & any(Cam(1)=='xyzr')
   st=Stack{j1,3};opt=Stack{j1,4};
   if ischar(opt)&isempty(intersect(opt,'xyzr')); opt=str2num(opt);end
   i2=find('xyzr'==Cam(1));
   if i2<4; r1=FEnode(:,4+i2);else; r1=sqrt(sum(FEnode(:,5:6).^2,2));end

   if ischar(opt);x=FEnode(:,5);y=FEnode(:,6);z=FEnode(:,7);r=r1;
     i4=find(eval(horzcat(Cam,st,opt)));
   elseif comstr(st,'==') i4=find(abs(r1-opt)<epsl);
   elseif comstr(st,'~=') i4=find(abs(r1-opt)>epsl);
   elseif comstr(st,'>=') i4=find(r1-opt>=-epsl);
   elseif comstr(st,'>')  i4=find(r1-opt>0);
   elseif strcmp(st,'<')  i4=find(r1-opt<0);
   else                   i4=find(r1-opt<=epsl);
   end
   if ~isempty(i4) i4=FEnode(i4,1); end

 elseif strcmp(Stack{j1,1},'}')
 else % give coordinates  with possible infinity

   opt=Stack{j1,4}; if ischar(opt) opt=str2num(opt);end
   i4=1:size(FEnode,1);
   if finite(opt(1,1)) i4=i4(find(abs(FEnode(i4,5)-opt(1))<epsl)); end
   if ~isempty(i4)& finite(opt(1,2))
         i4 = i4(find(abs(FEnode(i4,6)-opt(2))<epsl)); end
   if ~isempty(i4)& finite(opt(1,3))
         i4 = i4(find(abs(FEnode(i4,7)-opt(3))<epsl)); end
   i4=FEnode(i4,1);

 end %  u u u u u u u u u u u u u u u

 % bolean operations on node sets

 if isempty(out) out=i4;
 elseif strcmp(Bole,'&') % keep those present in all
   if isempty(i4) out=i4;
   else
    i4=find(sparse(i4,ones(size(i4)),ones(size(i4))));
    i4=sparse(i4,1,i4);if max(out)>length(i4) i4(max(out))=0;end
    out = nonzeros(i4(out));
   end
 elseif strcmp(Bole,'|') % keep those present in any
   out=[i4;out];   out=find(sparse(out,ones(size(out)),ones(size(out))));
 end
end % of j1 loop on arguments

if ~isempty(out)&out(end)==size(FEnode,1)&any(~finite(FEnode(end,5:7)))
  out=out(1:end-1);
end
if nargout==2 % return the nodes too
  NNode=sparse(FEnode(:,1),1,1:size(FEnode,1));
  if ~isempty(out); out1=FEnode(NNode(out),:);else out1=[];end
end

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
% [ind,elt,IndWithHeaders]=feutil('findelt',FEnode,FEelt,FEel0)
elseif comstr(Cam,'el') [CAM,Cam]=comstr(CAM,4);

%i4 current element set, out final element set, i5 operator positions
%i7 used groups

ModelStack={};
[carg,node,elt,el0,ModelStack]=get_nodeelt(varargin,carg,ModelStack);

if comstr(Cam,'stack');[CAM,Cam]=comstr(CAM,6);st2='return';else;st2=''; end
if isempty(elt)&isempty(el0)
 Stack={'','Empty model'}; out2=[]; nGroup=0;
elseif carg<=nargin&isa(varargin{carg},'cell')
  Stack=varargin{carg};carg=carg+1;
  [EGroup,nGroup]=getegroup(elt); out2='';
else % stack is not given
 if isempty(CAM)&carg<=nargin&ischar(varargin{carg})
   [CAM,Cam]=comstr(varargin{carg},1); carg=carg+1;
 end
 out2=''; Stack=BuildSelStack(CAM,node,elt,el0,1,varargin{carg:end});
end % is stack given
if strcmp(st2,'return'); out=Stack;return; end
if ~strcmp(Stack{end,1},'}') & ~strcmp(Stack{1,2},'Empty model')
  Stack{end+1,1}='}';
end
j1=0; out=[];
while j1<size(Stack,1)-1 % loop on stack - -  - - - - - - - - - - - - - -

 j1=j1+1; Bole=Stack{j1,1}; 
 st=Stack{j1,2};opts=Stack{j1,3};opt=Stack{j1,4};

 % Node Based selection
 if ischar(st)&~isempty(strfind(st,'Node'))
  if strcmp(Stack{j1,3},'{}')
    ind=SubStackIndex(Stack,j1+1);
    model=struct('Node',node,'Elt',elt,'El0',el0,'Stack',[]);
    model.Stack=ModelStack;
    if strcmp(opts,'{}') opt=feutil('findnode',model,Stack(ind,:)); end
    j1=max(ind);
  else; opt=Stack{j1,4};
  end
  if any(opt<0); NNode=(sparse(-opt+1,1,1:length(opt)));
  else NNode=(sparse(opt+1,1,1:length(opt)));end
 end
 % reindexing for options

 i4=[];
 if comstr(comstr(st,-27),'seledge')
    if j1>1  i7=[];for jGroup = 1:nGroup
      if any(out>EGroup(jGroup)&out<EGroup(jGroup+1)) i7(end+1)=jGroup;end
     end
     out2=sort([out;EGroup(i7)']); elt=elt(out2,:);
    end
    [elt,CAM]=feutil(['getedgeline' Stack{j1,4}],node,elt);
    out=find(finite(elt(:,1))); [EGroup,nGroup]=getegroup(elt);out2='edge';

 elseif comstr(comstr(st,-27),'set')

   elt=FeutilMatchSet(ModelStack,Stack(j1,:),elt);
   [EGroup,nGroup]=getegroup(elt); 
    if isempty(elt) error('No face/edge was selected');end
    out=find(finite(elt(:,1))); out2='set';

 elseif comstr(comstr(st,-27),'selface')

    if j1>1  
     i7=[];
     for jGroup = 1:nGroup
      if any(out>EGroup(jGroup)&out<EGroup(jGroup+1)) i7(end+1)=jGroup;end
     end
     out2=sort([out;EGroup(i7)']); elt=elt(out2,:);
    end
    if isempty(elt) error('No element to select'); end
    [elt,CAM]=feutil('getedgepatch',node,elt,[],Stack{j1,4});

    [EGroup,nGroup]=getegroup(elt); 
    if isempty(elt) error('No face was selected');end
    out=find(finite(elt(:,1))); out2='edge';

 else % standard element selection (loop on groups)

 [EGroup,nGroup]=getegroup(elt);

 for jGroup = 1:nGroup %loop on element groups
   [ElemF,EGID]= feutil('getelemf',elt(EGroup(jGroup),:),jGroup);EGID=EGID(1);
   [i2,i6]=fe_super('prop',ElemF);
   cEGI=EGroup(jGroup)+1:EGroup(jGroup+1)-1; i3 = [];

   switch comstr(st,-27)
   case 'group' % group
      i3=IEqTest(jGroup,opts,opt);
      if ~isempty(i3)
          i3 = [1:EGroup(jGroup+1)-EGroup(jGroup)-1];
          if i6(1)==1 i4(end+1,1)=EGroup(jGroup); end
      end
   case 'egid' % egid
      i3=IEqTest(EGID(1),opts,opt);
      if ~isempty(i3)
          i3 = [1:EGroup(jGroup+1)-EGroup(jGroup)-1];
          if i6(1)==1 i4(end+1,1)=EGroup(jGroup); end
      end
   case 'eltind' % eltind 
       if i6~=2 i3=IEqTest(EGroup(jGroup),opts,opt); 
       else i3=IEqTest(cEGI,opts,opt);end 
   case 'eltid' % eltind
       eltid=feutil('eltid;',elt);
       if i6~=2 i3=IEqTest(eltid(EGroup(jGroup)),opts,opt);
       else i3=IEqTest(eltid(cEGI),opts,opt); end
   case 'eltname' % eltname
      if strcmp(opts,'~=')
        if ~comstr(ElemF,opt)
          i3 = [1:EGroup(jGroup+1)-EGroup(jGroup)-1];
          if i6(1)==1 i4(end+1,1)=EGroup(jGroup); end
        end
      elseif isempty(opts)
        if comstr(ElemF,opt)
          i3 = [1:EGroup(jGroup+1)-EGroup(jGroup)-1];
          if i6(1)==1 i4(end+1,1)=EGroup(jGroup); end
        end
      end
   case 'facing' % facing Node angle

    i3=fe_super('face',ElemF);
    if length(opt)<4 opt(4)=0; end
    if size(i3,1)==1 
     [nor,cg]=feutil('getnormal',node, ...
      elt(EGroup(jGroup):EGroup(jGroup+1)-1,:));
     cg=opt(ones(size(cg,1),1),2:4)-cg;r3=sqrt(sum(cg.*cg,2));
     i3=find(any(nor,2));r3=sum(nor(i3,:).*cg(i3,:),2)./r3(i3);
     switch opts
     case '>';  i3=i3(find(r3>opt(1)));
     case '>=';  i3=i3(find(r3>=opt(1)));
     case '<';  i3=i3(find(r3<opt(1)));
     case '<=';  i3=i3(find(r3<=opt(1)));
     case '<>';  i3=i3(find(abs(r3)>opt(1)));
     case '><';  i3=i3(find(abs(r3)<opt(1)));
     otherwise; error('Not a supported opts value');
     end
     i3=i3-1;
    else
     warning('Facing selection can only be applied to flat elements')
     i3=[];
    end

   case 'matid' % mat
     if i2(1)~=0
      if i2(1)>size(elt,2); i3=[]; % No matid is given
      else
        i2=elt(cEGI,i2(1)); i2(find(rem(i2,1)))=0; % integer only
        i3=IEqTest(i2,opts,opt);
      end
     end
   case 'proid' % pro
     if i2(2)==0|i2(2)>size(elt,2);
     else;  
       i2=elt(cEGI,i2(2)); i2(find(rem(i2,1)))=0; % integer only
       i3=IEqTest(i2,opts,opt);
     end
   % with/without/in node single superelement
   case {'withnode','withoutnode','innode'} 

    if i6(1)==1
     i2=fe_super('node',ElemF);
     if max(max(i2))+1>length(NNode) NNode(max(max(i2))+1)=0;end
     i2 = reshape(full(NNode(i2+1)),size(i2,1),size(i2,2))';
     if     comstr(st,'WithNode') & any(i2)  i4(end+1,1)=EGroup(jGroup); 
     elseif comstr(st,'InNode')   &  all(i2)  i4(end+1,1)=EGroup(jGroup); 
     elseif comstr(st,'Without')  & ~any(i2) i4(end+1,1)=EGroup(jGroup);
     end
   % with/in node standard element
    else
     i2 = elt(EGroup(jGroup)+1:EGroup(jGroup+1)-1,fe_super('node',ElemF));
     if max(max(i2))+1>length(NNode); NNode(max(max(i2))+1)=0;end
     i2 = reshape(full(NNode(i2+1)),size(i2,1),size(i2,2));
     if     comstr(st,'WithNode') i3 = find(any(i2,2));  % with
     elseif comstr(st,'InNode')   i3 = find(all(i2,2));  % in
     elseif comstr(st,'Without')  i3 = find(~any(i2,2)); % without
     end
    end
   otherwise error(sprintf('''%s'' not a supported FindElt command',st)); end
   if ~isempty(i3)  i4=[i4;cEGI(i3(:))']; end
 end% of jGroup loop

 % bolean operations on element sets
 if isempty(out) out=i4;
 elseif Bole=='&'   % keep those present in all
   if isempty(i4) out=i4;
   else
    i4=find(sparse(i4,ones(size(i4)),ones(size(i4))));
    i4=sparse(i4,1,i4);if max(out)>length(i4) i4(max(out))=0;end
    out = nonzeros(i4(out));
   end
 elseif Bole=='|'   % keep those present in any
   out=[i4;out];   out=find(sparse(out,ones(size(out)),ones(size(out))));
 end

 end % or seledge

end % of j1 loop on Stack

if nargout>1 % return the elements too
  i7=[];for jGroup = 1:nGroup
   if any(out>EGroup(jGroup)&out<EGroup(jGroup+1)) i7(end+1)=jGroup;end
  end
  if comstr(out2,'edge') 
   out2=sort([out;EGroup(i7)']); out1=elt(out2,:);
   out=[]; out2=[1:size(elt,1)]';
  elseif nGroup==0 out2=[]; out1=[];
  else out2=sort([out;EGroup(i7)']); out1=elt(out2,:); end
elseif comstr(out2,'edge')
  warning('With Edge the selection is not a sub-part of the model');
end

else out='unknown'; end % subcommand selection - - - - - - - - - - - - - -

% ------------------------------------------------------------------------
% recombines a constraint equation so that MASTER DOFs can be defined
elseif comstr(Cam,'fixmpcmaster')

 c=varargin{carg};carg=carg+1;
 if carg<=nargin;  % an initial list of masters is provided
  out1=varargin{carg};carg=carg+1;out1=out1(:);
  [i1,i2]=find(c(:,out1)-speye(size(c,1)));
  if isempty(i1); out=c; return;end % the slave set is consistent
  i1=unique(i1);
  % masters that are not correct and dependent of incorrect ones
  ind=[];while ~isequal(i1,ind);
   ind=i1; i1=find(any(c(:,find(any(c(ind,:),1))),2));
  end
  % consider what is left separately
  [i2,i3]=setdiff(1:size(c,1),i1);
  out=struct('i0',i3,'slave0',out1(i3),'c0',c(i3,:));c=c(i1,:);
  try;
   out1(i3)=[];
   [l,u]=lu(c(:,out1),0);
   c1=u\ ( l\c);
   if norm(c1(:,out1)-speye(length(out1)),'inf')>1e-6
    fprintf('Independence problem on MPC');error('Go to catch');
   end
   c1(:,out1)=speye(length(out1));
   out1=[out.slave0;out1(:)]; out=[out.c0;c1];return;
  catch; out1=[];
  end
 else; out=[]; out1=[];
 end


 ind=find(any(c,1)); r1=c(:,ind)';
 [r2,i2]=max(abs(r1),[],1); % check that there is no obvious solution
 if length(unique(i2))==size(r1,2) % one max term per colum, all independent
  [l,u]=lu(r1(i2,:)); r1=(r1/u)/l;
  out1=ind(i2);
 else

  for j1=1:size(c,1); % select master and eliminate from constraints 
         [r2,i2]=max(abs(r1(:,j1))); out1(end+1)=ind(i2);
         r1(:,j1)=r1(:,j1)/r1(i2,j1);
         r2=r1(i2,:);r2(j1)=0;r1=r1-r1(:,j1)*r2;
         r2=zeros(size(r2));r2(j1)=1;r1(i2,:)=r2;
  end
 end
 c(:,ind)=r1';
 if ~isempty(out); 
  out1=[out.slave0;out1(:)]; out=[out.c0;c];
 else; out=c;if ~isempty(out1); out1=out1(:);end
 end
 

% ------------------------------------------------------------------------
elseif comstr(Cam,'geomdilat')

 node=varargin{carg};carg=carg+1;
 if isfield(node,'Node'); node=node.Node;end

 if nargin<carg; nu=0.3; else nu=varargin{carg};carg=carg+1; end

 out=struct('DOF',[node(:,1)+0.01;node(:,1)+0.02;node(:,1)+0.03], ...
            'def',[node(:,5);-nu*node(:,6);-nu*node(:,7)],'lab',{{'Dilatation'}}  ); 

% ------------------------------------------------------------------------
elseif comstr(Cam,'geomrb')

 model=varargin{carg};carg=carg+1;
 if ~isfield(model,'Node') error('You must provide a model for GEOMRB');end

 [node,bas,NNode]=feutil('getnodebas',model);
 i1=unique(node(:,1));i1=i1(:)'; 
 rdof=[i1+.01;i1+.02;i1+.03];rdof=rdof(:);

 r2=zeros(length(rdof),6); r2(1:3:end,1)=1; r2(2:3:end,2)=1;r2(3:3:end,3)=1;
 i1=fe_c(rdof,node(:,1)+.01,'ind');
 r2(i1+1,4)=-node(:,7);r2(i1+2,4)=node(:,6);
 r2(i1+2,5)=-node(:,5);r2(i1+0,5)=node(:,7);
 r2(i1+0,6)=-node(:,5);r2(i1+1,6)=node(:,5);
 out=struct('def',r2,'DOF',rdof,'name','RB modes');

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

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
% feutil('GetDof',model)
if comstr(Cam,'dof') [CAM,Cam]=comstr(CAM,4);

model=varargin{carg};carg=carg+1;
ModelStack={};
if isnumeric(model)
 if finite(model(1))&size(model,2)==7 % Node elt
  carg=carg-1;
  [carg,FEnode,FEelt,FEel0,ModelStack]=get_nodeelt(varargin,carg,ModelStack);
   model=struct('Node',FEnode,'Elt',FEelt);carg=carg+1;
 else % backward compatible elt only
   model=struct('Node',[],'Elt',model,'Stack',{{}});
 end
elseif isfield(model,'Elt')
 carg=carg-1;
 [carg,FEnode,FEelt,FEel0,ModelStack]=get_nodeelt(varargin,carg,ModelStack);
end

if ~isempty(Cam) % element selection
 [out,model.Elt]=feutil(['findelt' CAM],FEnode,FEelt,FEel0,varargin{carg:end});
end
if ~isfield(model,'Stack');model.Stack=ModelStack;end
[eltid,model.Elt]=feutil('eltidfix;',model.Elt);
[EGroup,nGroup]=getegroup(model.Elt);

% determination of DOFs in the unrestrained model
% i3 generic node set, i4 generic dof set, out2 word count

ndof=[];edof=[]; out2=0; nd=[];

for jGroup=1:nGroup
  [ElemF,i1]=feutil('getelemf',model.Elt(EGroup(jGroup),:),jGroup);EGID=i1(1);
  cEGI = EGroup(jGroup)+1:EGroup(jGroup+1)-1;
  try; i4=feval(ElemF,'dofcall');catch;i4=[];end
  if EGID<0  % display group
  elseif ischar(i4); 
  % there is a dof call strategy wich is assumed to return 
  % i2 (nodal DOF*100), i3 element DOF*1000
  % beware this must increment out2 (word count) too
   eval(i4); ndof=[ndof;i2(:)];edof=[edof;abs(i3(:))];
  elseif EGID>=0 % not a display group
   [i3,SEopt]=fe_super('node',ElemF,model.Stack);
   i4=fe_super('dof',ElemF,model.Stack);
   i4=i4(:)*100;i3=i3(:);
   if size(SEopt,2)<3|SEopt(1,3)==0 i6=length(i4)^2;else i6=SEopt(1,3);end
   if     SEopt(1,1)==2 out2 = out2 + length(cEGI)*i6;  % generic
   elseif SEopt(1,1)==1 out2 = out2 + i6; end % unique

   i1=find(i4<0); if ~isempty(i1) % element DOFs
      if SEopt(1,1)==2 & ~isempty(cEGI) % generic
        i2=eltid(cEGI);i2=i2(:)'; 
	i2=i2(ones(size(i1)),:)*1e3-i4(i1,ones(size(cEGI)))*10-1000;
      else i2=-i4(i1)*10+(999000+EGID*1000);end;% unique
      edof=[edof;i2(:)];
   end
   i1=find(i4>0);
   if (~isempty(i1)&strcmp(ElemF,'rigid'))|strcmp(ElemF,'celas')
   % master and slave DOFs in rigid & celas
     i2=[];
     [i3,i8,i7]=unique(model.Elt(cEGI,3:4),'rows');
     for j1=1:size(i3,1)
      i8=cEGI(find(i7==j1));
      if i3(j1,1)>0
       i4=[1:6];i5=[abs(sprintf('%i',i3(j1,1)))-48];
      elseif i3(j1,1)<0 & i3(j1,2)==0
       i4=[abs(sprintf('%i',-i3(j1,1)))-48];i5=i4;
      elseif i3(j1,1)<0 
       i4=[abs(sprintf('%i',abs(i3(j1,1))))-48];
       i5=[abs(sprintf('%i',abs(i3(j1,2))))-48];
      elseif i3(j1)==0; 
       sdtw('_nb','Missing DOFs for %s (elements %i:%i)',ElemF,cEGI([1 end]));
       i4=[];i5=[];
      end
      if ~isempty(i4);
       i4=model.Elt(i8,ones(1,length(i4)))*100+i4(ones(length(i8),1),:);
       i5=model.Elt(i8,2*ones(1,length(i5)))*100+i5(ones(length(i8),1),:);
       i2=[i2;i4(:);i5(:)];
      end
     end
     i2=unique(i2(find(i2>100))); ndof=[ndof;i2];
   elseif ~isempty(i1) % nodal DOFs
      if SEopt(1,1)==2 & ~isempty(cEGI) % generic
	i3=i3(:);nind(i3,1)=[1:length(i3)]';
        i4=[nind(fix(i4(i1)/100)) round(rem(i4(i1),100))]';
        if max(i3)>size(model.Elt,2);model.Elt(1,max(i3))=0;end
	i2 = model.Elt(cEGI,i3)*100;
        if strcmp(ElemF,'mass1') 
         i2=i2(:,i4(1,:)).*(model.Elt(cEGI,2:7)~=0)+i4(2*ones(size(i2,1),1),:);
        else
           i2=i2(:,i4(1,:))+i4(2*ones(size(i2,1),1),:);
        end
      else i2=round(i4(i1)); end
      i2=unique(i2(:)); ndof=[ndof;i2];
   end
  end % of EGID>=0
end
ndof=unique(round(ndof))/100; ndof=ndof(find(ndof>1));
out = [ndof; -unique(edof)/1000];
out1=[];

if carg<=nargin&~isempty(varargin{carg})
    if carg+1<=nargin&~isempty(varargin{carg+1})
      r1=varargin{carg};i1=fe_c(r1,out,'ind');
      out1=varargin{carg+1};out1=out1(i1,:);out=r1(i1);
    else if ~isempty(out) out=fe_c(out,varargin{carg},'dof'); end;out1=[];
    end    
end

if nargout>2 out3=model.Elt;out4=eltid;end

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
%	elt=feutil('getedge[Patch,Line]')
elseif comstr(Cam,'edge') [CAM,Cam]=comstr(CAM,5);

[carg,node,elt,el0,ModelStack]=get_nodeelt(varargin,carg,ModelStack);

if comstr(Cam,'patch') % 'getEdgePatch' - - - - - - - - - - - - - - -
[CAM,Cam] = comstr(CAM,6);

% SelEdgePatch replaces the selected group(s0 by a quad4 group containing
% edge patches
 i4=[];
 [EGroup,nGroup]=getegroup(elt); i6=[];
 % note : this will not work if two face sets have the same number of nodes

 [i1,elt]=feutil('eltidfix;',elt);
 FaceStack={};
 for jGroup = 1:nGroup
   
   [ElemF,i1,ElemP]= feutil('getelemf',elt(EGroup(jGroup),:),jGroup);
   r1=fe_super('face',ElemF);cEGI = EGroup(jGroup)+1:EGroup(jGroup+1)-1;
   i1=fe_super('prop',ElemF); 
   if isempty(r1)
   else
    % decompose in face subsets for variable node numbers
    r2={};while ~isempty(r1)&size(r1,2)>1;
     i2=find(r1(:,end-1)~=r1(:,end));
     if ~isempty(i2); r2{end+1}=r1(i2,:);r1(i2,:)=[];end;
     r1(:,end)=[];
    end

    if ~isempty(r1);r2{end+1}=r1;end
    for j1=1:length(r2);
     r1=r2{j1}; % loop on number of node per face
     if length(FaceStack)<size(r1,2); FaceStack{2,size(r1,2)}=[];end
     i5=ones(size(r1,1)*length(cEGI),1)*jGroup;
     if i1(1); i2=elt(cEGI,i1(1)*ones(1,size(r1,1)))';i5(:,2)=i2(:); end
     if i1(2); i2=elt(cEGI,i1(2)*ones(1,size(r1,1)))';i5(:,3)=i2(:); end
     if i1(3)&size(elt,2)>=i1(3)
            i2=elt(cEGI,i1(3)*ones(1,size(r1,1)))';i5(:,4)=i2(:);
     end
     i2=reshape(elt(cEGI,r1')',size(r1,2),size(r1,1)*length(cEGI));

     FaceStack{1,size(r1,2)}=[FaceStack{1,size(r1,2)};i2'];
     FaceStack{2,size(r1,2)}=[FaceStack{2,size(r1,2)};i5];
    end
   end
 end
 elt=[];
 for j1=1:size(FaceStack,2)
  if ~isempty(FaceStack{1,j1})
   i4=FaceStack{1,j1};
   [i2,i1]=sortrows(sort(i4,2)); % eliminate dupicated faces
   i3=find(~any(diff(i2),2));i1(i3)=0;i1(i3+1)=0;i1=i1(find(i1));
   switch size(i4,2)
   case 3
     elt(end+[1:length(i1)+1],1:6) = [Inf abs('tria3') ;
           i4(i1,:)  FaceStack{2,j1}(i1,2:4)];
   case 4
     elt(end+[1:length(i1)+1],1:7) = [Inf abs('quad4') 0 ;
           i4(i1,:)  FaceStack{2,j1}(i1,2:4)];
   case 6
     elt(end+[1:length(i1)+1],1:9) = [Inf abs('tria6') 0 0 0;
           i4(i1,:)  FaceStack{2,j1}(i1,2:4)];
   case 8
     el1=[Inf abs('quadb') zeros(1,5) ;
           i4(i1,1:8)  FaceStack{2,j1}(i1,2:4)];
      % deal with degenerated quadb
     ind=find(el1(:,6)==el1(:,7) & el1(:,6)==el1(:,8) ...
         & el1(:,6) & el1(:,7) & el1(:,8));
     i2=el1(ind,[1:6 9:end]); 
     el1(ind,:)=[];
     el1(end+[1:length(ind)+1],1:size(i2,2)) = ...
             [Inf abs('tria6') zeros(1,size(i2,2)-6);i2];
     elt(end+[1:size(el1,1)],1:size(el1,2)) = el1;
   otherwise error('Not supported get face');
   end
  end
 end



% 'getEdgeLine[ProId,MatId,All,In]'  - - - - - - - - - - - - - - - - - - - -
elseif comstr(Cam,'line')  

 [CAM,Cam] = comstr(CAM,5);
 % GetEdgeLine replaces the selected group by a beam1 group containing 
 % segments at the edge of the of the current group

 if comstr(Cam,'g') [CAM,Cam]=comstr(CAM,'group','%c');    opt=1;
 elseif comstr(Cam,'p') [CAM,Cam]=comstr(CAM,'proid','%c');opt=2;
 elseif comstr(Cam,'m') [CAM,Cam]=comstr(CAM,'matid','%c');opt=3;
 elseif comstr(Cam,'a') [CAM,Cam]=comstr(CAM,'all','%c');opt=5;
 elseif comstr(Cam,'i') opt=4; ind=varargin{carg};carg=carg+1;
     [a,elt]=feutil('eltidfix;',elt);
 else opt=0; end

 i4=[];[EGroup,nGroup]=getegroup(elt); i6=[];
 for jGroup = 1:nGroup
   ElemF= feutil('getelemf',elt(EGroup(jGroup),:),jGroup);
   cEGI = EGroup(jGroup)+1:EGroup(jGroup+1)-1;LD2(jGroup,1)=0;
   [NodeLine,SEopt]=fe_super('edge',ElemF);NodeLine=NodeLine';
     
   if ~isempty(cEGI) i3=elt(cEGI,NodeLine(:))';
      % case of a single superelement
   else  [i3]=fe_super('node',ElemF);i3=i3(NodeLine);     end
   i5=ones(size(NodeLine,2)*size(i3,2),1)*jGroup;
   i1=fe_super('prop',ElemF);
   if i1(1)&size(elt,2)>=i1(1)
    i2=elt(cEGI,i1(1)*ones(1,size(NodeLine,2)))';i5(:,2)=i2(:); 
   end
   if i1(2)&size(elt,2)>=i1(2) 
     i2=elt(cEGI,i1(2)*ones(1,size(NodeLine,2)))';i5(:,3)=i2(:); 
   end
   if i1(3)&size(elt,2)>=i1(3)
            i2=elt(cEGI,i1(3)*ones(1,size(NodeLine,2)))';i5(:,4)=i2(:);
   end

   if ~all(i3) error('Some nodes are equal to zero'); end
   i3=reshape(i3,size(NodeLine,1),size(NodeLine,2)*length(cEGI))';
   i4(end+1:end+size(i3,1),1:size(i3,2))=i3;
   i6(end+1:end+size(i5,1),1:size(i5,2))=i5;
 end % loop on groups

 if opt==5 % Single instance of all Edges
  [i3,i1]=unique(sort(i4,2),'rows');
  i4=i4(i1,:);i6=i6(i1,:);

 elseif opt==4 % EdgesInNodes

   nind=zeros(max(node(:,1)),1);nind(ind)=1:length(ind);
   i3=find(all(nind(i4),2)); 
   i4=i4(i3,:); i6=i6(i3,:); 
   [i4,i1]=sortrows(sort(i4,2));  i6=i6(i1,:);
   i2=find(~any(diff(i4),2)); i6(i2,1)=i6(i2,4); % EltNumber in col 1
   nind=[1:size(i4,1) 0]; nind(i2)=0; j1=1;
   while any(nind)
    i3=i2+j1; if ~isempty(i3) i5=find(nind(i3)); else i5=[];end
    if ~isempty(i5)   j1=j1+1;i6(i2(i5),j1)=i6(i3(i5),4); nind(i3(i5))=0;end
   end
   i6=[i2(:) i6(i2,1:j1)]; i4=i4(i2,:);

 elseif opt % some gpm selection

  ind=find(sparse(i6(:,opt)+1,1,1))-1; % selection values
  i3=[]; i5=[];
  for j1=1:length(ind) % edge of each group/mat/prop
    in2=find(i6(:,opt)==ind(j1));
    [i1,i2]=single_line(i4(in2,:),i6(in2,:));
    i3=[i3;i1];i5=[i5;i2];
  end
  % now eliminate duplications & create new number of interface boundaries
  [i3,i1]=sortrows(i3);   i2=find(~any(diff(i3),2));
  i5(i1(i2),opt)=sort([i5(i1(i2),opt) i5(i1(i2+1),opt)],2)* ...
          [1;10^ceil(log10(length(ind)))];
  in2=1:size(i3,1);in2(i2+1)=0;i2=find(in2);
  i3=i3(i2,:);i5=i5(i1(i2),:);
  i5(:,2)=i5(:,opt); % new numbers as materials
  i4=i3;i6=i5;
 else % standard selection of edges
  [i4,i6]=single_line(i4,i6);
 end
 
 if isempty(i4)
   warning('GetEdge: selected group has no edges'); elt=[];
 else 
    if size(i4,2)<3 i4(:,3)=0;end
    i1=find(i4(:,3)==0&i4(:,1));elt=[];
    if ~isempty(i1)
      elt = [Inf abs('beam1')];
      elt(2:length(i1)+1,1:1+size(i6,2))=[i4(i1,1:2) i6(i1,2:end)];
      i4(i1,:)=0;
    end
    i1=find(all(i4,2));
    if ~isempty(i1)
      elt(end+1,1:6) = [Inf abs('beam3')];
      elt(2:length(i1)+1,1:2+size(i6,2))=[i4(i1,[1:3]) i6(i1,2:end)];
      i4(i1,:)=0;
    end

    if any(any(i4))
     warning('Edges with 2 or 3 nodes are the only supported');
    end

 end

 elt(find(finite(elt(:,1))&elt(:,1)==elt(:,2)),:)=[];%eliminate nodes

else error('Not a known SelEdge command');
end % of SelEdgeLine

[CAM,Cam]=comstr(CAM,1); 
if ~isempty(CAM)&CAM(1)=='{'&CAM(end)=='}' [CAM,Cam]=comstr(CAM(2:end-1),1);end

out=elt; out1=CAM;

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
%	[ElemF,opt]= feutil('getelemf',elt(EGroup(jGroup),:),jGroup)
elseif comstr(Cam,'elemf') [CAM,Cam]=comstr(CAM,6);

   st=varargin{carg};carg=carg+1; 
   i1 = [min([find(~st) find(st==32)]) length(st)+1];
   out = setstr(st(2:i1(1)-1));
   if nargout>1 
    out1 = st(i1(1)+1:length(st)); 
    if carg<=nargin r1=varargin{carg};carg=carg+1; else r1=0; end
    if isempty(out1) out1=r1; elseif out1(1)==0 out1=r1(1); end
   end
   if nargout>2    out2=fe_super('parent',out);  end

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 % this should be maintained for edge selection
% feutil('GetLine',node,elt)
elseif comstr(Cam,'line') [CAM,Cam]=comstr(CAM,5);

[carg,node,elt,ModelStack]=get_nodeelt(varargin,carg,ModelStack);

if carg>nargin nInf=0;else nInf=varargin{carg};carg=carg+1;end
if ~isempty(elt)
  EGroup = [find(elt(:,1)==Inf)' size(elt,1)+1];nGroup = length(EGroup)-1;
else	nGroup=0;EGroup=[];elt([1 83])=[1 0];  end

if nGroup==0 % trace line plot - - - - - - - - - - - - - - - - - - - -

   for j1 = 1:size(elt,1);
        i1=elt(j1,[83:82+elt(j1,1)]);
	LD2(j1,1:length(i1)+1) = [length(i1)+1 i1];
   end

elseif nGroup>0 % element plot - - - - - - - - - - - - - - - - - - - -

LD2(nGroup,1)=0; out=struct('Group',[]);out.Group={};out.Line={};

for jGroup = 1:nGroup
   ElemF= feutil('getelemf',elt(EGroup(jGroup),:),jGroup);
   cEGI = EGroup(jGroup)+1:EGroup(jGroup+1)-1;
   [NodeLine,SEopt]=fe_super('line',ElemF);NodeLine=NodeLine(:)';

   if (~isempty(NodeLine)&~isempty(cEGI)) | SEopt(1,1)==1

     if length(NodeLine)>1
      NodeLine=[NodeLine(1:length(NodeLine)-1);NodeLine(2:length(NodeLine))];
     else NodeLine=[NodeLine;NodeLine];end
     i1 = find(NodeLine(1,:)&NodeLine(2,:));NodeLine=NodeLine(:,i1);
     NodeLine=NodeLine(:);
     
     if ~isempty(cEGI) i3=elt(cEGI,NodeLine)'; % case of a unique superelement
     else  [i3]=fe_super('node',ElemF);i3=i3(NodeLine,1); end

     if ~isempty(i3)

        i3=reshape(i3,2,prod(size(i3))/2);
        i3=sort(i3(:,find(all(i3))));% ignoring nodes set to zero
	[i1,i2]=find(sparse(i3(1,:),i3(2,:),i3(1,:)));
        i1=[i1 i2 zeros(size(i2))]'; i1=i1(:);
        LD2(jGroup,1:length(i1)) = [length(i1) i1(1:length(i1)-1)'];
        out.Group{end+1}=sprintf('Group %i',jGroup);
        out.Line{end+1}=i1(1:length(i1)-1)';
     elseif isempty(i3)        LD2(jGroup,1) = [0];
     end % all(i3) the nodes are in FEnode
else disp(sprintf('No ''line''  plot for group %i (%s)',jGroup,ElemF));end

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
end;end % jGroup and type of nGroup

if ~isempty(strfind(Cam,'struct'))
 return;
elseif nargout>2
 i1 =LD2(:,2:end);i1=i1(:);out2=zeros(max(i1)+1,1);out2(i1+1)=1;
 out2=find(out2)-1;out2=out2(2:end);
end
if ~isempty(node) % if node is defined give in indices rather than numbers
 nInf=size(node,1);
 if finite(node(nInf,5))
   node=[node;max(node(:,1))+1 0 0 0 Inf Inf Inf];nInf=nInf+1;
 else node(nInf,1)=max(node(1:nInf-1,1))+1;end
 NNode=sparse(node(:,1),1,1:size(node,1));
 for j1=1:size(LD2,1); if LD2(j1,1)>0
  i1 =  LD2(j1,2:LD2(j1,1)); i2=find(i1==0); i1(i2)=i1(i2)+node(nInf,1);
  LD2(j1,1:length(i1)+1)=[length(i1)+1 reshape(full(NNode(i1)),1,length(i1))];
 end;end
 if nargout>2 out2=[full(NNode(out2))';nInf];end
end
out = LD2; out1=nInf;

% end of GetLine
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
% feutil('GetPatch',node,elt)
elseif comstr(Cam,'patch') [CAM,Cam]=comstr(CAM,6);

model=[];
[carg,node,elt,el0,ModelStack]=get_nodeelt(varargin,carg,ModelStack);
if isempty(model); model=struct('Stack',{ModelStack});end

if carg>nargin nInf=size(node,1);else nInf=varargin{carg};carg=carg+1;end
[nElt,i1]=size(elt);nInf=size(node,1);
if nElt>0; [EGroup,nGroup]=getegroup(elt);% elements have been defined
else	nGroup=0;EGroup=[];elt([1 83])=[1 0];  
end

if comstr(Cam,'new') % 'GetPatchNew returns a feplot selection

 out=struct('opt',[0 0 0 0 0],'selelt','', ...
     'vert0',[],'cna',[],'fs',[],'ifs',[],'f1',[],'if1',[],'f2',[],'if2',[],...
     'fsProp',[],'f2Prop',[],'f1Prop',[]);
 out.fsProp={'FaceColor',[1 1 1],'edgecolor',[0 0 0],'facelighting','phong'};
 out.f2Prop={'FaceColor','none','edgecolor',[0 0 0],'linewidth',2};
 out.f1Prop={'FaceColor','none','edgecolor',[0 0 0],'marker','o'};

end
if nGroup==0 % trace line plot - - - - - -
 i2 =elt; elt=[]; for j1=1:size(i2,1);
  i1=[i2(j1,82+[1:i2(j1,1)-1]);i2(j1,82+[2:i2(j1,1)])];
  i1=i1(:,find(~any(i1==nInf)&~any(i1==0)))';
  elt=[elt;Inf abs('beam1'); i1 zeros(size(i1,1),4)];
 end
 EGroup = [find(elt(:,1)==Inf)' size(elt,1)+1];nGroup = length(EGroup)-1;
end

LD2=[];i7=[];

for jGroup = 1:nGroup
   [ElemF,i1]= feutil('getelemf',elt(EGroup(jGroup),:),jGroup);
   cEGI = EGroup(jGroup)+1:EGroup(jGroup+1)-1;
   EGID = i1(1); i7(end+1)=EGID;
   [NodeSurf,SEopt]=fe_super('patch',ElemF,model);

% new version of PATCH matrix corresponds to patch faces property
if SEopt(1,1)==3  % locally built patch
  i2=NodeSurf;i7 = ones(1,size(i2,1)); % needs revising
elseif ~isempty(NodeSurf)&(SEopt(1,1)==1|~isempty(cEGI))  i2=[];
  if SEopt(1,1)==2 % generic superelement
   i2=[]; i7=[];for j1=1:size(NodeSurf,1) 
       i1 = NodeSurf(j1,:);
       i1 = reshape(elt(cEGI,i1),length(cEGI),length(i1));
       i2(end+[1:size(i1,1)],1:size(i1,2))=i1; i7 = [i7;cEGI(:)];
   end % of loop on Patches linked to the element
  else % single superelement
    if ~isempty(setdiff(NodeSurf(:),node(:,1)));  % XXX FESUPER CHECK
     i2=reshape(node(NodeSurf,1),size(NodeSurf,1),size(NodeSurf,2));
    else; i2=NodeSurf;
    end
    i7 = ones(1,size(i2,1));
  end
  % for volumes elimination of duplicated patches within same group
  if size(i2,2)>2 & size(NodeSurf,1)>1
   i1 = sort(i2,2); [i1,i3]=sortrows(i1);
   i4=find(~any(diff(i1),2));
   i5=ones(size(i2,1),1);i5(i3(i4))=0;i5(i3(i4+1))=0; i5=find(i5);
   i2=i2(i5,:);i7=i7(i5);
  end
elseif isfield(NodeSurf,'fs')
else 
   disp(sprintf('No ''patch'' plot for group %i (%s)',jGroup,ElemF));
   i2=[]; i7=[];
end

 if comstr(Cam,'new') % returns a feplot selection
   if     size(i2,2)==1  % Points
    out.f1(end+[1:size(i2,1)],1)=i2;
    out.if1(end+[1:length(i7)],1)=i7;
   elseif size(i2,2)==2 % Lines
    i3=find(any(i2==0,2)|i2(:,2)==i2(:,1));
    if ~isempty(i3) % Show degenerate as points
      out.f1(end+[1:length(i3)],1)=i2(i3,1);
      out.if1(end+[1:length(i3)],1)=i7(i3);
      i2(i3,:)=[]; i7(i3,:)=[];
    end
    if ~isempty(i2) % if some non degenerates
     out.f2(end+[1:size(i2,1)],1:2)=i2; out.if2(end+[1:length(i7)])=i7;
    end
   elseif size(i2,2)>2 % standard faces
     if ~isempty(out.fs) & size(out.fs,2)<size(i2,2) 
       out.fs(:,end+1:size(i2,2))= ...
          out.fs(:,end)*ones(1,size(i2,2)-size(out.fs,2));
     end
     % filling of zeros xxx needs a systematic retest
     i1=find(~i2); while ~isempty(i1)
       if any(i1<size(i2,1)) error('First node in face cannot be 0');end
       i2(i1)=i2(i1-size(i2,1)); i1=find(~i2); 
     end
     if size(i2,2)<size(out.fs,2)
       i2(:,end+1:size(out.fs,2))=i2(:,end)*ones(1,size(out.fs,2)-size(i2,2));
     end
     out.fs(end+[1:size(i2,1)],1:size(i2,2))=i2;
     out.ifs(end+[1:length(i7)])=i7;
   end  


 else % need something for old/feplot compatibility
   i2(:,5)=i7;i2(:,6)=jGroup;
   LD2(end+[1:size(i2,1)],1:size(i2,2))=[i2];
 end
end % jGroup  - - - - -

if comstr(Cam,'new') % returns a feplot selection
  NNode=sparse(node(:,1)+1,1,1:size(node,1));

  % get nodes used by the faces
  i2=[out.fs(:);out.f1(:);out.f2(:)];
  out.Node=find(sparse(i2+1,1,i2))-1;
  if max(out.Node)>length(NNode); NNode(max(out.Node)+1)=0; end
  i2 =NNode(out.Node+1);
  i1=find(i2<1);if ~isempty(i1)
    fprintf('%i %i %i %i %i %i\n',out.Node(i1));
    error('NodId in .Elt but not in .Node');
  end
  out.vert0=node(NNode(out.Node+1),5:7);
  % renumber faces
  NNode=sparse(out.Node,1,1:size(out.Node,1));
  %NNode(out.Node)=[1:length(out.Node)]';

  if ~isempty(out.fs) 
   out.opt(1,3)=1;i1=find(out.fs);out.fs(i1)=NNode(out.fs(i1));
   if ~isempty(find(~any(out.fs,2)))
     error('Some faces have zero nodes');
   end
  end
  if size(out.fs,1)==size(out.vert0,1) % duplicate face for color control
    out.fs=out.fs([1:end end],:); out.ifs=out.ifs([1:end end]);
  end

  if ~isempty(out.f1)
   out.opt(1,5)=1;i1=find(out.f1);  out.f1(i1)=NNode(out.f1(i1));
   if ~isempty(find(~any(out.fs,2)))
     error('Some point faces have zero nodes');
   end
  end
  if ~isempty(out.f2) 
   out.opt(1,4)=1;i1=find(out.f2);  out.f2(i1)=NNode(out.f2(i1));
   if ~isempty(find(~any(out.f2,2)))
     error('Some 2 node faces have zero nodes');
   end
   i1=find(~out.f2(:,2));if ~isempty(i1) out.f2(i1,2)=out.f2(i1,1);end
   i1=find(~out.f2(:,1));if ~isempty(i1) out.f2(i1,1)=out.f2(i1,2);end
  end
  try;if isfield(varargin{2},'IndWithHeaders')
    i2=varargin{2}.IndWithHeaders;
    if isfield(out,'ifs') out.ifs=i2(out.ifs); end 
    if isfield(out,'if1') out.if1=i2(out.if1); end 
    if isfield(out,'if2') out.if2=i2(out.if2); end 
  end;end
else  % not new
  if ~isempty(node) % if node is defined in indices rather than numbers
   nInf=size(node,1); if finite(node(nInf,5)) nInf=nInf+1; end
   NNode(node(:,1)+1)=[1:length(node(:,1))]';
   i2=find(LD2(:,1:4));LD2(i2)=NNode(LD2(i2)+1);
  end
  i2=cumsum(~finite(elt(:,1)));
  out=[1 1 size(elt,1) zeros(1,size(LD2,2)-3);LD2];
end


% end of GetPatch
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
% feutil('GetEG',node,NNode,elt_head,elt_elt)
elseif comstr(Cam,'eg') [CAM,Cam]=comstr(CAM,3);

node  = varargin{carg};carg=carg+1;
NNode = varargin{carg};carg=carg+1;
elt   = varargin{carg};carg=carg+1;
[ElemF,EGID]= feutil('getelemf',elt(1,:),1);cEGI=2:size(elt,1);

% nodes linked to group

[i2,SEopt]=fe_super('node',ElemF); if SEopt(1,1)==2 i2 = elt(cEGI,i2);end
i2=i2(:); out.Node = node(NNode(find(sparse(i2+1,1,i2))),:);

NNode=[]; NNode(out.Node(:,1)+1)=1:size(out.Node,1); nInf=size(out.Node,1)+1; 
out.Elt=elt;

% Basic Patch representation of Group - - - - - - - - - - - - - - - - - -

LD2=[1 1 size(elt,1)];
[Faces,SEopt]=fe_super('patch',ElemF);
if isempty(Faces) % in case there is no patch definition 
  [NodeLine,SEopt]=fe_super('line',ElemF);NodeLine=NodeLine(:)';
  if length(NodeLine)>1
  NodeLine=[NodeLine(1:length(NodeLine)-1);NodeLine(2:length(NodeLine))];
  else NodeLine=[NodeLine;NodeLine];end
  i1 = find(NodeLine(1,:)&NodeLine(2,:));Faces=NodeLine(:,i1)';
end

% Patch matrix : Object vertex matrix 
if ~isempty(Faces)&(SEopt(1,1)==1|~isempty(cEGI))

  if SEopt(1,1)==2 % generic element or superelement
   i2=[]; i7=[];for j1=1:size(Faces,1) 
       i1 = Faces(j1,:);
       i1 = reshape(elt(cEGI,i1),length(cEGI),length(i1));
       i2(end+[1:size(i1,1)],1:size(i1,2))=i1; i7 = [i7 cEGI];
   end % of loop on Patches linked to the element
  else % single superelement
    i2=fe_super('node',ElemF); i2=i2(:,1); i1=find(Faces);
    Faces(i1)=i2(Faces(i1));  i2=Faces; i7 = ones(1,size(i2,1));
  end
  % elimination of duplicated patches within same group
  if  size(i2,2)>2  & ~any(strcmp(fe_super('parent',ElemF), ...
    {'beam3','tria3','tria6','quad4','quadb'}))
   i1 = sort(i2,2); [i1,i3]=sortrows(i1);
   i4=find(~any(diff(i1),2));i5=ones(size(i2,1),1);i5(i3(i4))=0;i5(i3(i4+1))=0;
   i5=find(i5);i2=i2(i5,:);i7=i7(i5);
  elseif size(i2,2)==1 i2(:,2)=0;
  elseif size(i2,2)==2 i2(:,3)=0;
  end
end
i1=find(i2);i2(i1)=NNode(i2(i1)+1);  % switch to node indices

% treat faces with nodes at zero here

% output the various face matrices
i1 = sum((finite(i2)&i2)');

i3 = find(i1>2);  if ~isempty(i3) out.ps=i2(i3,:); out.es=i7(i3);
out.opt(1,1)=1; end
i3 = find(i1==2); if ~isempty(i3) out.p2=i2(i3,1:2); out.e2 = i7(i3); 
out.opt(1,2)=1;end
i3 = find(i1==1); if ~isempty(i3) out.p1=i2(i3,1); out.e1 = i7(i3); 
if any(out.p1==0) error(1); end
out.opt(1,3)=1;else out.opt(1,3)=0; end
if ~any(i1)
  warning(sprintf('Group %i (%s) is not displayed',varargin{5},ElemF));
end
out.opt(1,5)=0;

% end of GetEG
% 'GetNormal [,Map][Elt,Node]' - - - - - - - - - - - - - - - - - - - - - -
elseif comstr(Cam,'normal') [CAM,Cam]=comstr(Cam,7);

  if strfind(Cam,'map'); 
   CAM(strfind(Cam,'map')+[0:2])=''; [CAM,Cam]=comstr(CAM,1); st='map'; 
  else st='';end

  % allow for selection
  if isfield(varargin{carg},'Elt')&carg<nargin&ischar(varargin{carg+1}) 
   model=varargin{carg};carg=carg+1; node=model.Node;
    elt=feutil(['selelt' varargin{carg}],model);
  else
    [carg,node,elt,ModelStack]=get_nodeelt(varargin,carg,ModelStack);
  end
  [EGroup,nGroup]=getegroup(elt);
  NNode=[];NNode(node(:,1)+1)=1:size(node,1);
  out=zeros(size(elt,1),3); out1=zeros(size(elt,1),3);
  out2=zeros(size(elt,1),1); out3=zeros(size(elt,1),3);

  if comstr(Cam,'node') r3=spalloc(size(elt,1),max(node(:,1)),1);
  else r3=[]; end

  for jGroup = 1:nGroup %loop on groups
     [ElemF,i1,ElemP]=feutil('getelemf',elt(EGroup(jGroup),:),jGroup);
     i1=fe_super('node',ElemP);
     cEGI = [EGroup(jGroup)+1:EGroup(jGroup+1)-1]';
     if comstr(ElemP,'quad4')|comstr(ElemP,'quadb') | ...
        comstr(ElemP,'tria3')|comstr(ElemP,'tria6')|comstr(ElemP,'quad9')
      for jElt=1:length(cEGI) %loop on elements of group
        i2=NNode(elt(cEGI(jElt),i1)+1);
        r1 = node(i2(find(i2)),5:7);
        r2 = mean(r1,1); out1(cEGI(jElt),1:3)=r2;
        r1=r1-r2(ones(size(r1,1),1),:);
        [bas,r1,out2(cEGI(jElt),1)]= basis(r1);
        out (cEGI(jElt),1:3)=bas(:,3)'; % normal
        out3(cEGI(jElt),1:3)=bas(:,1)'; % x direction
      end % loop on elements of group
     % orient 2-d elements by pointing along axis
     %elseif comstr(ElemP,'beam3') 
     % [i2,i3]=min(std(node(:,5:7)));r1=zeros(1,3);r1(i2)=1; % plane normal
     elseif any(strcmp(ElemP,{'beam1','bar1'}))

      z0=std(node(:,5:7));[z0,ind]=min(z0);z0=full(sparse(1,ind,1,1,3));
      for jElt=1:length(cEGI) %loop on elements of group
        i2=NNode(elt(cEGI(jElt),i1)+1);
        r1 = node(i2(find(i2)),5:7);
        if size(r1,1)>1;
          r2 = mean(r1,1); out1(cEGI(jElt),1:3)=r2;
          r1=diff(r1,[],1); out2(cEGI(jElt))=norm(r1);
          p=sp_util('basis',r1,z0); %basis(-diff(r1,[],1),[ 0 0 0],1);
          out(cEGI(jElt),1:3)=p(:,3)';
        end
      end % loop on elements of group
     else
      for jElt=1:length(cEGI) %loop on elements of group
        i2=NNode(elt(cEGI(jElt),i1)+1);
        r1 = node(i2(find(i2)),5:7);
        r2 = mean(r1,1); out1(cEGI(jElt),1:3)=r2;
      end % loop on elements of group
     end
     % centroid
     if ~isempty(r3) r3=r3+sparse(cEGI(:,ones(size(i1))),elt(cEGI,i1), ...
         out2(cEGI(:,ones(size(i1)))),size(elt,1),max(node(:,1)));
     end
  end % loop on groups

  if comstr(st,'map')  % New Normal map format

    if ~isempty(r3) % MAP at nodes

     i1=find(any(r3)); r3=r3(:,i1);   r3=r3*diag(sum(r3).^(-1));
     r3=out'*r3;r3=r3*diag(sparse(sum(r3.^2).^(-.5)));
     out=struct('normal',r3','ID',i1(:),'opt',2);

    else % MAP at elements
     eltid=feutil('eltidfix;',elt);
     i1=find(eltid & any(out,2));
     out=struct('normal',out(i1,:),'ID',eltid(i1),'vertex',out1(i1,:),'opt',1);
    end
   
  else  % old nor/cg outputs

    if ~isempty(r3) % determine normal on each node
     i1=find(any(r3)); r3=r3(:,i1);   r3=r3*diag(sum(r3).^(-1));
     r3=out'*r3;r3=r3*diag(sparse(sum(r3.^2).^(-.5)));
     out=[i1(:) i1(:) r3'];
    end

  end

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
% cg=feutil('GetCG',model) get element center
elseif comstr(Cam,'cg')

  % allow for selection
  if isfield(varargin{carg},'Elt')&carg<nargin&ischar(varargin{carg+1}) 
   model=varargin{carg};carg=carg+1; node=model.Node;
   elt=feutil(['selelt' varargin{carg}],model);
  else
    [carg,node,elt,ModelStack]=get_nodeelt(varargin,carg,ModelStack);
  end
  [EGroup,nGroup]=getegroup(elt);
  NNode=[];NNode(node(:,1)+1)=1:size(node,1);
  out=zeros(size(elt,1),1); 
  for jGroup = 1:nGroup %loop on groups
     [ElemF,i1,ElemP]=feutil('getelemf',elt(EGroup(jGroup),:),jGroup);
     i1=fe_super('node',ElemP);
     cEGI = [EGroup(jGroup)+1:EGroup(jGroup+1)-1]';
     for jElt=1:length(cEGI) %loop on elements of group
        i2=NNode(elt(cEGI(jElt),i1)+1);
        r1 = node(i2(find(i2)),5:7);out(cEGI(jElt),1:3)=mean(r1,1);
     end % loop on elements of group
  end % loop on groups

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
% [node,bas,NNode]=feutil('GetNodeBas',model) in global coordinate system
elseif comstr(Cam,'nodebas')

model=varargin{carg};carg=carg+1;
if isfield(model,'bas') &~isempty(model.bas)
  [out,out1]=basis(model.Node,model.bas);
  if nargout>3 % Return local to global transformation
    out3=basis('trans l',out1,out);
  end
else out=model.Node; out1=[]; out3=[];
end
if nargout>2 out2=sparse(out(:,1),1,1:size(out,1));end

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
% node=feutil('GetNode (sel)',model) in global coordinate system
elseif comstr(Cam,'node')

[out1,out]=feutil(['findnode' CAM(5:end)],varargin{2:end});

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
% pl=feutil('Getpl (MatId)',model)
elseif comstr(Cam,'pl');[CAM,Cam]=comstr(CAM,3);

i1=comstr(CAM,-1); model=varargin{carg};carg=carg+1;
if isempty(i1);i1=feutil('matid',model);i1=unique(i1);i1=i1(find(i1));end
out=model.pl(find(ismember(model.pl(:,1),i1)),:);

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
% il=feutil('Getil (ProId)',model)
elseif comstr(Cam,'il');[CAM,Cam]=comstr(CAM,3);

i1=comstr(CAM,-1); model=varargin{carg};carg=carg+1;
if isempty(i1);i1=feutil('proid',model);i1=unique(i1);i1=i1(find(i1));end
out=model.il(find(ismember(model.il(:,1),i1)),:);


% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
else sdtw('''Get%s'' unknown',CAM); end % subcommand selection - - - - - - - 

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
% [model,def]=feutil('DefLocalToGlobal',model,def) deformation in global
elseif comstr(Cam,'deflocaltoglobal')

model=varargin{carg};carg=carg+1;
def=varargin{carg};carg=carg+1;
if ~isfield(model,'DOF'); 
   model.DOF=feutil('getdof',model); 
end
[node,bas,NNode,r1]=feutil('getnodebas',model);
[r1,mdof]=basis('trans',bas,node,def.DOF);
i1=find(any(r1,2));r1=r1(i1,:);mdof=mdof(i1);
def.def=r1*def.def; def.DOF=mdof; 
model.Node(:,3)=0;

out=model; out1=def;


% ------------------------------------------------------------------------
elseif comstr(Cam,'info')  [CAM,Cam]=comstr(CAM,5);

if length(Cam)==0 % guess info about what
 r1=varargin{carg};
 if isstruct(r1)|(~isempty(r1)&~finite(r1(1))) 
   Cam=['elt' Cam];CAM=['elt' CAM];
 end
end

if comstr(Cam,'elt') [CAM,Cam]=comstr(CAM,4); % 'InfoElt' - - - - 

  elt = varargin{carg};carg=carg+1;
  if isstruct(elt)&isfield(elt,'Elt') elt=elt.Elt; end
  if isempty(elt)|ischar(elt) st=sprintf('\nModel %s is empty\n',CAM);
  else
   [EGroup,nGroup]=getegroup(elt);
   st=sprintf('\nModel %s with %i element(s) in %i group(s)\n', ...
        CAM,size(elt,1)-nGroup,nGroup);

   table={};
   for jGroup = 1:nGroup %loop on element groups
    [ElemF,i1]= feutil('getelemf',elt(EGroup(jGroup),:),jGroup);
    [i3,SEopt]=fe_super('node',ElemF); i4 =fe_super('prop',ElemF);
    i3 = EGroup(jGroup+1)-EGroup(jGroup)-1;
    table{jGroup,1}=sprintf('%i',jGroup);table{jGroup,6}='';
    table{jGroup,4}=ElemF;
    if SEopt(1,1)==2
     table{jGroup,3}=sprintf('%i',i3);
     if i1(1)<0
      if i1(1)==-1  table{jGroup,2}='test wire frame';
      elseif i1(1)==-2 table{jGroup,2}='test/FEM node links';
      elseif i1(1)==-3 table{jGroup,2}='rotation interpolation links';
      else             table{jGroup,2}=sprintf('EGID %i',i1(1)); end
     elseif i1(1)~=jGroup; table{jGroup,2}=sprintf('%i',i1(1));
     else table{jGroup,2}=''; end
     if i4(1)~=0 & i4(1)<=size(elt,2) % MatId
      i5 = elt(EGroup(jGroup)+1:EGroup(jGroup+1)-1,i4(1));
      i5=unique(fix(i5));
      table{jGroup,5}=sdtw('_clip 20 1','%i ',i5);
     end
     if i4(2)~=0 & i4(2)<=size(elt,2) % proid
      i5=EGroup(jGroup)+1:EGroup(jGroup+1)-1;
      if ~isempty(i5) i5 = elt(i5,i4(2));
       i5=unique(fix(i5));%i5 = find(sparse(fix(i5+1),1,i5+1))-1; 
       table{jGroup,6}=sdtw('_clip 20 1','%i ',i5);
      end
     end
    elseif SEopt(1,1)==1 % single superelement
     if i1(1)~=jGroup table{jGroup,2}=sprintf('%i',i1(1)); end
     if SEopt(1,2)==0 table{jGroup,5}='single SE';
     else
      table{jGroup,5}='single SE';
      table{jGroup,6}=sprintf('%i',SEopt(1,2));
     end
    else
     table{jGroup,5}='Element Not Found';
    end
   end % of jGroup loop
   if nargout<2 % text display

     st2='Group ';i1=ones(size(table,1),1);
     st1=[st2(i1,:) char(table{:,1}) st2(i1,end)];
     i2=size(st1,2);
     for j1=1:size(table,1)
      if ~isempty(table{j1,2})
       st2=sprintf('(%s) ',table{j1,2});
       st1(j1,i2+[1:length(st2)])=st2;
      end
     end
     st2=' : '; st1=[st1 st2(i1,:) char(table{:,3})  ...
         st2(i1,1) char(table{:,4})];
     i2=size(st1,2);
     for j1=1:size(table,1)
      st2='';
      if ~isempty(table{j1,5})
       st2=sprintf('  MatId %s ',table{j1,5});
      end
      if ~isempty(table{j1,6})
       st2=sprintf('%s  ProId %s ',st2,table{j1,6});
      end
      if ~isempty(st2) st1(j1,i2+[1:length(st2)])=st2;end
     end
     st1(find(st1==0))=' ';
     st1=cellstr(st1);st1=sprintf('\n%s',st1{:});
     st=[st st1];

   end
  end % model is empty

 if nargout==0 disp(st);
 elseif nargout==1 out=st; 
 elseif nargout==2 out=st; out1=table; end

% 'InfoNode' - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
elseif comstr(Cam,'node')

  [carg,FEnode,FEelt,ModelStack]=get_nodeelt(varargin,carg,ModelStack);

  opt = comstr(comstr(CAM,5),[-1 1]);
  [EGroup,nGroup]=getegroup(FEelt);
  [i1,i2]=find(any(FEelt==opt(1),2));
  i1=sort([i1;EGroup(find(diff(EGroup)==1))']);
  i3=find(FEnode(:,1)==opt(1));
  fprintf(1,'\n')
  if isempty(i3)
   fprintf(1,'node %i not found !!!\n\n',opt(1));
  else
   out=sprintf('node [%i,%i %i %i,%.12g %.12g %.12g]\n\n',FEnode(i3,:));
   fprintf(1,out);
  end
  jGroup=0;
  for j1=1:length(i1)
    if max(find(EGroup<=i1(j1)))~=jGroup
      jGroup = max(find(EGroup<=i1(j1)));
      ElemF= feutil('getelemf',FEelt(EGroup(jGroup),:),jGroup);
      i3=fe_super('node',ElemF); 
    end
    if isempty(find(FEelt(i1(j1),:)==opt(1)))
    elseif ~any(min(find(FEelt(i1(j1),:)==opt(1)))==i3)
    elseif EGroup(jGroup)==i1(j1) % single superelement
      [i2]=fe_super('node',ElemF);
      if ~isempty(i2) & any(i2==opt(1))
       fprintf(1,'e %i (g %i : %s) single superelement\n',i1(j1),jGroup,ElemF);
      end
    else
     fprintf(1,'e %i (g %i : %s), ',i1(j1),jGroup,ElemF);
     i2=FEelt(i1(j1),:);fprintf(1,' %i',i2(1:length(i3)));
     fprintf(1,' %g',i2(length(i3)+1:end));
     fprintf(1,'\n')
    end    
  end


else out='unknown'; end % 'info' subcommand selection - - - - - - - - - - - -

% ------------------------------------------------------------------------
elseif comstr(Cam,'join')   [CAM,Cam]=comstr(CAM,5);

elt=varargin{carg};
[EGroup,nGroup]=getegroup(elt);
if comstr(Cam,'group')
   opt=comstr(CAM(6:length(CAM)),-1);st=[];
   if ~isempty(opt) opt=round(opt(:));opt=opt(find(opt>0&opt<=nGroup));end
   if length(opt)==1 disp(['Cannot join a single group' 7]); end
   st = feutil('getelemf',elt(EGroup(opt(1)),:),opt(1));
   i4=EGroup(opt(1));
else st = setstr(CAM); opt=[0 0];i4=[];
  if isempty(st)
      st = feutil('getelemf',elt(EGroup(1),:),1);
  end
end
if length(opt)>1
  i3=[]; i2=[];
  for jGroup = 1:nGroup %loop on element groups
   if strcmp(feutil('getelemf',elt(EGroup(jGroup),:),jGroup),st)
    if opt(1)==0 | any(opt==jGroup)
     i2=[i2 EGroup(jGroup)+1:EGroup(jGroup+1)-1];
    else i3=[i3 EGroup(jGroup):EGroup(jGroup+1)-1]; end
    if isempty(i4) i4=EGroup(jGroup); end
   elseif any(opt==jGroup)
      disp(sprintf(['Cannot join group %i : different element type' 7],jGroup))
      i3=[i3 EGroup(jGroup):EGroup(jGroup+1)-1];
   else % not a group to join
     i3=[i3 EGroup(jGroup):EGroup(jGroup+1)-1];
   end
  end % of loop on jGroup
end

if isempty(i2) disp(sprintf('Found no groups of ''%s'' elements',st));
else
  i2=[i3(find(i3<min(i2))) i4 i2 i3(find(i3>min(i2)))];
  elt=elt(i2,:);
end
out=elt;


% ------------------------------------------------------------------------
elseif comstr(Cam,'matcomp')  

r1 = varargin{carg};carg=carg+1; r2 = varargin{carg};carg=carg+1;
i1 = find(r1); i2=find(r2); 
if all(size(i1)==size(i2)) & all(i1==i2)
else  st='different sparsity'; end

d=r1-r2; 
[i3,i4] = find(r1-r2); for j1=1:length(i3)
 i5(j1)= (d(i3(j1),i4(j1))^2/r1(i3(j1),i3(j1))*r1(i4(j1),i4(j1)));
end
if ~isempty(i3) st=sprintf('%g ',max(i5));
else st='matrices are identical'; end
if nargout==1 out=st; else fprintf('%s\n',st); clear out; end

% ------------------------------------------------------------------------
elseif comstr(Cam,'nnode')  

r1=varargin{carg};carg=carg+1;
if isfield(r1,'Node'); r1=r1.Node;end
out=sparse(r1(:,1),1,1:size(r1,1));

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

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
% ObjectBeamLine
if comstr(Cam,'beamline') [CAM,Cam] = comstr(CAM,9);

if strfind(Cam,'-egid') % Get The EGID if possible
  i1=strfind(Cam,'-egid'); 
  [EGID,i2,i5,i3]=sscanf(CAM(i1+5:end),'%i',1);
  EGID=[0 EGID];
  CAM(i1:i1+i3+3)='';[CAM,Cam] = comstr(CAM,1);
else; EGID=[]; end
if isempty(Cam)|Cam(1)~='-' st='beam1';
else [st,i1,i2,i3]=sscanf(CAM(2:end),'%s',1);CAM=CAM(i3+1:end);
end

opt = comstr(CAM,[-1]); 

if isempty(opt)&carg<=nargin opt=varargin{carg};carg=carg+1; end

opt=opt(:);opt(end+1)=0;
if length(opt)==0 error('You need at least two nodes for a beam line');end
if carg<=nargin i2=varargin{carg};carg=carg+1;else i2=[1 1];end
if ~isnumeric(i2)|size(i2,1)~=1 i2=[1 1];end

elt=[Inf abs(st) EGID];
elt(2:length(opt),1:2+length(i2))=[opt(1:length(opt)-1) opt(2:length(opt)) ...
       ones(length(opt)-1,1)*i2];
i1=find(elt(2:end-1,1)==0&elt(3:end,2)==0 ...
        &elt(2:end-1,2)==elt(3:end,1));elt(i1+1,1)=elt(i1+1,2);
i1=find(elt(:,1)&elt(:,2)); elt=elt(i1,:);

out=elt;

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
% model=feutil('ObjectCircle x y z r nx ny nz Nseg',model)
elseif comstr(Cam,'circle');[CAM,Cam]=comstr(CAM,7);
 
 r1=comstr(Cam,[-1]);
 if carg<=nargin&isstruct(varargin{carg}); model=varargin{carg};carg=carg+1;
 else; model=struct('Node',[],'Elt',[]);
 end

 if length(r1)~=8; error(' x y z r nx ny nz Nseg must be provided');end

 r2=linspace(0,2*pi,r1(8)+1);  p=basis(r1(5:7),[0 0 0],0);p=p(:,[2 3 1]);
 r2=r1(4)*p*[cos(r2);sin(r2);zeros(size(r2))];
 r2=r2'+ones(size(r2,2),1)*r1(1:3);

 [model.Node,i2]=feutil('addnode',model.Node,r2);
 i2=model.Node(i2,1);
 i1=[i2(1:end-1) i2(2:end) ones(size(i2,1)-1,1)*[1 1 0 0]];
 model.Elt(end+[1:size(i1,1)+1],1:6)=[Inf abs('beam1');i1];
 out=model;

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
% model=feutil('ObjectDisk x y z r nx ny nz NsegT NsegR',model)
elseif comstr(Cam,'disk');[CAM,Cam]=comstr(CAM,5);
 
 r1=comstr(Cam,[-1]);
 if carg<=nargin&isstruct(varargin{carg}); model=varargin{carg};carg=carg+1;
 else; model=struct('Node',[],'Elt',[]);
 end
 if length(r1)~=9; error(' x y z r nx ny nz NsegT NsegR must be provided');end
 r3=linspace(0,r1(4),r1(9)+1)';
 r2=linspace(0,2*pi,r1(8)+1);  p=basis(r1(5:7),[0 0 0],0);p=p(:,[2 3 1]);
 r2=[cos(r2);sin(r2);zeros(size(r2))];
 i1=length(r3)*length(r2);
 r2=[reshape(r3*p(1,:)*r2,i1,1)+r1(1) reshape(r3*p(2,:)*r2,i1,1)+r1(2) ...
     reshape(r3*p(3,:)*r2,i1,1)+r1(3)];

 [model.Node,i2]=feutil('addnode',model.Node,r2);
 i2=model.Node(i2,1);
 i1=length(r3);i1=[1:i1-1;2:i1;i1+2:2*i1;i1+1:2*i1-1];i1=i1(:);
 i3=[0:r1(8)-1]*length(r3);
 i1=i1(:,ones(r1(8),1))+i3(ones(size(i1,1),1),:);
 i1=reshape(i2(i1),4,prod(size(i1))/4)';i1(:,5:6)=1;
 model.Elt(end+[1:size(i1,1)+1],1:6)=[Inf abs('quad4');i1];
 out=model;

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
% model=feutil('ObjectAnnulus x y z r1 r2 nx ny nz NsegT NsegR',model)
elseif comstr(Cam,'annulus');[CAM,Cam]=comstr(CAM,8);
 
 r1=comstr(Cam,[-1]);
 if carg<=nargin&isstruct(varargin{carg}); model=varargin{carg};carg=carg+1;
 else; model=struct('Node',[],'Elt',[]);
 end
 if length(r1)~=10; error(' x y z r1 r2 nx ny nz NsegT NsegR must be provided');end
 r3=linspace(r1(4),r1(5),r1(10)+1)';
 r2=linspace(0,2*pi,r1(9)+1);  p=basis(r1(6:8),[0 0 0],0);p=p(:,[2 3 1]);
 r2=[cos(r2);sin(r2);zeros(size(r2))];
 i1=length(r3)*length(r2);
 r2=[reshape(r3*p(1,:)*r2,i1,1)+r1(1) reshape(r3*p(2,:)*r2,i1,1)+r1(2) ...
     reshape(r3*p(3,:)*r2,i1,1)+r1(3)];

 [model.Node,i2]=feutil('addnode',model.Node,r2);
 i2=model.Node(i2,1);
 i1=length(r3);i1=[1:i1-1;2:i1;i1+2:2*i1;i1+1:2*i1-1];i1=i1(:);
 i3=[0:r1(9)-1]*length(r3);
 i1=i1(:,ones(r1(9),1))+i3(ones(size(i1,1),1),:);
 i1=reshape(i2(i1),4,prod(size(i1))/4)';i1(:,5:6)=1;
 model.Elt(end+[1:size(i1,1)+1],1:6)=[Inf abs('quad4');i1];
 out=model;

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
% model=feutil('ObjectCylinder x1 y1 z1 x2 y2 z2 R NsegR NsegT',model)
elseif comstr(Cam,'cylinder');[CAM,Cam]=comstr(CAM,9);
 
 r1=comstr(Cam,[-1]);
 if carg<=nargin&isstruct(varargin{carg}); model=varargin{carg};carg=carg+1;
 else; model=struct('Node',[],'Elt',[]);
 end

 if length(r1)~=9; 
  error('x1 y1 z1 x2 y2 z2 R NsegR NsegT must be provided');
 end

 r2=linspace(0,2*pi,r1(8)+1);
 p=basis(r1(4:6)-r1(1:3),[0 0 0],0);p=p(:,[2 3 1]);
 r2=r1(7)*p*[cos(r2);sin(r2);zeros(size(r2))];
 r2=r2'+ones(size(r2,2),1)*r1(1:3);

 % build nodes and topology
 r3=zeros(size(r2,1)*r1(9),3);r4=ones(size(r2,1),1)*(r1(4:6)-r1(1:3))/r1(9);
 i2=size(r2,1);ind=[1:i2]';r3(ind,:)=r2;i1=ones((i2-1)*(r1(9)),6);
 for j1=1:r1(9); 
  i1((1:i2-1)+(i2-1)*(j1-1),[1 4 2 3])=[ ...
    ind(1:end-1)+i2*(j1-1) ind(1:end-1)+i2*j1 ...
    ind(2:end)+i2*(j1-1) ind(2:end)+i2*j1];
  r3(ind+j1*size(r2,1),:)=r2+r4*j1; 
 end

 [model.Node,i2]=feutil('addnode',model.Node,r3);
 i2=model.Node(i2,1); i1(:,1:4)=reshape(i2(i1(:,1:4)),size(i1,1),4);
 model.Elt(end+[1:size(i1,1)+1],1:6)=[Inf abs('quad4');i1];
 out=model;


% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
% ObjectWireFrame
elseif comstr(Cam,'wireframe');
 
elt=[];  i2=[1 1];
while carg<=nargin

 opt=varargin{carg};carg=carg+1;
 opt=opt(:);opt(end+1)=0;
 if isempty(elt); elt=[Inf abs('beam1') 0 -1];
 else; elt(end+1,:)=elt(1,:);
 end

 opt=[opt(1:length(opt)-1) opt(2:length(opt)) ...
       ones(length(opt)-1,1)*i2];
 i1=find(opt(1:end-1,1)==0&opt(2:end,2)==0 ...
        &opt(1:end-1,2)==opt(2:end,1));opt(i1,1)=opt(i1,2);
 opt=opt(find(opt(:,1)&opt(:,2)),:);
 elt(end+[1:size(opt,1)],1:size(opt,2))=opt;

end
out=elt;

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
% ObjectMass
elseif comstr(Cam,'mass'); 

if length(Cam)>4 & ~isempty(strfind(Cam,'egid'))
 i2=strfind(Cam,'egid');
 [i1,i3,i4,i5]=sscanf(CAM(i2+4:end),'%d',1);
 [CAM,Cam]=comstr(CAM([1:i2-1 i2+3+i5:end]),1);
else i1=0;
end

opt = comstr(comstr(CAM,5),[-1]);
if isempty(opt)&carg<=nargin opt=varargin{carg};carg=carg+1; end
opt=opt(:);opt=opt(find(opt));

out=[Inf abs('mass1') 0]; if i1 out=[out i1];end
out(1+[1:length(opt)],1) = opt(:);

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
% [model,node,carg]=feutil('objectget',carg,varargin{:});
elseif comstr(Cam,'get'); [CAM,Cam]=comstr(CAM,4);

 if isstruct(varargin{carg+2}); model=varargin{carg+2};carg=carg+1;
 else; model=struct('Node',[],'Elt',[]);
 end

 node=varargin{carg+2};carg=carg+1; 
 if size(node,2)==7; node=node(:,5:7);
 elseif size(node,2)==1 % nodeID given
  node=feutil('getnode',model,node); node=node(:,5:7);
 elseif size(node,2)==3
 else; error('Bad format for nodes');
 end
 out=model; out1=node; out2=carg;

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
% model=feutil('objectbeam',model,[Corners],dvx,dvy);
elseif comstr(Cam,'beam'); [CAM,Cam]=comstr(CAM,5);

 i3=comstr(Cam,[-1 1 1]);
 [model,node,carg]=feutil('objectget',carg,varargin{:});

 dv1=varargin{carg};carg=carg+1; 
 if length(dv1)==1 opt(1)=dv1+1;dv1=linspace(0,1,opt(1));
 else; opt(1)=length(dv1);
 end
 [r1,i1]=feutil('objectrefcube',dv1);r1(:,4)=0;
 rule=integrules('bar1',r1); 
 [model.Node,i2]=feutil('addnode',model.Node,rule.N*node);
 i2=model.Node(i2,1); i1=reshape(i2(i1),size(i1,1),size(i1,2));
 i1(:,3)=i3(1);i1(:,4)=i3(2);i1(1,6)=0;

 model.Elt(end+[1:size(i1,1)+1],1:6)=[Inf abs('beam1');i1];
 out=model;

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
% model=feutil('objectquad',model,[Corners],dvx,dvy);
elseif comstr(Cam,'quad'); [CAM,Cam]=comstr(CAM,5);

 i3=comstr(Cam,[-1 1 1]);
 [model,node,carg]=feutil('objectget',carg,varargin{:});

 dv1=varargin{carg};carg=carg+1; 
 if length(dv1)==1; opt(1)=dv1+1;dv1=linspace(0,1,opt(1));
 else; opt(1)=length(dv1);
 end
 dv2=varargin{carg};carg=carg+1; 
 if length(dv2)==1 opt(2)=dv2+1;dv2=linspace(0,1,opt(2));
 else; opt(2)=length(dv2);
 end

 [r1,i1]=feutil('objectrefcube',dv1,dv2); % Reference nodes and connectivity
 r1(:,3:4)=.5;rule=integrules('quad4',r1*2-1); 
 [model.Node,i2]=feutil('addnode',model.Node,rule.N*node);
 i2=model.Node(i2,1); i1=reshape(i2(i1),size(i1,1),size(i1,2));
 i1(:,5)=i3(1); i1(:,6)=i3(2);
 model.Elt(end+[1:size(i1,1)+1],1:6)=[Inf abs('quad4');i1];
 out=model;

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
% ObjectHexa
% model=feutil('objecthexa',[Oxyz;0Axyz;OBxyz;OCxyz],dvx,dvy,dvz);
elseif comstr(Cam,'hexa'); [CAM,Cam]=comstr(CAM,5);

 i3=comstr(Cam,[-1 1 1]);
 [model,node,carg]=feutil('objectget',carg,varargin{:});

 if size(node,1)==4 %[Oxyz 0Axyz OBxyz OCxyz]
  node=node(ones(8,1),:)+ ...
  [0 0 0;node(2,:);node(2,:);0 0 0;0 0 0;node(2,:);node(2,:);0 0 0]+ ...
  [0 0 0;0 0 0;node([3 3],:);zeros(2,3);node([3 3],:)]+...
  [zeros(4,3);node([4 4 4 4],:)];
 end

 dv1=varargin{carg};carg=carg+1; 
 if length(dv1)==1 opt(1)=dv1+1;dv1=linspace(0,1,opt(1));
 else; opt(1)=length(dv1);
 end
 dv2=varargin{carg};carg=carg+1; 
 if length(dv2)==1 opt(2)=dv2+1;dv2=linspace(0,1,opt(2));
 else; opt(2)=length(dv2);
 end
 dv3=varargin{carg};carg=carg+1; 
 if length(dv3)==1 opt(3)=dv3+1;dv3=linspace(0,1,opt(3));
 else; opt(3)=length(dv3);
 end

 [r1,i1]=feutil('objectrefcube',dv1,dv2,dv3); % Reference nodes and connectivity
 r1(:,4)=.5; rule=integrules('hexa8 node',r1*2-1); 
 [model.Node,i2]=feutil('addnode',model.Node,rule.N*node);
 i2=model.Node(i2,1); i1=reshape(i2(i1),size(i1,1),size(i1,2));
 i1(:,9)=i3(1); i1(:,10)=i3(2);
 model.Elt(end+[1:size(i1,1)+1],1:10)=[Inf abs('hexa8') 0 0 0 0;i1];
 out=model;


% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
% ObjectRefCube : generates nodes and connectivity of a reference cube
elseif comstr(Cam,'refcube'); 

 dv1=varargin{carg};carg=carg+1; 
 if length(dv1)==1 opt(1)=dv1;dv1=linspace(0,1,opt(1));
 else; opt(1)=length(dv1);
 end
 if carg>nargin dv2=[]; else; dv2=varargin{carg};carg=carg+1; end
 if length(dv2)==1 opt(2)=dv2;dv2=linspace(0,1,opt(2));
 else; opt(2)=length(dv2);
 end
 if carg>nargin dv3=[]; else; dv3=varargin{carg};carg=carg+1; end
 if length(dv3)==1 opt(3)=dv3;dv3=linspace(0,1,opt(3));
 else; opt(3)=length(dv3);
 end
 if length(unique(dv1))~=length(dv1);warning('degenerate r division');end
 if length(unique(dv2))~=length(dv2);warning('degenerate s division');end
 if length(unique(dv3))~=length(dv3);warning('degenerate t division');end

if isempty(dv2) % 1D generation of a beam

  out=[dv1(:)*[1 0 0]];
  out1=[1:size(out,1)-1;2:size(out,1)]';

elseif isempty(dv3) % 2D generation of a quad4  - - - - - - - - - - -

     [r1,r2]=ndgrid(dv1,dv2);
     out=[r1(:) r2(:)*[1 0]];
     i4 = [1:opt(1)-1];i4=i4(ones(opt(2)-1,1),:)';
     i5=[0:opt(2)-2]*(opt(1));i5=i5(ones(opt(1)-1,1),:);
     i4=i4(:)+i5(:); i4 =[i4 i4+opt(1) i4+opt(1)+1 i4+1];
     out1=i4;

else  % 3D generation of an hexa8  - - - - - - - - - -
 [r1,r2,r3]=ndgrid(dv1,dv2,dv3);
 out=[r1(:) r2(:) r3(:)];

     % row of elements
     i4 = [1:opt(1)-1;2:opt(1)]';
     i4=[i4 i4(:,[2 1])+opt(1)];i4=[i4 i4+opt(1)*opt(2)];
     % face of elements
     i5=[0:opt(1):opt(1)*(opt(2)-1)-1];i5=i5(ones(size(i4,1),1),:);i5=i5(:);
     i6=[1:size(i4,1)]';i6=i6(:,ones(opt(2)-1,1)); i6=i6(:);
     i4=i4(i6,:)+i5(:,ones(size(i4,2),1));
     % volume of elements
     i5=[0:opt(1)*opt(2):opt(1)*opt(2)*(opt(3)-1)-1];
     i5=i5(ones(size(i4,1),1),:);i5=i5(:);
     i6=[1:size(i4,1)]';i6=i6(:,ones(opt(3)-1,1)); i6=i6(:);
     i4=i4(i6,:)+i5(:,ones(size(i4,2),1));
 out1=i4;

end

else out='unknown';end % subcommand selection - - - - - - - - - - - - - - -

% ------------------------------------------------------------------------
% optimisation functions
elseif comstr(Cam,'optim') [CAM,Cam] = comstr(CAM,6);

[carg,FEnode,FEelt,FEel0,ModelStack]=get_nodeelt(varargin,carg,ModelStack);

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
% OptimModel: removes unused nodes from FEnode
if comstr(Cam,'model')

  [EGroup,nGroup]=getegroup(FEelt);
  %NNode(FEnode(:,1))=[1:length(FEnode(:,1))]';
  NNode=sparse(FEnode(:,1),1,1:length(FEnode));
  disp(sprintf('Model %s with %i element groups','FEelt',nGroup));

  i2 =[]; % nodes used
  for jGroup = 1:nGroup %loop on element groups
   ElemF= feutil('getelemf',FEelt(EGroup(jGroup),:),jGroup);
   i3=fe_super('node',ElemF);
   i3 = FEelt(EGroup(jGroup)+1:EGroup(jGroup+1)-1,i3);
   if ~all(all(i3)) i3=i3(find(i3));end
   i2 = [i2;i3(:)];
  end % of jGroup loop
  % indices in FEnode of the retained nodes
  i2 = NNode(full(find(sparse(i2,ones(size(i2)),ones(size(i2))))));
  i2 = i2(:); if any(~finite(FEnode(end,:))) i2=[i2;size(FEnode,1)];end
  FEnode = FEnode(i2,:);
  out=FEnode; out1=FEelt;

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
% OptimNodeNum: optimizes model node numbering
elseif comstr(Cam,'nodenum')

  [EGroup,nGroup]=getegroup(FEelt);
  NNode(FEnode(:,1))=[1:length(FEnode(:,1))]';
  i2 =[]; i4=size(FEnode,1); 
  for jGroup = 1:nGroup %loop on element groups
   ElemF= feutil('getelemf',FEelt(EGroup(jGroup),:),jGroup);
   i3=fe_super('node',ElemF);
   i3 = FEelt(EGroup(jGroup)+1:EGroup(jGroup+1)-1,i3);
   i5 = find(i3); i3(i5)=NNode(i3(i5));
   i5=[];for j1=1:size(i3,2)
    i5=[i5 (i3(:,ones(size(i3,2)-j1+1,1)*j1)-1)*i4+i3(:,j1:end)];
   end; i2=[i2;i5(:)];
  end % of jGroup loop

  k=reshape(sparse(i2,1,1,i4^2,1),i4,i4); k = k+k';
  i3=symrcm(k); 
  %subplot(211);spy(k); subplot(212);spy(k(i3,i3))
  FEnode=FEnode(i3,:);
  out=FEnode; out1=FEelt;

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
% OptimEltCheck: checks the conformity of elements (twisted things, etc).
% check/fix ill-conditionned hexa8 / penta6 / quad4
% Contributed by Mathieu Corus
elseif comstr(Cam,'eltcheck')
  
  if nargin==1
    test_fe_check;
    return;
  elseif nargin==2
    model=varargin{2}; if ~isstruct(model) error('Please input a structure with .Node and .Elt fields'); end;
    bypass=0;
  elseif nargin==3
    if isnumeric(varargin{2}) & isnumeric(varargin{3})
      model.Node=varargin{2};
      model.Elt=varargin{3};
      bypass=0;
    else
      model=varargin{2}; if ~isstruct(model) error('Please input a structure with .Node and .Elt fields'); end;
      bypass=varargin{3};
    end;
  elseif nargin==4
    if isnumeric(varargin{2}) & isnumeric(varargin{3}) & isnumeric(varargin{4})
      model.Node=varargin{2};
      model.Elt=varargin{3};
      bypass=varargin{4};
    else
      error('Bad argument types');
    end;
  else
    error('Too many input arguments');
  end;

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

%-- loop on hexa8 / penta6 / quad4 groups

Gr_inf=find(model.Elt(:,1)==Inf);
Gr_inf=[Gr_inf;length(model.Elt(:,1))+1];
ind_new_deg_quad=[];
new_deg_quad=[];
for i1=1:length(Gr_inf)-1
  if comstr(comstr(char(model.Elt(Gr_inf(i1),2:end)),-27),'hexa8')
    elt_type='hexa8';ind_elt=1:8;
    for j1=(Gr_inf(i1)+1):(Gr_inf(i1+1)-1)
      model.Elt(j1,ind_elt)=elt_check(model.Elt(j1,ind_elt),elt_type,model,NNode,bypass);
    end;
  elseif comstr(comstr(char(model.Elt(Gr_inf(i1),2:end)),-27),'penta6')
    elt_type='penta6';ind_elt=1:6;
    for j1=(Gr_inf(i1)+1):(Gr_inf(i1+1)-1)
      model.Elt(j1,ind_elt)=elt_check(model.Elt(j1,ind_elt),elt_type,model,NNode,bypass);
    end;
  elseif comstr(comstr(char(model.Elt(Gr_inf(i1),2:end)),-27),'quad4')
    elt_type='quad4';ind_elt=1:4;
    for j1=(Gr_inf(i1)+1):(Gr_inf(i1+1)-1)
      nodes=elt_check(model.Elt(j1,ind_elt),elt_type,model,NNode,bypass);
      if size(nodes,1)==2 
        ind_new_deg_quad=[ind_new_deg_quad;j1];
        new_deg_quad=[new_deg_quad ; nodes];
      else
        model.Elt(j1,ind_elt)=nodes;
      end;
    end;
  end;
end;

%-- add new lines in model.Elt if necessary --%

if ~isempty(ind_new_deg_quad)
  ind_new_deg_quad=ind_new_deg_quad+(0:(length(ind_new_deg_quad)-1))';
  for i1=1:length(ind_new_deg_quad);
    model.Elt=model.Elt([1:ind_new_deg_quad(i1) ind_new_deg_quad(i1):end],:);
    model.Elt(ind_new_deg_quad(i1)+(0:1),1:4)=new_deg_quad(2*i1-1:2*i1,:);
  end;
end;
out=model;


else out='unknown'; end % subcommand selection - - - - - - - - - - - - - -


% ------------------------------------------------------------------------
% 'Orient' - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
elseif comstr(Cam,'orient'); [CAM,Cam]=comstr(CAM,7);

  [carg,node,elt,el0,ModelStack]=get_nodeelt(varargin,carg,ModelStack);
  [EGroup,nGroup]=getegroup(elt);
  NNode=[];NNode(node(:,1)+1)=1:size(node,1);

  % give element groups and/or nodes
  RunOpt=struct('neg',1,'Group',[]);
  i1=strfind(Cam,'-neg'); if ~isempty(i1);
    RunOpt=struct('neg',-1); CAM(i1+[0:3])='';[CAM,Cam]=comstr(CAM,1); 
  end
  i1=strfind(Cam,'nlist'); if ~isempty(i1);
    RunOpt.GiveList=2; CAM(i1+[0:4])='';[CAM,Cam]=comstr(CAM,1); 
    RunOpt.List=[];
  else; RunOpt.GiveList=0;
  end
  i1=strfind(Cam,'gauss'); if ~isempty(i1);
    RunOpt.GiveList=1; CAM(i1+[0:4])='';[CAM,Cam]=comstr(CAM,1); 
  end

  i1=strfind(Cam,'n'); if isempty(i1) i1=length(CAM)+1; end
  RunOpt.Group=comstr(Cam(1:i1-1),[-1]); i4=0;
  RunOpt.Warped=[];

  if isempty(RunOpt.Group) % orientation based on k-matrix sign - - - - - - 

  i5=zeros(size(elt,1),1);

  for jGroup = 1:nGroup %loop on groups
     [ElemF,i1,ElemP]=feutil('getelemf',elt(EGroup(jGroup),:),jGroup);
     i1=fe_super('node',ElemP);
     cEGI = [EGroup(jGroup)+1:EGroup(jGroup+1)-1]';
     %i2(DOF) i3(node)
     switch ElemP
      case {'tria3','t3p'}; ElemP=ElemF; [in2,in1]=tria3('flip');  i2=6; i3=3;
      case {'tria6','t6p'}; ElemP=ElemF; [in2,in1]=tria6('flip');  i2=12;i3=6;
      case {'quad4','q4p'}; ElemP=ElemF; [in2,in1]=quad4('flip');  i2=8; i3=4;
      case {'quadb','q8p'}; ElemP=ElemF; [in2,in1]=quadb('flip');  i2=16; i3=8;
      case {'quad9','q9a'}; ElemP=ElemF; [in2,in1]=quad9('flip');  i2=18;i3=9;
      case 'q5p';                        [in2,in1]=q5p('flip');    i2=10;i3=5;
      case {'tetra4','flui4'} ;          [in2,in1]=tetra4('flip'); i2=12;i3=4;
      case 'tetra10';                    [in2,in1]=tetra10('flip');i2=30;i3=10;
      case {'penta6','flui6'};           [in2,in1]=penta6('flip'); i2=18;i3=6;
      case 'penta15' ;                   [in2,in1]=penta15('flip');i2=45;i3=15;
      case {'hexa8','flui8'} ;           [in2,in1]=hexa8('flip');  i2=24;i3=8;
      case {'hexa20'} ;               [in2,in1]=hexa20('flip');  i2=60;i3=20;
      case {'beam1','bar1','mass1','mass2'}
      otherwise 
       sdtw('Orientation test not implemented for %s',ElemP);
     end
     % this should be rewritten using calls of BuildNDN 
     if any(strcmp(ElemP,...
             {'penta6','penta15','tetra4','tetra10','hexa8','hexa20'}))
      if RunOpt.GiveList==1; EltConst=integrules(ElemP); 
      else; EltConst=integrules(ElemP,'node'); 
      end; % if change keep same for GiveList
      for jElt=1:length(cEGI) %loop on elements of group
        nodeE = node(NNode(elt(cEGI(jElt),i1)+1),[5:7 4]);%xxx 4? also penta?
        of_mk('buildndn',3,EltConst,nodeE);
        %k=of_mk(ElemP,int32([i2*(i2+1)/2*[1 1] 0 0 0 0 0 0 0]), ...
        %  int32([0 0 i2 i3 1 0 1 2]),[7800 0 2e11 .3 0],nodeE);
        if all(EltConst.jdet>=0); %k(1)<0  
        elseif RunOpt.GiveList==2
         RunOpt.List(end+[1 2],[1:length(i1)+1])= ...
          [cEGI(jElt) elt(cEGI(jElt),i1);0 EltConst.jdet(:)'];
        elseif all(EltConst.jdet<=0)
         elt(cEGI(jElt),in1)=elt(cEGI(jElt),in2);i5(cEGI(jElt))=1; 
        else; 
         RunOpt.Warped(end+1,1:size(elt,2))=elt(cEGI(jElt),:);
        end
      end % loop on elements of group
     elseif any(strcmp(ElemP,{'t3p','t6p','q4p','q5p'}))
      for jElt=1:length(cEGI) %loop on elements of group
        r1 = node(NNode(elt(cEGI(jElt),i1)+1),[5:7 4]);
        if ~any(r1(:,2))&any(r1(:,3)) r1=r1(:,[1 3 2 4]); end
        k=of_mk(ElemP,int32([i2*(i2+1)/2 0 0 0 0 0 0 0 0]), ...
           int32([0 0 i2 i3 1 0 1 2]),[7800 2e11 .3],r1);
        if k(1)<0;
          elt(cEGI(jElt),in1)=elt(cEGI(jElt),in2);i5(cEGI(jElt))=1;
        end
      end % loop on elements of group
     end % elements to run
  end % loop on groups
  if ~isempty(RunOpt.Warped);
   disp(RunOpt.Warped);error('There were some warped volumes, list above');
  end
  if RunOpt.GiveList==2 out=RunOpt.List;
  else; 
   out=elt; out1=find(i5); i5=nnz(i5);
   if i5&~any(Cam==';') fprintf('%i elements reoriented\n',i5);end
  end
  % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  else % surface orientation based on nodes given ...

  out1=[];
  onode = comstr(Cam(i1+1:end),[-1]);
  if isempty(onode)&carg<=nargin & isa(varargin{carg},'double')
    onode=varargin{carg};carg=carg+1;
    if isempty(onode)&carg<=nargin
      onode=varargin{carg};carg=carg+1;
    end
    if size(onode,2)==7 onode=onode(:,5:7);end
  elseif rem(prod(size(onode)),3)==0 
    onode=reshape(onode,3,prod(size(onode))/3)';
  end
  fprintf('Orientation nodes \n');disp(onode)
  [EGroup,nGroup]=getegroup(elt);
  NNode(node(:,1))=[1:length(node(:,1))]';
  for jGroup = RunOpt.Group(find(RunOpt.Group<=nGroup)) %loop on groups
   [ElemF,i1,ElemP]=feutil('getelemf',elt(EGroup(jGroup),:),jGroup);
   i1=fe_super('node',ElemP);
   cEGI = [EGroup(jGroup)+1:EGroup(jGroup+1)-1]';
   if comstr(ElemP,'quad4')|comstr(ElemP,'quadb')| ...
      comstr(ElemP,'tria3')
      for jElt=1:length(cEGI) %loop on elements of group
        i2=elt(cEGI(jElt),i1);i2=NNode(i2(find(i2)));
        r1 = node(i2,5:7); r2 = mean(r1);r1=r1-r2(ones(size(r1,1),1),:);
        [bas,r1] = basis(r1);
        out1(end+1,1:6)=[r2 bas(:,3)'];
        r1 = onode-r2(ones(size(onode,1),1),:);
        [r2,i3]=sort(sum(r1.^2,2));

        % change orientation of surfaces based on normal
        if r1(i3(1),:)*bas(:,3)*RunOpt.neg<0  
	  if     comstr(ElemP,'quadb')
            elt(cEGI(jElt),[1:8])=elt(cEGI(jElt),[3 2 1 4 6 5 8 7]);
          elseif comstr(ElemP,'quad4')
            if ~diff(elt(cEGI(jElt),[3:4])) % degenerate
             elt(cEGI(jElt),[1 3 4])=elt(cEGI(jElt),[3 1 1]);
            else elt(cEGI(jElt),[1 3])=elt(cEGI(jElt),[3 1]);end
          else % tria3
            elt(cEGI(jElt),[1 3])=elt(cEGI(jElt),[3 1]);
          end
          i4 = i4+1; out1(end,4:6)=[-bas(:,3)'];
        end
      end
   elseif comstr(ElemP,'beam3')
      for jElt=1:length(cEGI) %loop on elements of group
        i2=elt(cEGI(jElt),i1);i2=NNode(i2(find(i2)));
        r1 = node(i2(1:2),5:7);
        r2 = mean(r1);r1=r1-r2(ones(size(r1,1),1),:);
        bas = basis(r1(2,:),[0 0 0],1);
        out1(end+1,1:6)=[r2 bas(:,1)'];
        r1 = onode-r2(ones(size(onode,1),1),:);
        [r2,i3]=sort(sum(r1'.^2));
        % change orientation of lines based on axis
        if r1(i3(1),:)*bas(:,1)*RunOpt.neg<0
            elt(cEGI(jElt),1:2)=elt(cEGI(jElt),[2 1]);
            out1(end,4:6)=bas(:,1)';  i4 = i4+1;
        end
      end
   elseif any(strcmp(ElemP,{'beam1','bar1'}))
   else
       error(sprintf('I don''t know how to orient %s',ElemP));
   end % supported
  end % loop on groups
  disp(sprintf('%i  element(s) reoriented',i4));
  out=elt;

  end %of orientation type

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

  [carg,FEnode,FEelt,FEel0,ModelStack]=get_nodeelt(varargin,carg,ModelStack);
  [EGroup,nGroup]=getegroup(FEelt);

  r1=[.5 .5 0 0;0 .5 .5 0;0 0 .5 .5;.5 0 0 .5];
  i2 =[]; % i2 first group nodes, i3 second group nodes
  for jGroup = 1:nGroup %loop on element groups
   ElemF= feutil('getelemf',FEelt(EGroup(jGroup),:),jGroup);
   if strcmp(ElemF,'quad4')
     FEelt(EGroup(jGroup),1:7)=[Inf abs('quadb') 0];
     for j1=EGroup(jGroup)+1:EGroup(jGroup+1)-1
       NNode(FEnode(:,1))=[1:length(FEnode(:,1))]';
       r2 = r1*FEnode(NNode(FEelt(j1,1:4)),5:7);
       [FEnode,i2]=feutil('AddNode',FEnode,r2);
       r2=[FEelt(j1,1:4) FEnode(i2(:),1)' FEelt(j1,5:max(find(FEelt(j1,:))))];
       FEelt(j1,1:length(r2))=r2;
     end
   end
  end % of jGroup loop
  out=FEnode;out1=FEelt;

% ------------------------------------------------------------------------
% q42q9
elseif comstr(Cam,'quad42quad9')  [CAM,Cam] = comstr(CAM,7);

  [carg,FEnode,FEelt,FEel0,ModelStack]=get_nodeelt(varargin,carg,ModelStack);
  [EGroup,nGroup]=getegroup(FEelt);

  r1=[.5 .5 0 0;0 .5 .5 0;0 0 .5 .5;.5 0 0 .5;.25 .25 .25 .25];
  i2 =[]; % i2 first group nodes, i3 second group nodes
  for jGroup = 1:nGroup %loop on element groups
   ElemF= feutil('getelemf',FEelt(EGroup(jGroup),:),jGroup);
   if strcmp(ElemF,'quad4')
     FEelt(EGroup(jGroup),1:7)=[Inf abs('quad9') 0];
     for j1=EGroup(jGroup)+1:EGroup(jGroup+1)-1
       NNode(FEnode(:,1))=[1:length(FEnode(:,1))]';
       r2 = r1*FEnode(NNode(FEelt(j1,1:4)),5:7);
       [FEnode,i2]=feutil('AddNode',FEnode,r2);
       r2=[FEelt(j1,1:4) FEnode(i2(:),1)' FEelt(j1,5:max(find(FEelt(j1,:))))];
       FEelt(j1,1:length(r2))=r2;
     end
   end
  end % of jGroup loop
  out=FEnode;out1=FEelt;

% ------------------------------------------------------------------------
% q42q5
elseif comstr(Cam,'quad42q5p')  [CAM,Cam] = comstr(CAM,7);

  [carg,FEnode,FEelt,FEel0,ModelStack]=get_nodeelt(varargin,carg,ModelStack);
  [EGroup,nGroup]=getegroup(FEelt);

  r1=[.25 .25 .25 .25];
  i2 =[]; % i2 first group nodes, i3 second group nodes
  for jGroup = 1:nGroup %loop on element groups
   ElemF= feutil('getelemf',FEelt(EGroup(jGroup),:),jGroup);
   if strcmp(ElemF,'quad4')
     FEelt(EGroup(jGroup),1:7)=[Inf abs('q5p') 0 0 0];
     for j1=EGroup(jGroup)+1:EGroup(jGroup+1)-1
       NNode(FEnode(:,1))=[1:length(FEnode(:,1))]';
       r2 = r1*FEnode(NNode(FEelt(j1,1:4)),5:7);
       [FEnode,i2]=feutil('AddNode',FEnode,r2);
       r2=[FEelt(j1,1:4) FEnode(i2(:),1)' FEelt(j1,5:max(find(FEelt(j1,:))))];
       FEelt(j1,1:length(r2))=r2;
     end
   end
  end % of jGroup loop
  out=FEnode;out1=FEelt;

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

 if comstr(Cam,'elt') % 'RemoveElt' - - - - - - - - - - - - - -

  [carg,FEnode,FEelt,FEel0,ModelStack]=get_nodeelt(varargin,carg,ModelStack);
  [i4,out1]=feutil(['Find elt' CAM(4:end)],FEnode,FEelt,FEel0, ...
     varargin{carg:end});

  out=FEelt(setdiff(1:size(FEelt,1),i4),:);
  % remove empty groups
  i1=find(~finite(out(:,1)));i1=intersect(i1,i1-1);
  if ~finite(out(end,1));i1(end+1,1)=size(out,1);end
  if ~isempty(i1);out(i1,:)=[];end

 else out='unknown'; end % subcommand selection - - - - - - - - - - - - - -

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

model=varargin{carg};carg=carg+1;
if   carg>nargin % renumber nodes from 1 to N
  i1=1:size(model.Node,1); NNode=sparse(model.Node(:,1)+1,1,i1);
elseif carg+1<=nargin% .Node already numbered only fix NNode
   i1=varargin{carg};carg=carg+1;
   NNode=varargin{carg};carg=carg+1;
else % provide new node numbers in i1
  i1=varargin{carg};carg=carg+1;NNode=sparse(model.Node(:,1)+1,1,i1); 
end


% deal with local bases
% renumbering of reference to nodes would be good
if isfield(model,'bas')  
        [node,model.bas]=basis(model.Node,model.bas);
end

% renumber nodes in elements
[EGroup,nGroup]=getegroup(model.Elt);

for jGroup = 1:nGroup %loop on element groups
   ElemF= feutil('getelemf',model.Elt(EGroup(jGroup),:),jGroup);
   i3=fe_super('node',ElemF); ind=EGroup(jGroup)+1:EGroup(jGroup+1)-1;
   i2 = model.Elt(ind,i3); % nodes
   model.Elt(ind,i3) =reshape(full(NNode(i2+1)),size(i2,1),size(i2,2));
end % of jGroup loop

if ~isempty(i1)
 model=stack_set(model,'info','OrigNumbering',int32([model.Node(:,1) i1(:)]));
 model.Node(:,1)=i1;
 if length(unique(i1))~=length(i1) % If some nodes are duplicated eliminate
   i2=[]; i2(flipud(i1))=length(i1):-1:1; 
   model.Node=model.Node(i2(find(i2)),:);
 end
end

% deal with deformations

if isfield(model,'DOF')
  model.DOF=rem(model.DOF,1)+full(NNode(fix(model.DOF)+1));
end
if isfield(model,'def')&isfield(model.def,'DOF')
  model.def.DOF=rem(model.def.DOF,1)+full(NNode(fix(model.def.DOF)+1));
end

out=model;

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

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
% 'RepeatSel'      nITE tx ty tz Orig# AnleDeg nx ny nz ite0
if comstr(Cam,'sel')  [CAM,Cam]=comstr(CAM,4);

[carg,FEnode,FEel0,r1,ModelStack]=get_nodeelt(varargin,carg,ModelStack);

if carg<=nargin & all(size(varargin{carg})==4)
   i6=varargin{carg};carg=carg+1; opt(1)=-1;
else
  if ~isempty(strfind(Cam,'o'))
    i1=strfind(Cam,'o');
    opt=comstr(CAM([1:i1-1 i1+1:length(CAM)]),[-1 1 0 0 0  0 0 0   0   1 0 0]);
    i5=opt(5:7); opt=[opt(1:4) 0 opt(8:10)];
  else
    opt=comstr(CAM,[-1  1   0  0  0   0      0     1  0  0]);
    i3=find(FEnode(:,1)==opt(5));
    if isempty(i3) i5(1,1:3)=[0 0 0];
    else i5(1,1:3)=[FEnode(i3(1),5:7)];end
  end

  % Repeats element groups in FEel0
  elt=[];
  % i5 [origin;tx ty tz;rx ry rz];
  i5(2,1:3)=opt(2:4);
  i3=opt(7:9); if norm(i3)~=0 i3=i3/norm(i3); else opt(6)=0; end;i5(3,:)=i3;
  i5(4:6,1:3)=[0 -i5(3,3) i5(3,2);i5(3,3) 0 -i5(3,1);-i5(3,2) i5(3,1) 0];
  opt(6)=opt(6)/180*pi;
    if opt(6)~=0 %rotations included
        i6=[eye(3) i5(1,:)';0 0 0 1]* ...
           [cos(opt(6))*eye(3)+i5(3,:)'*i5(3,:)*(1-cos(opt(6)))+ ...
	    sin(opt(6))*i5(4:6,1:3) zeros(3,1);zeros(1,3) 1]* ...
           [eye(3) -i5(1,:)';0 0 0 1];
    else i6=eye(4); end
    i6 = i6*[eye(3) i5(2,:)';0 0 0 1]; %translations
end

[EGroup,nGroup]=getegroup(FEel0);
elt=[]; r3=zeros(1,nGroup); i0=0;

for jGroup=1:nGroup;

  ElemF= feutil('getelemf',FEel0(EGroup(jGroup),:),jGroup);
  elt = [elt;FEel0(EGroup(jGroup),:)];
  celle = FEel0(EGroup(jGroup):EGroup(jGroup+1)-1,:);

  % define celle node numbers in terms of indices in FEnode
  iNode=fe_super('node',ElemF);iNode = iNode(:)';
  NNode(FEnode(:,1))=[1:length(FEnode(:,1))]';

  % i1 indices in FEnode of all the needed nodes, i2 nodes without repeat
  i1 = celle(2:size(celle,1),iNode);
  i1 = reshape(NNode(i1(:)),size(i1,1),size(i1,2));
  i2 = find(sparse(i1(:),1,i1(:))); ci2(i2) = 1:length(i2);
  % loop on nITE
  j2 = 1:opt(1); if opt(1)==-1 j2 = 2; end % this is for translations
  % r1 contains current cell positions as a quadrivector
  r1=[FEnode(i2,5:7) ones(length(i2),1)]; if opt(1)==-1 r1 = (i6*r1')'; end
  if length(opt)>=10&opt(10)<0 % % move back if needed
    for j1=1:abs(opt(10));r1=(i6\r1')';end
  end  % move back if needed
  i0=i0(end)+[1:size(r1,1)];
  for j1 = j2-1
    [FEnode,i4]=feutil('AddNode',FEnode,r1(:,1:3));
    r3(i0,j1+1)=FEnode(i4,1); % set node group value to cell number
    i4=reshape(FEnode(i4(ci2(i1)),1),size(i1,1),size(i1,2));
    i5 = celle(2:size(celle,1),:); i5(:,iNode) = i4; 
    elt = [elt; i5];
    r1 = (i6*r1')'; % modified cell nodes
  end

end % jGroup of loop on groups
out=struct('Node',FEnode,'Elt',elt,'Stack',[]);
out.Stack={'info','CellNodes',unique(r3,'rows')};

else sdtw('Repeat%s unknown',CAM); end % subcommand selection - - - - - -

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

 [r1,out]=feutil(['find' CAM],varargin{2:end});

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

% 'SetGroup ...'
if comstr(Cam,'group') [CAM,Cam] = comstr(CAM,6);

[carg,FEnode,elt,el0,ModelStack]=get_nodeelt(varargin,carg,ModelStack);
if carg<=nargin st1=varargin{carg};carg=carg+1; else st1='elt';end
[EGroup,nGroup]=getegroup(elt);

RunOpt=struct('MatId',[],'ProId',[],'EGID',[]);
i2=[strfind(Cam,'mat') length(Cam)+1];i1(1)=i2(1);
i2=[strfind(Cam,'pro') strfind(Cam,'sec') length(Cam)+1];i1(2)=i2(1);
i2=[strfind(Cam,'eg') length(Cam)+1];i1(3)=i2(1);
i2=[strfind(Cam,'name') length(Cam)+1];i1(4)=i2(1);
i1(5)=length(Cam)+1;

% select group
if ~isletter(Cam(1)) opt=comstr(Cam(1:min(i1)-1),[-1]);else opt=[];end
if isempty(opt & min(i1)~=1)  % group selected by name 
  st = comstr(CAM(1:min(i1)-1),1);
  if strcmp(st,'all') opt=1:nGroup;
  else
   for jGroup=1:nGroup
      [ElemF]= feutil('getelemf',elt(EGroup(jGroup),:));
      if comstr(ElemF,st) opt(end+1)=jGroup; end
   end
  end
end
opt=fix(opt);opt=opt(find(opt<=nGroup));%contains selected groups
if isempty(opt) error(['SetGroup: not a group in ' st1]); end


i2=i1(4)+4:min([i1(find(i1>i1(1))) i1(end)])-1; % mat
if i1(4)<i1(5) % name
 [i5,i2,i3,i4]=sscanf(CAM(i1(4)+4:length(CAM)),'%s',1);
 CAM(i2)=' ';
else i5=[];end


i2=i1(1)+3:min([i1(find(i1>i1(1))) i1(end)])-1; % mat
if ~isempty(i2) 
 i3=comstr(comstr(Cam(i2),'id','%c'),[-1 0]);RunOpt.MatId=fix(i3(1));
else i3=0;end
i2=i1(2)+3:min([i1(find(i1>i1(2))) i1(end)])-1;% pro
if ~isempty(i2) 
 i4 = comstr(comstr(Cam(i2),'id','%c'),[-1 0]);RunOpt.ProId=fix(i4(1));
else i4=0;end
i2=i1(3)+2:min([i1(find(i1>i1(3))) i1(end)])-1;% egid
if ~isempty(i2) 
  i8=comstr(comstr(Cam(i2),'id','%c'),[-1 0]);RunOpt.EGID=fix(i8(1));
end

for jGroup = opt %loop on element groups
   [ElemF,i7]= feutil('getelemf',elt(EGroup(jGroup),:),jGroup);
   [i6,SEopt]=fe_super('prop',ElemF);i6=[i6(:);0;0;0];
   cEGI = [EGroup(jGroup)+1:EGroup(jGroup+1)-1]';
   if ~isempty(RunOpt.MatId);
    if i6(1); elt(cEGI,i6(1))=ones(size(cEGI))*RunOpt.MatId;
    else; warning(sprintf('MatID not defined for %s',ElemF));end
   end
   if ~isempty(RunOpt.ProId);
    if i6(2); elt(cEGI,i6(2))=ones(size(cEGI))*RunOpt.ProId; 
    elseif SEopt(1)==1; fesuper(sprintf('set %s ProID %i',ElemF,RunOpt.ProId))
    else; warning(sprintf('ProID not defined for %s',ElemF));end
   end
   if ~isempty(RunOpt.EGID) % egid is given
     i7(1)=RunOpt.EGID;elt(EGroup(jGroup),1:2+length(ElemF)+length(i7))=...
      [Inf abs(ElemF) 0 i7];
   end
   if ~isempty(i5)
    i7 = elt(EGroup(jGroup),2+length(ElemF):end);
    if isempty(i7)  elt(EGroup(jGroup),2:end)=0;
                    elt(EGroup(jGroup),2:length(i5)+1)=abs(i5);
    else elt(EGroup(jGroup),2:length(i5)+2+length(i7))=[abs(i5) 0 i7];
    end
   end
end

st=['For group(s) ' sprintf('%i ',opt)];
if ~isempty(RunOpt.MatId); st=[st sprintf('  MatID set to %i  ',i3)]; end
if ~isempty(RunOpt.ProId); st=[st sprintf('  ProID set to %i  ',i4)]; end
if ~isempty(RunOpt.EGID); st=[st sprintf('  EGID  set to %i  ',i8)]; end
if ~isempty(i5) st=[st sprintf('  Name set to %s  ',i5)]; end
out1=st;
out=elt;

else error('Not a valid command');end

% ------------------------------------------------------------------------
elseif comstr(Cam,'s')  CAM=comstr(CAM,'string','%c');Cam=comstr(CAM,-27);

%StringDof
if comstr(Cam,'d') CAM=comstr(CAM,'dof','%s');Cam=comstr(CAM,-27);

    for j1=99:-1:13 st{j1}=sprintf('.%i',j1);end
    st([1:12 19 20 21])={'x','y','z','\theta x','\theta y','\theta z', ...
        '-x','-y','-z','-\theta x','-\theta y','-\theta z','p','T','V'}';
    r1 = varargin{carg};carg=carg+1; r1=r1(:);
    if isa(r1,'cell')
      out=zeros(length(r1),1); for j1=1:length(r1)
        [i2,st1] = comstr(r1{j1},'','%i');
        i3=strmatch(st1,st); if isempty(i3)|i3(1)>100; i3=99;end
        out(j1,1)=i2+i3(1)/100;
      end
    else
      out=cell(length(r1),1); for j1=1:length(r1)
        i2=abs(r1(j1));i2 = [fix(i2) round(rem(i2,1)*100)];
        if r1(j1)<0
           out{j1}=sprintf('e%i.%3i',i2(1),round(rem(-r1(j1),1)*1000));
        elseif i2(1)<1&i2(2)>0; out{j1}=sprintf('%s',st{i2(2)});
        elseif i2(2)==0; out{j1}=sprintf('%i*',i2(1));
        elseif i2(1)<0;out{j1}=sprintf('e%i.%3i',abs(i2));
        else; out{j1}=sprintf('%i%s',i2(1),st{i2(2)}); end
      end
    end
% StringIODOF FLOR:180:+Z / FRNT:15:+Z
elseif comstr(Cam,'iodof')
 
else out='unknown'; end % subcommand selection - - - - - - - - - - - - - -

% ------------------------------------------------------------------------
elseif comstr(Cam,'trace2elt') CAM=comstr(CAM,'trace2elt','%s');Cam=comstr(CAM,-27);

   elt=varargin{carg};carg=carg+1;
   r1=[Inf abs('beam1') 0 -1];r2=[Inf abs('mass2') 0 -1]; 
   out1={};
   for j1 = 1:size(elt,1)
        i2=[elt(j1,83:end-1);elt(j1,84:end)];
        i2 = i2(:,find(all(i2)))';
        i4=[0 elt(j1,83:end) 0];
        i3=find(i4(1:end-2)==0&i4(2:end-1)~=0&i4(3:end)==0);
        if ~isempty(i3) i2=[i2;i4(i3+1)'*[1 1]]; end
        i3=find(i2(:,1)~=i2(:,2)); if ~isempty(i3)
          r1(end+[1:size(i3,1)],1:4)=[i2(i3,:) ones(size(i3))*[j1 elt(j1,2)]];
        end
        i3=find(i2(:,1)==i2(:,2)); if ~isempty(i3)
          r2(end+[1:size(i3,1)],[1 14 15])=[i2(i3,1) ...
                                        ones(size(i3))*[j1 elt(j1,2)]];
        end
        out1{j1,1}=deblank(char(elt(j1,3:82)));
   end; 

   if size(r1,1)>1 out = [r1]; else out=[]; end
   if size(r2,1)>1 if size(out,2)<15 out(1,15)=0; end;out = [out;r2]; end
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
elseif comstr(Cam,'cvs')

 out='$Revision: 1.142 $  $Date: 2006/05/16 17:26:30 $';

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
elseif comstr(Cam,'rmfield') % Remove fields from a structure if they exist

r1=varargin{carg}; carg=carg+1;
st=varargin{carg}; carg=carg+1;
if ~isa(st,'cell') st={st varargin{carg:end}};carg=nargin+1;end
st=intersect(fieldnames(r1),st);
if ~isempty(st); out=rmfield(r1,st);else; out=r1;end

% ------------------------------------------------------------------------
elseif exist('feutilb','file') 

 if nargout<2 out=feutilb(varargin{:});
 elseif nargout==2 [out,out1]=feutilb(varargin{:});
 else [out,out1,out2]=feutilb(varargin{:});end

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


% ------------------------------------------------------------------------
% ------------------------------------------------------------------------
% ------------------------------------------------------------------------
function [st,i1,carg]=ineqarg(Cam,carg,varargin);

   st=''; Cam=comstr(Cam,1);
   if ~isempty(Cam)&any(Cam(1)=='~<>=') st=Cam(1);Cam=Cam(2:end);end;
   if ~isempty(Cam)&any(Cam(1)=='=') st=[st Cam(1)];Cam=Cam(2:end);end;

   i1=comstr(Cam,-1);
   if isempty(i1)&carg<=nargin-2
    i1=varargin{carg};carg=carg+1;i1=i1(:)';
   elseif Cam(1)=='>' | Cam(1)=='<'
    st(2)=Cam(1);i1=comstr(Cam(2:end),-1);
   end

% ------------------------------------------------------------------------
% integer inequality tests
function [i4]=IEqTest(val,opts,opt);

   if ischar(opt) opt=comstr(opt,-1);end

   if comstr(opts,'~=') 
     %i0=-min(min(opt),min(val))+1; NNode(i0+opt)=1:length(opt);
     %val=val+i0; if max(val)>length(NNode) NNode(max(val))=0; end
     %i4=NNode(val);i4=find(i4==0);
     i4=find(~ismember(val,opt));
   elseif comstr(opts,'>=') i4=find(val>=opt(1));
   elseif comstr(opts,'<=') i4=find(val<=opt(1));
   elseif comstr(opts,'>')  i4=find(val>opt(1));
   elseif comstr(opts,'<')  i4=find(val<opt(1));
   else % ==
     %i0=-min(min(opt),min(val))+1;NNode(i0+opt)=1:length(opt);
     %val=val+i0;if max(val)>length(NNode) NNode(max(val))=0; end
     %i4=NNode(val); i4=find(i4);
     i4=find(ismember(val,opt));
   end

% ------------------------------------------------------------------------
% WHEN CHANGING CHANGE IN FEUTILB
function [carg,node,elt,el0,ModelStack]=get_nodeelt(varg,carg,ModelStack);

  node=varg{carg};carg=carg+1;
  if isstruct(node)&isfield(node,'Node')

   if isfield(node,'mdl'); assignin('caller','model',node.mdl);
    if isfield(node.mdl,'Stack'); ModelStack=node.mdl.Stack;end
   elseif isfield(node,'Stack'); ModelStack=node.Stack;
   end

   if ~isfield(node,'Elt') error('Fields .Node & .Elt must be defined');end
   elt=node.Elt; 
   if isfield(node,'El0') el0=node.El0; else el0=[];end
   if isfield(node,'bas')&~isempty(node.bas)
        [node,bas]=basis('nodebas',node); %node.Node,node.bas);
   else node=node.Node; end

  elseif ~isempty(node)&size(node,2)~=7
    error('Node matrices must have 7 columns');
  else
   elt = varg{carg};carg=carg+1;
   if nargout<4
   elseif carg<=length(varg)&(isempty(varg{carg})|~finite(varg{carg}(1)))
      el0 = varg{carg};carg=carg+1;
   else el0=[]; end
  end


% ------------------------------------------------------------------------
% return non repeated lines in a list
function [list,prop]=single_line(list,prop);

% eliminate degenerate lines
i1=find(list(:,1)==list(:,2));list(i1,:)=[];prop(i1,:)=[];

% flip lines to have lowest number node first

% fill in zeros
for j1=size(list,2)-1:-1:2 
 i1=find(~any(list(:,j1+1:end),2)&list(:,1)>list(:,j1));
 if ~isempty(i1) list(i1,1:j1)=list(i1,j1:-1:1);end
end

i3=1;
while ~isempty(i3)
   [r2,i1]=sortrows(sort(list,2));
   i2=1:size(list,1);i3=find(~any(diff(r2),2));
   i2(i3)=0;i2(i3+1)=0;i2=find(i2);
   list=list(i1(i2),:);prop=prop(i1(i2),:);
end


% ------------------------------------------------------------------------
% ------------------------------------------------------------------------
function [Stack,carg]=BuildSelStack(CAM,node,elt,el0,carg,varargin)

Args=varargin;

% end of edge selection - - - - - - - - - - - - - - - - - - - - - - - - - -

i7=[]; CAM=[CAM ' '];out=[]; opt=[];
[EGroup,nGroup]=getegroup(elt);

ind=[find(CAM=='&'|CAM=='|'|CAM=='{') length(CAM)+1]; j1=0;
Stack={};

% ind keeps track of what has been read
while j1<=length(CAM)  % loop on the selected arguments to build stack
try;
 [st,Cam]=comstr(CAM(j1+1:ind(1)-1),1); 
 if ind(1)<length(CAM)&CAM(ind(1))=='{'
  [st1,CAM(ind(1):end)]=comstr(CAM(ind(1):end),-8,'{}');
  ind = ind(find(ind>=ind(1)+length(st1)));
 else st1=''; end
 if j1 Stack{end+1,1}=CAM(j1); else Stack{1}='{'; end
 % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 if comstr(Cam,'between')  
   [opt,st,Cam]=comstr(Cam,'between','%i');
   if isempty(opt)&carg<=nargin & ~ischar(Args{carg})
    opt=Args{carg}; carg=carg+1;
   end
   Stack(end,2:4)={'between',[],opt};
 % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 elseif comstr(Cam,'egid')  [st,Cam]=comstr(st,5);
   [opts,opt,carg]=ineqarg(Cam,carg,Args{:});i1 = 5;
   if isempty(opt)&carg<=nargin & ~ischar(Args{carg})
    opt=Args{carg}; carg=carg+1;
   end
   Stack(end,2:4)={'EGID',opts,opt};
 % EltName - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 elseif comstr(Cam,'eltname')
   [opt,Cam]=comstr(st,8); i1 = 7; opts='';
   if comstr(opt,'~=') opts='~='; opt=comstr(opt,3);end
   Stack(end,2:4)={'EltName',opts,opt};

 % EltId - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 elseif comstr(Cam,'eltid') [st,Cam]=comstr(Cam,'eltid','%c');
   [opts,opt,carg]=ineqarg(Cam,carg,Args{:});i1 = 9;
   if isempty(opt)&carg<=nargin & ~ischar(Args{carg})
    opt=Args{carg}; carg=carg+1;
   end
   Stack(end,2:4)={'EltId',opts,opt};
 % EltInd - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 elseif comstr(Cam,'eltind') [st,Cam]=comstr(Cam,'eltind','%c');
   [opts,opt,carg]=ineqarg(Cam,carg,Args{:});i1 = 8;
   if isempty(opt)&carg<=nargin & ~ischar(Args{carg})
    opt=Args{carg}; carg=carg+1;
   end
   Stack(end,2:4)={'EltInd',opts,opt};

 % Facing- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 elseif comstr(Cam,'facing') [st,Cam]=comstr(Cam,'facing','%c');
   [opts,opt,carg]=ineqarg(Cam,carg,Args{:});
   if isempty(opt)&carg<=nargin & ~ischar(Args{carg})
    opt=Args{carg}; carg=carg+1;
   end
   Stack(end,2:4)={'Facing',opts,opt};

 % GID - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 elseif comstr(Cam,'gid') [st,Cam]=comstr(st,4);

   st1='GID'; 
   [opts,opt,carg]=ineqarg(Cam,carg,Args{:});
   if isempty(opt)&carg<=nargin & ~ischar(Args{carg})
    opt=Args{carg}; carg=carg+1;
   end
   Stack(end,2:4)={st1,opts,opt};

 % Group - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 elseif comstr(Cam,'group') [st,Cam]=comstr(st,6);

   st1='Group'; 
   if comstr(Cam,'al') opt=1:nGroup; opts='==';
   elseif comstr(Cam,'a') 
      st1='GroupA'; [opts,opt,carg]=ineqarg(Cam(2:end),carg,Args{:});
   else [opts,opt,carg]=ineqarg(Cam,carg,Args{:}); end
   if isempty(opt)&carg<=nargin & ~ischar(Args{carg})
    opt=Args{carg}; carg=carg+1;
   end
   Stack(end,2:4)={st1,opts,opt};

 elseif comstr(Cam,'in')  [st,Cam]=comstr(st,3);% - - - - - - - - - - - -
 % InNode - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  if comstr(Cam,'node')

   [st,Cam]=comstr(st,5); opt = round(comstr(st,[-1]));
   if     isempty(opt)&isempty(st1)&carg<=nargin & ~ischar(Args{carg})
     opt=Args{carg}; carg=carg+1;   Stack(end,2:4)={'InNode','',opt};
   elseif isempty(opt)&(~isempty(st1)|~isempty(st))
    if isempty(st1) st1=st; else st1=st1(2:end-1); end
    %opt=feutil(['findnode' st1(2:end-1)],node,elt,el0);
    Stack(end,2:4)={'InNode','{}',''};
    [Stack2,carg]=BuildSelStack(st1,node,elt,el0,carg,Args{:});
    Stack(end+[1:size(Stack2,1)],1:4)=Stack2;
   elseif ~isempty(opt)&isempty(st1)
    Stack(end,2:4)={'InNode','',opt};
   elseif ~isempty(opt)&~isempty(st1)
    error('InNode : Recursion is only valid with string argument');
   else warning('Incomplete InNode Selection');
   end

 % 'InElt' - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  elseif comstr(Cam,'elt')

   if isempty(st1) error('InElt requires specification of a selection');end
   if  Cam(end)=='t' Stack(end,2:3)={'inelt','{}'};
   elseif  Cam(end)=='0' Stack(end,2:3)={'inel0','{}'};
   else error('Not a valid FindNode InElt Command');
   end
   [Stack2,carg]=BuildSelStack(st1(2:end-1),node,elt,el0,carg,Args{:});
   Stack(end+[1:size(Stack2,1)],1:4)=Stack2;
  else
    error(sprintf('in %s : improper feutil selection command',Cam));
  end

 % MatId  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 elseif comstr(Cam,'mat')  [st,Cam]=comstr(Cam,'matid','%c');
   [opts,opt,carg]=ineqarg(Cam,carg,Args{:});i1 = 4;
   if isempty(opt)&carg<=nargin & ~ischar(Args{carg})
    opt=Args{carg}; carg=carg+1;
   end
   Stack(end,2:4)={'MatId',opts,opt};

 % NodeId - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 elseif comstr(Cam,'nodeid') [st,Cam]=comstr(st,7);

   [st,opt,carg]=ineqarg(Cam,carg,Args{:});
   Stack(end,2:4)={'nodeid',st,opt};

 % 'NotIn' - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  elseif comstr(Cam,'notin')

   if isempty(st1) error('NotIn requires specification of a selection');end
   Stack(end,2:3)={'notin','{}'};
   [Stack2,carg]=BuildSelStack(st1(2:end-1),node,elt,el0,carg,Args{:});
   Stack(end+[1:size(Stack2,1)],1:4)=Stack2;

 % Plane - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 elseif comstr(Cam,'plane') [st,Cam]=comstr(st,6);

  st1='==';if any(Cam(1)=='~<>=') st1=Cam(1);Cam=Cam(2:end);end;
  if ~isempty(Cam)&any(Cam(1)=='=') st1=[st1 Cam(1)];Cam=Cam(2:end);end;
  [st,Cam]=comstr(Cam,1);
  if comstr(Cam,'o')    opt=comstr(Cam(2:length(Cam)),[-1 0 0 0 0 0 0]);
  elseif isempty(Cam)&carg<=nargin  opt=Args{carg};carg=carg+1;
  else    opt =  comstr(Cam,[-1 0 0 0 0]);i1=find(node(:,1)==opt(1));
          if isempty(i1) i1=[0 0 0];else i1=node(i1(1),[5:7]); end
          opt=[i1 opt(2:4)];
  end

  Stack(end,2:4)={'plane',st1,opt};

 % ProId  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 elseif comstr(Cam,'pro') [st,Cam]=comstr(Cam,'proid','%c');
   [opts,opt,carg]=ineqarg(Cam,carg,Args{:});i1 = 3;
   if isempty(opt)&carg<=nargin & ~ischar(Args{carg})
    opt=Args{carg}; carg=carg+1;
   end
   Stack(end,2:4)={'ProId',opts,opt};
 % Rad - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 elseif comstr(Cam,'rad') [st,Cam]=comstr(st,4); st='';

   if any(Cam(1)=='~<>=') st=Cam(1);Cam=Cam(2:end);end;
   if ~isempty(Cam)&any(Cam(1)=='=') st=[st Cam(1)];Cam=Cam(2:end);end;

   if isempty(Cam)&carg<=nargin
    i1=Args{carg};carg=carg+1;i1=i1(:)'; i2 = length(i1);
   else [i1,i2,i3,i6]=sscanf(Cam,'%g %g %g %g');i1=i1'; end
   Stack(end,2:4)={'rad',st,i1};

 % SelEdge - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 elseif comstr(Cam,'seledge')
  [st,Cam]=comstr(Cam,'seledge','%c');
  Stack(end,2:4)={'SelEdge',[],st};
 % SelFace  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 elseif comstr(Cam,'selface')
  [st,Cam]=comstr(Cam,'selface','%c');
  Stack(end,2:4)={'SelFace',[],st};
 % SelName  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 elseif comstr(Cam,'selname'); 

   Cam(1:7)='';[st,i1,i2,i3]=sscanf(Cam,'%s',1);Cam=Cam(i3:end);
   if st(1)=='"'; st(1)=[];% with blanks
    i2=min(find(st=='"'));
    if ~isempty(i2); 
     st=st(1:i2-1);Cam=comstr(horzcat(st(i2+1:end),Cam),1);
     if ~isempty(Cam)&Cam(1)=='"';Cam=comstr(Cam(2:end),1);end
    else
     i2=min(find(Cam=='"')); st=horzcat(st(1:end),Cam(1:i2-1));
     Cam=comstr(Cam(i2+1:end),1);
    end
   end
   Stack(end,2:4)={'Sel',[],st};
 % SetName  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 elseif comstr(Cam,'setname'); 

   Cam(1:7)='';[st,i1,i2,i3]=sscanf(Cam,'%s',1);Cam=Cam(i3:end);
   if st(1)=='"'; st(1)=[];% with blanks
    i2=min(find(st=='"'));
    if ~isempty(i2); 
     st=st(1:i2-1);Cam=comstr(horzcat(st(i2+1:end),Cam),1);
     if ~isempty(Cam)&Cam(1)=='"';Cam=comstr(Cam(2:end),1);end
    else
     i2=min(find(Cam=='"')); st=horzcat(st(1:end),Cam(1:i2-1));
     Cam=comstr(Cam(i2+1:end),1);
    end
   end
   Stack(end,2:4)={'Set',[],st};

 % Set  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 elseif comstr(Cam,'set') [st,Cam]=comstr(Cam,'set','%c');
   [opts,opt,carg]=ineqarg(Cam,carg,Args{:});i1 = 3;
   if isempty(opt)&carg<=nargin & ~ischar(Args{carg})
    opt=Args{carg}; carg=carg+1;
   end
   Stack(end,2:4)={'Set',opts,opt};

 % WithoutNode - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 elseif comstr(Cam,'without') [st,Cam]=comstr(st,8); % - - - - - - - - - -
  if ~comstr(Cam,'node')
    error(sprintf('without %s : improper feutil(''findelt'') command',Cam));
  end
  [st,Cam]=comstr(st,5); opt = round(comstr(st,[-1]));   
  if     isempty(opt)&isempty(st1)&carg<=nargin & ~ischar(Args{carg})
    opt=Args{carg}; carg=carg+1;Stack(end,2:4)={'WithoutNode','',opt};
  elseif isempty(opt)&~isempty(st1)
   Stack(end,2:4)={'WithoutNode','{}',''};
   [Stack2,carg]=BuildSelStack(st1(2:end-1),node,elt,el0,carg,Args{:});
   Stack(end+[1:size(Stack2,1)],1:4)=Stack2;
  elseif ~isempty(opt)&isempty(st1)
   Stack(end,2:4)={'WithoutNode','',opt};
  elseif ~isempty(opt)&~isempty(st1)
   error('WithoutNode : Recursion is only valid with string argument');
  end

 elseif comstr(Cam,'with')  [st,Cam]=comstr(st,5); % - - - - - - - - - -
  if ~comstr(Cam,'node')
    error(sprintf('with %s : improper feutil(''findelt'') command',Cam));
  end
  [st,Cam]=comstr(st,5); opt = round(comstr(st,[-1])); 
  if     isempty(opt)&isempty(st1)&carg<=nargin & ~ischar(Args{carg})
    opt=Args{carg}; carg=carg+1;    Stack(end,2:4)={'WithNode','',opt};
  elseif isempty(opt)&~isempty(st1)
   %opt=feutil(['findnode' st1(2:end-1)],node,elt,el0);
   Stack(end,2:4)={'WithNode','{}',''};
   [Stack2,carg]=BuildSelStack(st1(2:end-1),node,elt,el0,carg,Args{:});
   Stack(end+[1:size(Stack2,1)],1:4)=Stack2;
  elseif ~isempty(opt)&isempty(st1)
   Stack(end,2:4)={'WithNode','',opt};
  elseif ~isempty(opt)&~isempty(st1)
   error('WithNode : Recursion is only valid with string argument');
  end

 % xyzr - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 elseif ~isempty(Cam) & any(Cam(1)=='xyzr') % FindNode xyz

   i2=find('xyzr'==Cam(1)); st1='xyzr';   [st,Cam]=comstr(st,2); st='';
   if any(Cam(1)=='~<>=') st=Cam(1);Cam=Cam(2:end);end;
   if ~isempty(Cam)&any(Cam(1)=='=') st=[st Cam(1)];Cam=Cam(2:end);end;
   i1=comstr(Cam,[-1]); 

   if isempty(i1)&~isempty(Cam); i1=Cam;
   elseif isempty(i1)&carg<=nargin  i1=Args{carg};carg=carg+1;
   end
   if isempty(i1) error('xyzr ><= ... must specify value'); end
   if ~ischar(i1); i1=i1(1);end 
   Stack(end,2:4)={st1(i2),st,i1};

 % give coordinates  with possible infinity - - - - - - - - - - - - - - - -
 elseif ~isempty(Cam)&~all(Cam==' ') 

   if isempty(st)&carg<=nargin; opt=Args{carg};carg=carg+1;
     Stack(end,[2 4])={'Position',opt};

   elseif ~isempty(st) &any(comstr(st(1),-27)=='0123456789.I+-') 
     opt=comstr(st,[-1]);
     if length(opt)==3|any(rem(opt,1)~=0);
       opt(end+1:3)=Inf;Stack(end,[2 4])={'Position',opt};
     else
       Stack(end,[2 3 4])={'NodeId','==',opt(:)};
     end
   else error(sprintf('''%s'' is not a valid FindNode/FindElt command',st));
   end
   

 end

 catch
  if j1<2
   error(sprintf('''%s ...'' is not a valid selection format', ...
     CAM(1:min(j1+20,size(CAM,2)))));
  else
   error(sprintf(''' ... %s'' is not a valid selection format', ...
     CAM(j1:min(j1+20,size(CAM,2)))));
  end
 end

 j1=ind(1);ind=ind(2:end); 

end % of j1 loop on arguments
Stack{end+1,1}='}';

% ------------------------------------------------------------------------
% ------------------------------------------------------------------------
function ind=SubStackIndex(Stack,i0);

   i1=find(strcmp(Stack(:,1),'{'));i2=find(strcmp(Stack(:,1),'}'));

   i4=[1 i0]; 
   while i4(1)~=0 
       i3=min(i1(find(i1>i4(2))));
       i5=min(i2(find(i2>i4(2))));
       if i3<i5 % opening first
         i4(1)=i4(1)+1; i4(2)=i3;
       else % i5 is closing
         i4(1)=i4(1)-1; i4(2)=i5;
       end
   end
   ind=i0:i4(2);

% ------------------------------------------------------------------------
% ------------------------------------------------------------------------
function [epsl,CAM,Cam]=test_epsl(epsl,CAM,Cam);

if isempty(epsl);try;epsl=sp_util('epsl');catch;epsl=sdtdef('epsl'); end;end 
i1=strfind(Cam,'epsl'); if ~isempty(i1);
  [epsl,i2,i3,i4]=sscanf(Cam(i1+4:end),'%g');
  CAM(i1+[0:i4+2])=''; [CAM,Cam] = comstr(CAM,1);
end

% ------------------------------------------------------------------------
% ------------------------------------------------------------------------
% Sets are documented in doc('model'), one supports SelN and Set
function i4=FeutilMatchSet(ModelStack,Stack,elt);


j1=1;i4=[];

if isnumeric(Stack{j1,4}) % Selection by ID number in SelV

 % Build the list of ID values
 i1=zeros(size(ModelStack,1),2);
 for j2=1:size(ModelStack,1)
   try; i1(j2)=ModelStack{j2,3}.ID(1);end
   i2=strmatch(comstr(ModelStack{j2,1},-27),{'set','seln','sele'},'exact');
   if ~isempty(i2) i1(j2,2)=i2(1);end
 end
 i2=find(i1(:,2)==1&i1(:,1));
 opt=IEqTest(i1(i2,1),Stack{j1,3},Stack{j1,4});

 [i4,opt]=intersect(i1(i2,1),opt);opt=i2(opt);
elseif ischar(Stack{j1,4});
 opt=strmatch(comstr(Stack{j1,4},-27),comstr(ModelStack(:,2),-27));
 if isempty(opt); error('Name is not a set'); end
 if ~any(strcmp(comstr(ModelStack{opt,1},-27),{'set','seln','sele'})); 
  error('Name is not a set');
 end
end
if isempty(opt); error('Set selection failed');end

 if nargin==2 % Match a set pointing to nodes
  i4=[]; for j2=1:length(opt)
    r1=ModelStack{opt(j2),3}; if isfield(r1,'data');r1=r1.data(:,1);end
    i4=[i4;r1];
  end
 elseif length(opt)==1&ischar(ModelStack{opt,3}) % char selection

  node=evalin('caller','node');
  if comstr(ModelStack(opt,1),'seln'); % THIS DOES NOT WORK ?
   i4=feutil(horzcat('findnode',ModelStack{opt,3}),node,elt);
  else % 
   i4=feutil(horzcat('selelt',ModelStack{opt,3}),node,elt);
  end

 else % Matching faces based on set information

 EltId=feutil('eltidfix;',elt);assignin('caller','EltId',EltId);
 [EGroup,nGroup]=getegroup(elt);
 i4=[];
 for jGroup=1:nGroup
  [ElemF,i1,ElemP]= feutil('getelemf',elt(EGroup(jGroup),:),jGroup);
  if any(strcmp(ElemF, ...  % 2-D elements
    {'q4p','t3p','q5p','q8p','q8a','q9a','t3a','t3p','t6a','t6p'}))
   FaceIndex=feval(ElemP,'edge')';
   if size(FaceIndex,1)==2 i4(end+1,1:6)=[Inf abs('beam1')];
   else; i4(end+1,1:6)=[Inf abs('beam3')];end
  else FaceIndex=feval(ElemP,'face')';
   if size(FaceIndex,1)==4 i4(end+1,1:6)=[Inf abs('quad4')];
   else; error('Not implemented');
   end
  end
  if length(opt)>1 error('Multiple set selection not implemented');end


  nind=sparse(EltId+1,1,1:size(elt,1));
  ind=full(nind(ModelStack{opt,3}.data(:,1)+1));

  r1=ind(:,ones(size(FaceIndex,1),1))+ ...
   (FaceIndex(:,ModelStack{opt,3}.data(:,2))-1)'*size(elt,1);

  i4(end+[1:size(r1,1)],1:size(r1,2))=elt(r1);
 end

  
 end




%-------------------------------------------%
%-- Additionnal functions to check meshes --%
%-- Mathieu Corus -- %
%-------------------------------------------%

%-- Check Penta6 / Hexa8 element --%

function t_node=elt_check(t_node,elt_type,model,NNode,bypass);

switch comstr(elt_type,-27)
  case 'hexa8'
    
    t_node3=unique(t_node);
    %-- check if hexa8 is correct --%
    rls=[0 0 0 1
      1 0 0 1
      0 1 0 1
      1 1 0 1
      0 0 1 1
      1 0 1 1
      0 1 1 1
      1 1 1 1
      0.5 0.5 0 1
      0.5 0.5 1 1
      0.5 0 0.5 1
      0.5 1 0.5 1
      0 .5 0.5 1
      1 0.5 0.5 1];
    rls(:,1:3)=rls(:,1:3)*2;
    rls(:,1:3)=rls(:,1:3)-ones(14,3);
    rules=integrules('hexa8',rls); % might want to add points
    of_mk('buildndn',3,rules,model.Node(NNode(t_node),5:7));
    if bypass==0    
      if all(rules.jdet>0) return; end; %-- if hexa8 is correct, do nothing --%
    end;
    disp(['Processing Hexa8 with nodes ' num2str(t_node(:)') ' ...']);
    pause(.01)
    %-- if not, find the best hexa8 possible --%
    
    % get convex hull for point cloud --%
    [hull,v]=convhulln(model.Node(NNode(t_node),5:7));
    vec1=model.Node(NNode(t_node(hull(:,2))),5:7)-model.Node(NNode(t_node(hull(:,1))),5:7);
    vec2=model.Node(NNode(t_node(hull(:,3))),5:7)-model.Node(NNode(t_node(hull(:,2))),5:7);
    nor=cross(vec1,vec2); for i1=1:length(nor) nor(i1,:)=nor(i1,:)./norm(nor(i1,:)); end;
    
    %-- True Hexa8 element --%
    %-- not degenerated element --%
    if length(t_node3)==8
      % get the planest "face" 
      tt=unique(hull(:)); if length(tt)<8 disp(['Warning : Hexa8 element with nodes ' num2str(t_node(:)') ...
            ' is not convex']); end; 
      ori=[];
      for i1=1:length(nor)-1 ori(i1,1:length(nor)-i1)=nor(i1,:)*nor(i1+1:end,:)'; end;
      [val,x1]=max(ori); 
      max_val=find(val==max(val));
      t_node2=[]; deter=[];
      for k1=1:length(max_val) %-- if maximum is not unique, get best hexa8 --%                   
        
        [val,x1]=max(ori);
        y1=max_val(k1);x1=x1(y1);y1=y1+x1;
        face1=hull(x1,:)';face2=hull(y1,:)';
        n1=fe_c(face1,face2,'dof',2);n2=fe_c(face2,face1,'dof',2);
        if length(n1) == 1 & length(n2) == 1  
          face=[fe_c(face1,face2,'dof',1)' n1 n2];face=face([1 3 2 4]);  
          
          %-- check face orientation --%                   
          oth_no=fe_c([1:8]',face(:),'dof',2);
          mid_node_f=sum(model.Node(NNode(t_node(face)),5:7))/size(face(:),1);
          mid_node_o=sum(model.Node(NNode(t_node(oth_no)),5:7))/size(oth_no(:),1);
          ve=model.Node(NNode(t_node(face([2 4]))),5:7)-...
            model.Node(NNode(t_node(face([1 1]))),5:7);
          if cross(ve(1,:),ve(2,:))*(mid_node_o-mid_node_f)' <= 0 face=face(end:-1:1); end;
          
          %-- find opposite face --%
          vec1=model.Node(NNode(t_node(oth_no)),5:7)-...
            model.Node(NNode(t_node(face(ones(length(oth_no),1)))),5:7);
          for i1=1:length(vec1) vec1(i1,:)=vec1(i1,:)./norm(vec1(i1,:)); end;
          ind=[];                
          for i1=1:length(hull)
            if isempty(fe_c(face',hull(i1,:)','ind')) ind=[ind i1]; end;
          end;
          
          if length(ind)==2
            face1=hull(ind(1),:)';face2=hull(ind(2),:)';
            n1=fe_c(face1,face2,'dof',2);n2=fe_c(face2,face1,'dof',2);
            face2=[fe_c(face1,face2,'dof',1)' n1 n2];face2=face2([4 2 3 1]);
            
            %-- check face orientation --%                   
            oth_no=fe_c([1:8]',face2(:),'dof',2);
            mid_node_f=sum(model.Node(NNode(t_node(face2)),5:7))/size(face2(:),1);
            mid_node_o=sum(model.Node(NNode(t_node(oth_no)),5:7))/size(oth_no(:),1);
            ve=model.Node(NNode(t_node(face2([2 4]))),5:7)-...
              model.Node(NNode(t_node(face2([1 1]))),5:7);
            if cross(ve(1,:),ve(2,:))*(mid_node_o-mid_node_f)' >= 0 face2=face2(end:-1:1); end;
            
            deter=[];
            for i1=1:4
              rules=integrules('hexa8',rls);
              of_mk('buildndn',3,rules,...
                model.Node(NNode(t_node([face face2([i1:end 1:i1-1])])),5:7));
              if ~all(rules.jdet>=0)  deter=[deter -1]; 
              else deter=[deter sum(rules.jdet)]; end;
            end;
            
            [val,ind]=max(deter);
            t_node2(k1,:)=t_node([face face2([ind:end 1:ind-1])]);
          end; %-- if length(ind)==2
        end; %-- if length(n1) == 1 & length(n2) == 1
      end; %--for k1=1:length(max_val)
      
      t_node2=t_node2(find(sum(abs(t_node2'))),:);
      deter=[];
      for i1=1:size(t_node2,1)
        rules=integrules('hexa8',rls);
        of_mk('buildndn',3,rules,...
          model.Node(NNode(t_node2(i1,:)),5:7));
        if ~all(rules.jdet>=0)  deter=[deter -1]; else deter=[deter sum(rules.jdet)]; end;
      end;
      [val,ind]=max(deter); 
      if val > 0 t_node=t_node2(ind,:); else disp(['Warning : failed to find adequate Hexa8 with nodes '  ...
            num2str(t_node(:)')]); 
        disp(['     -> Try to build Penta6 elements instead']);
        disp(' ');
      end;
      
      
      
    elseif length(t_node3)==7 disp(['Warning : 7 nodes degenerated Hexa8 is not supported - element with nodes '...
          num2str(t_node(:)')]);
    elseif length(t_node3)==6 %-- assemble Penta6 as degenerated Hexa8 --%
      t_node=hex_pen_check('penta6',t_node3,model,NNode);
      t_node=t_node([1 2 3 3 4 5 6 6]);
    else error('Not supported yet')
    end; %--length(t_node3)==8
    
  case 'penta6'
    
    t_node3=unique(t_node);    
    %-- check if penta6 is correct --%
    alt_node=[];
    if length(t_node)==8
      alt_node=t_node;
    elseif length(t_node)==6
      alt_node=t_node([1 2 3 3 4 5 6 6]);
    end;
    if ~isempty(alt_node)
      rules=penta6b('constants',[],[],[]);
      of_mk('buildndn',3,rules,model.Node(NNode(alt_node),5:7));
      if all(rules.jdet>=0) return; end; %-- if penta6 is correct, do nothing --%
    end;
    %-- if not, find the best hexa8 possible --%
    % get convex hull for point cloud --%
    hull=convhulln(model.Node(NNode(t_node3),5:7));
    vec1=model.Node(NNode(t_node3(hull(:,2))),5:7)-model.Node(NNode(t_node3(hull(:,1))),5:7);
    vec2=model.Node(NNode(t_node3(hull(:,3))),5:7)-model.Node(NNode(t_node3(hull(:,2))),5:7);
    nor=cross(vec1,vec2); for i1=1:length(nor) nor(i1,:)=nor(i1,:)./norm(nor(i1,:)); end;
    
    % get the best opposite faces --% 
    tt=unique(hull(:)); if length(tt)~=6 error(['Penta6 element with nodes ' num2str(t_node(:)') ...
          'is not valid']); end; 
    ori=[];
    for i1=1:length(nor)-1 ori(i1,1:length(nor)-i1)=nor(i1,:)*nor(i1+1:end,:)'; end;
    [val,x1]=min(ori); 
    %-- check for all "opposite" faces --%
    min_val=find(val<0);
    t_node2=[];
    
    for k1=1:length(min_val) %-- if maximum is not unique, get best hexa8 --%                   
      
      [val,x1]=min(ori);
      y1=min_val(k1);x1=x1(y1);y1=y1+x1;
      face=hull(x1,:)';face2=hull(y1,:)';
      
      %-- check faces orientation --%    
      
      oth_no=fe_c([1:6]',face(:),'dof',2);
      mid_node_f=sum(model.Node(NNode(t_node3(face)),5:7))/size(face(:),1);
      mid_node_o=sum(model.Node(NNode(t_node3(oth_no)),5:7))/size(oth_no(:),1);
      ve=model.Node(NNode(t_node3(face([2 3]))),5:7)-...
        model.Node(NNode(t_node3(face([1 1]))),5:7);
      if cross(ve(1,:),ve(2,:))*(mid_node_o-mid_node_f)' <= 0 face=face(end:-1:1); end;
      face=face([1 2 3 3])';
      
      oth_no=fe_c([1:6]',face2(:),'dof',2);
      mid_node_f=sum(model.Node(NNode(t_node3(face2)),5:7))/size(face2(:),1);
      mid_node_o=sum(model.Node(NNode(t_node3(oth_no)),5:7))/size(oth_no(:),1);
      ve=model.Node(NNode(t_node3(face2([2 3]))),5:7)-...
        model.Node(NNode(t_node3(face2([1 1]))),5:7);
      if cross(ve(1,:),ve(2,:))*(mid_node_o-mid_node_f)' >= 0 face2=face2(end:-1:1); end;
      
      deter=[];faceb=[];
      for i1=1:3
        if i1==1 faceb(i1,:)=face2([i1:end end])'; 
        else faceb(i1,:)=face2([i1:end 1:i1-1 i1-1])';
        end;
        
        rules=penta6b('constants',[],[],[]);
        of_mk('buildndn',3,rules,...
          model.Node(NNode(t_node3([face faceb(i1,:)])),5:7));
        if ~all(rules.jdet>=0)  deter=[deter -1]; 
        else deter=[deter sum(rules.jdet)]; 
        end;
      end;
      [val,ind]=max(deter);
      t_node2(k1,:)=t_node3([face faceb(ind,:)]);
    end; %--k1=1:length(min_val)
    
    t_node2=t_node2(find(sum(abs(t_node2'))),:);
    deter=[];
    for i1=1:size(t_node2,1)
      rules=penta6b('constants',[],[],[]);
      of_mk('buildndn',3,rules,...
        model.Node(NNode(t_node2(i1,:)),5:7));
      if ~all(rules.jdet>=0)  deter=[deter -1]; else deter=[deter sum(rules.jdet)]; end;
    end;
    [val,ind]=max(deter); 
    if val > 0 t_node=t_node2(ind,:); else warning(['Failed to find adequate Penta6 with nodes '  ...
          num2str(t_node(:)')]);
    end;
    t_node=t_node([1 2 3 5 6 7]);
  case 'quad4'
    t_node4=unique(t_node);
    if length(t_node4)<4 return; end;
    t_node4=unique(t_node);
    sens_rot=0;  %-- check convexity --%
    compt=0;
    vec=model.Node([1 1 1]*NNode(t_node(1)),5:7)-model.Node(NNode(t_node(2:4)),5:7);
    si=[cross(vec(1,:),vec(2,:))*cross(vec(2,:),vec(3,:))'];
    ve=diff(model.Node(NNode(t_node([1:end 1])),5:7),1,1);
    ve=[ve;ve(1,:)];for i1=1:size(ve,1) ve(i1,:)=ve(i1,:)./norm(ve(i1,:)); end;
    for k1=1:4 
      n(k1,:)=cross(ve(k1,:),ve(k1+1,:)); 
    end;
    nn=n*n';
    ii=mod(find(nn(1,:)<0)+1,4); %-- find inflexion point

    if length(ii)==1  %-- non convex element
        if ii==0 ii=ii+4; end;
        disp(['Warning : Quad4 element with nodes ' num2str(t_node(:)') ' is not convex']);
        disp('      -> It will be divided into two degenerated quad4 elements')
        t_node=t_node([ii:4 1:(ii-1)]);
        t_node=[[t_node(1:3) t_node(3)];[t_node([1 3 4]) t_node(4)]];
    else
      while sens_rot==0
        vec=model.Node([1 1 1]*NNode(t_node(1)),5:7)-model.Node(NNode(t_node(2:4)),5:7);
        si=[cross(vec(1,:),vec(2,:))*cross(vec(2,:),vec(3,:))'];
        if si>0 sens_rot=1; 
        elseif si<0 t_node=t_node(randperm(4)); 
        elseif si==0 t_node=t_node([2 3 4 1]);
        end;
        compt=compt+1;
        if compt>=1000 
          disp(['Warning : failed to generate a proper Quad4 element with nodes ' ...
              num2str(t_node(:)')]);break; 
        end;
      end;
    end;

  end; %-- switch lower(elt_type)
  
%-- Check Quad4 / Penta6 / Hexa8 OK --%
