function [out,out1]=femesh(varargin);

%FEMESH	GUI commands for finite element mesh generation
%
%	Syntax: femesh
%		femesh CommandString
%		femesh('Command', ... )
%
%	The following a gives a summary of available commands. Refer to the
%       HTML manual for details (type doc('femesh') at the MATLAB prompt).
%
%	; (command chaining), q (exit command mode) (see COMMODE for details)
%
%	Add [FEel(i) FEel(j), Sel]      : combine models
%	Add [Node, NodeNew]             : add new nodes
%	Add Test (NodeShift)    : femesh('addtest (n0)',node,ldraw)
%	Divide [(div1 div2 div3), InGroups, Group (i)]
%       DispFlag [on/off]  : level of infos returned
%	Extrude (nRep tx ty tz)
%	Find Dof (Element Selectors) : see FindElt for selectors
%	FindElt [[with,in] (node), mat (i), pro (i), group (i), egid(i)]
%		+ possibility to use &,|,{} and trailing numeric arguments
%	Find Node [Group (i), Groupa (i), Plane (Orig# nx ny nz), (x y z), ...
%                 rad <= r x y z, [x,y,z][>,<,==,>=,<=](i)] 
%               + possibility to combine with &| and trailing numeric arguments
%	Get Node [...] same as Find but returns nodes rather than node numbers
%	Info [ , FEel(i), node(i)]
%	Join [ElemtName, Group (i)]
%       Model  return model structure
%	ObjectBeamLine (i)
%	ObjectHoleInPlate CtrN# Edge1N# Edge2N# r1 r2 nDiv1 nDiv2 nQuadrant
%	ObjectMass m xxx
%	Optim [Model, Nodenum]
%	Orient (i) n (nx ny nz) 
%	Plot [Elt,El0]
%	quad2tria, quad42quadb, hexa82hexa20, tria32tria6, penta62penta15
%                hexa2tetra, hexa2penta
%	RemoveEl[t,0] (FindElt Selectors)
%	RepeatSel (nITE tx ty tz Orig# AnleDeg nx ny nz)
%	Refine beam (l)
%	Rev (nRep Orig# AngleDeg nx ny nz tx ty tz)
%	RotateSel (Orig# AnleDeg nx ny nz)
%	Sel Edge [,{El0 selection}] [Line, Patch] [,(element selectors)]
%	SelElt (calls FindElt and puts result in FEel0)
%	SelGroup (i)
%	SetGroup[ ,a] (i) [Mat (i), Pro (i), Name (s)]
%	SymSel (Orig# nx ny nz)
%	TransSel (tx ty tz)
%	UnJoin (i j) or femesh('unjoin','EltSel1','EltSel2')
%
%	Some commands can be used with FEMESH input arguments as detailed in
%	the manual.
%
%	FEMESH uses a number of global variables
%	  FEnode/FEn0/FEn1    main/selected/alternate node set
%	  FEelt/FEel0/FEel1   main/selected/alternate model description matrix
%       FEMESH automatically checks that these variables are declared as
%	global in your base and caller workspaces.
%
%       See also help   feplot, fecom, fe_mk, feutil. 
%                doc    fem, femesh
%		 demos  gartfe, d_truss, d_ubeam, ...

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

persistent FE FEInfo DispFlag
global FEnode FEn0 FEn1 FEelt FEel0 FEel1

if nargin==1 & comstr(varargin{1},'cvs')
 out='$Revision: 1.112 $  $Date: 2006/05/19 12:09:41 $'; return;
end

epsl=sdtdef('epsl');
if isempty(DispFlag) DispFlag=1; end

% Provide a model to fill in the global variables 
if nargin>0 & ~ischar(varargin{1})
   
     evalin('caller', ...
        'iigui({''FEnode'',''FEn0'',''FEn1'',''FEelt'',''FEel0'',''FEel1''},''caller'');');
   r1=varargin{1};
   if isfield(r1,'Node')&isfield(r1,'Elt')
     FE=r1;
     FEnode=FE.Node;
     FEelt=FE.Elt;
   else warning('Input argument to femesh must be a model');
   end
   FEn0=[];FEn1=[];FEel0=[]; FEel1=[]; carg=2;
   st={'bas','Stack','pl','il'};

% Rather obsolete INIT call (never documented)
elseif nargin==1&ischar(varargin{1})&strcmp(varargin{1},'init')

  warning('This call is obsolete and will be removed in the future');
  [CAM,Cam]=comstr(varargin{1},5); if isempty(CAM) CAM='FE';end
  if evalin('base',sprintf('exist(''%s'',''var'')',CAM))
   FE=evalin('base',CAM);
  elseif ~isfield(FE,'BaseVar')
   FE=struct('Node',[],'Elt',[],'bas',[],'pl',[],'il',[]);
  end
  FE.BaseVar=CAM;
  assignin('base',FE.BaseVar,FE);
  return;

% Just return the model 
elseif nargin==0 & nargout==1

  out=FE;  out.Elt=FEelt; 
  if ~isempty(FEel0);out.El0=FEel0;
  else; out=feutil('rmfield',out,'El0');
  end
  out.Node=FEnode; return;

else
   carg=1;
   if isfield(FE,'BaseVar') % new format with persistent variable
     FE=evalin('base',FE.BaseVar);
     FEnode=FE.Node;FEelt=FE.Elt;
     if isfield(FE,'El0') FEel0=FE.El0; else FEel0=[]; end
     if isfield(FE,'N0') FEn0=FE.N0; else FEn0=[]; end
   else  % Standard attempt to use the local values of caller as global
    if sdtdef('isdeployed'); st='commode';
    else;  
      r1=dbstack;if length(r1)>1; [wd,st]=fileparts(r1(2).name); else st='';end
    end
    if nargin==0 | (~comstr(st,'commode'))
     evalin('caller', ...
        'iigui({''FEnode'',''FEn0'',''FEn1'',''FEelt'',''FEel0'',''FEel1''},''caller'');');
    end
   end
end  % How to deal with global(s)

if nargin<carg & nargout==0 
 return; 
elseif isstruct(varargin{carg})
elseif ischar(varargin{carg}) 
 [CAM,Cam]=comstr(varargin{carg},1);carg=carg+1;
 if ~comstr(Cam,'dispflag') & (DispFlag==1)
  DispFlag=varargin{nargin}; if ischar(DispFlag)&DispFlag(end)==';' DispFlag=0;
  elseif CAM(end)==';'  DispFlag=0; else DispFlag=1; end
 end
end

if ~isempty(FEel0)&finite(FEel0(1)) FEel0=[];end

out='done'; yPlot=0; CAM0=CAM;

% ------------------------------------------------------------------------
if comstr(Cam,';')

  commode('femesh',CAM(2:length(CAM)))

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

if isempty(Cam) out=DispFlag; return; end
if      comstr(Cam,'on')  DispFlag=1.;
elseif  comstr(Cam,'off') DispFlag=0.;
else    error('Command required: outputinfo on/off');
end

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

% ('addnode',OldNode,NewNode) - - - - - - - - - - - - - - - - - - -
if comstr(Cam,'node')

 if carg<=nargin r1=varargin{carg};carg=carg+1;
  if carg<=nargin r2=varargin{carg};carg=carg+1; end
 else 
  if DispFlag; disp('Adding FEn0 to FEnode'); end
  r1 = FEnode; r2=FEn0; 
 end

 [out,out1]=feutil(['add' Cam],r1,r2);
 if nargin<3 FEnode = out; out = 'done'; end

%Add element models - - - - - - - - - - - - - - - - - - -
elseif comstr(Cam,'feel')

i1 = max([size(FEelt,2) size(FEel0,2) size(FEel1,2) ]);
if ~isempty(FEelt)&size(FEelt,2)<i1 FEelt(1,i1)=0; end
if ~isempty(FEel0)&size(FEel0,2)<i1 FEel0(1,i1)=0; end
if ~isempty(FEel1)&size(FEel1,2)<i1 FEel1(1,i1)=0; end

i1=1; i3=[];while ~isempty(i1)
   if Cam(i1+4)=='t'      i2='FEelt';
   elseif Cam(i1+4)=='0'  i2='FEel0';
   elseif Cam(i1+4)=='1'  i2='FEel1'; end
   if i1==1 i3=[i2 '= [' i2]; else i3 = [i3 ';' i2]; end
   i1 = i1+4;
   i1 = i1+strfind(Cam(i1:length(Cam)),'feel')-1;
end
i3 = [i3 '];']; eval(i3);

%'AddSel' selection to main model - - - - - - - - - - - - - - - - - - -
elseif comstr(Cam,'sel')

   i1 = max(size(FEelt,2),size(FEel0,2));
   if isempty(FEelt) FEelt=FEel0;
   elseif size(FEel0,1)>1
    FEelt(size(FEelt,1)+[1:size(FEel0,1)],1:size(FEel0,2))=FEel0;
   else % check for single superelement xxx
    disp('AddSel: empty selection');
   end

   if DispFlag; femesh('infoFEelt'); end

% -----------------------------------------------------------------------
%'AddTest' 
% femesh('addtest (NodeIDShift)',node,ldraw);
elseif comstr(Cam,'test') [CAM,Cam]=comstr(Cam,5);

r1 = varargin{carg};carg=carg+1;
if ~isfield(r1,'Node') r1=struct('Node',r1);end
if ~isfield(r1,'Elt')  r1.Elt=varargin{carg};carg=carg+1;end

model=struct('Node',FEnode,'Elt',FEelt);

[model,i1]=feutil(['AddTest' CAM],model,r1);
FEnode=model.Node; FEelt=model.Elt;
if nargout>0 out=i1; end % return new nodes

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

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

% DivideInGroups  - - - - - - - - - - - - - - - - - - - - - - - -
if comstr(Cam,'in')

FEel0=feutil('divideingroups',FEnode,FEel0,[]);
 
% ('divide by face')  - - - - - - - - - - - - - - NOT DOCUMENTED
% divides by finding sharp edges in a group of quad4 elements
elseif comstr(Cam,'by')  [opt,CAM,Cam]=comstr(CAM,'by','%i');

[EGroup,nGroup]=getegroup(FEel0);
NNode(FEnode(:,1))=[1:length(FEnode(:,1))]'; elt=[];

for jGroup = 1:nGroup %loop on element groups
   [ElemF,i1,ElemP]= feutil('getelemf',FEel0(EGroup(jGroup),:),jGroup);
   if comstr(ElemP,'quad4')
    cEGI = EGroup(jGroup)+1:EGroup(jGroup+1)-1;
    elt(end+1,1:6) = [Inf abs('beam1')];
    i1=[1:4 1];i2=[];i3=[];i4=[];for j2 = 1:length(i1)-1
        i2=[i2 [1 length(NNode)+1]*sort(NNode(FEel0(cEGI,i1(j2+[0 1]))'))];
        i3=[i3 1:length(cEGI)];i4(j2*length(cEGI)+[-length(cEGI)+1:0])=j2;
    end
    r3=sparse(i3,i2,i4);
    
    i3=find(i2==1)'; elt(end+[1:length(i3)],1:2)= ...
         [rem(i3,length(NNode)+1) round(i3/length(NNode)+1)];

    for jElt=1:length(cEGI)
      [r2,r4]=basis(FEnode(NNode(FEel0(cEGI(jElt),1:4)),5:7));
      r1(jElt,1:3)=r2(:,3)';
    end
    for j2=find(any(r3))  % loop on edges
       i2=find(r3(:,j2));
       if length(i2)==1 | ( length(i2)>1 & abs(r1(i2(1),:)*r1(i2(2),:)')<.5)
         elt(end+1,1:2)=FEel0(cEGI(i2(1)),i1(r3(i2(1),j2)+[0 1]));
       end
    end
   end % type of element
end % of loop on groups
FEel0=elt;
 

% ('divide group (i) (element selectors)')  - - - - - - - - - - - - - -
elseif comstr(Cam,'group') 


 FEelt=feutil(['divide' CAM],FEnode,FEelt,FEel0,varargin{2:end});

% 'Divide div1 div2 div3 - - - - - - - - - - - - - - - - - - - - - -
else

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

dv1=linspace(0,1,opt(1))';dv2=linspace(0,1,opt(2))';dv3=linspace(0,1,opt(3))';

i1=0;while carg<=nargin&isa(varargin{carg},'double')
 r1=varargin{carg};carg=carg+1; i1=i1+1; if isempty(r1) r1=[0 1]; end
 if size(r1,1)>1
    r1=r1';
    if size(r1,2)>0 dv1=r1(1:max(find(r1(:,1))),1); else dv1=[0;1];end
    if size(r1,2)>1 dv2=r1(1:max(find(r1(:,2))),2); else dv2=[0;1];end
    if size(r1,2)>2 dv3=r1(1:max(find(r1(:,3))),3); else dv3=[0;1];end
    warning('Use multiple arguments for each direction of irregular divide');
    break
 else
    if i1==1 dv1=r1(:); elseif i1==2 dv2=r1(:); elseif i1==3 dv3=r1(:); end
 end
end

opt(1:3) = [length(dv1) length(dv2) length(dv3)];
[EGroup,nGroup]=getegroup(FEel0);
NNode(FEnode(:,1))=[1:length(FEnode(:,1))]';
elt=[];node=FEnode;

for jGroup = 1:nGroup %loop on element groups
   [ElemF,i1,ElemP]= feutil('getelemf',FEel0(EGroup(jGroup),:),jGroup);
   cEGI = EGroup(jGroup)+1:EGroup(jGroup+1)-1;
   iNode=fe_super('node',ElemP); % node indices

   if     strcmp(ElemP,'beam1')      
    if DispFlag;     disp('beam1 type division');    end
   elseif strcmp(ElemP,'quad4')
    if DispFlag;disp(sprintf('quad4 type division in %i by %i',opt(1:2)-1));end
   elseif strcmp(ElemP,'quadb')
    if DispFlag;disp(sprintf('quadb type division in %i by %i',opt(1:2)-1));end
   elseif strcmp(ElemP,'quad9')
    if DispFlag;disp(sprintf('quad9 type division in %i by %i',opt(1:2)-1));end
   elseif strcmp(ElemP,'hexa8')
    if DispFlag;disp(sprintf('hexa8 division in %i by %i by %i',opt(1:3)-1));end
   elseif strcmp(ElemP,'hexa20')
    if DispFlag;disp(sprintf('hexa20 division in %i by %i by %i',opt(1:3)-1));end
   elseif strcmp(ElemP,'tria3')
    if DispFlag;disp(sprintf('tria3 type division %i elt per edge',opt(1)));end
   end
   elt(size(elt,1)+1,1:size(FEel0,2)) = FEel0(EGroup(jGroup),:);

   for jElt=cEGI %loop on elements of current group

     if strcmp(ElemP,'beam1') % case for line divide
     % interpolate nodes
     i7 = FEel0(jElt,3:end); i2 = FEnode(NNode(FEel0(jElt,iNode)),:);
     i2 = ones(size(dv1))*i2(1,[5:7])+dv1*(i2(2,5:7)-i2(1,5:7));
     [node,i2] = femesh('AddNode',node,i2);i2=node(i2,1);

     %define new elements
     i4 = ones(opt(1)-1,1)*FEel0(jElt,:); 
     i4(:,iNode(1))=i2(1:opt(1)-1); i4(:,iNode(2))=i2(2:opt(1));

     % case for quad4 divide - - - - - - - - - - - - - - - - - - - -
     elseif strcmp(ElemP,'quad4')

     i7 = FEel0(jElt,5:end);  opt=opt-1;
     % i2 and i3 isoparametric coordinates of divided mesh
     if RunOpt.NewQuadDivide;
      i2 = dv1(:,ones(size(dv2)))*2-1; i3 = dv2(:,ones(size(dv1)))'*2-1;
     else
      i3 = dv1(:,ones(size(dv2)))*2-1; i2 = dv2(:,ones(size(dv1)))'*2-1;
     end
     xi = [-1 -1 0;1 -1 0;1 1 0;-1 1 0];
     n = (1+i2(:)*xi(:,1)').*(1+i3(:)*xi(:,2)')/4; % isoparametric shape fcn
     n = n*FEnode(NNode(FEel0(jElt,1:4)),5:7);
     [node,i2] = femesh('AddNode',node,n);i2=node(i2,1);
     % i4 node indices in i2
     i4 = [1:opt(1)];i4=i4(ones(opt(2),1),:)';
     i5=[0:opt(2)-1]*(opt(1)+1);i5=i5(ones(opt(1),1),:);
     i4=i4(:)+i5(:); i4 =[i4 i4+opt(1)+1 i4+opt(1)+2 i4+1];
     i4=reshape(i2(i4),size(i4,1),size(i4,2));	% Renumber
     i4=[i4 i7(ones(size(i4,1),1),:)];		% Add property numbers
     opt=opt+1;

     % case for quadb divide - - - - - - - - - - - - - - - - - - - -
     elseif strcmp(ElemP,'quadb')
     opt=opt-1;  i7 = FEel0(jElt,9:end);
     % basic cell
     xi = [-1 -1 0;1 -1 0;1 1 0;-1 1 0;0 -1 0;1 0 0;0 1 0;-1 0 0];
     coef=xi(:,1:2)+1;i1 = 2./opt;

     % w points in isoparametric coordinates
     % xi isoparametric coordinates of nodes
     % na shape functions are w points
     w=[coef(:,1)/opt(1) coef(:,2)/opt(2)]-1; % basic cell
     i2 = [1:size(w,1)]'; i2=i2(:,ones(opt(1),1)); i2=i2(:);
     i3 = [0:opt(1)-1]*i1(1);i3=i3(ones(size(w,1),1),:);i3=i3(:);
     w=w(i2,:);w(:,1)=w(:,1)+i3;				% repeat in x
     i2 = [1:size(w,1)]'; i2=i2(:,ones(opt(2),1)); i2=i2(:);
     i3 = [0:opt(2)-1]*i1(2);i3=i3(ones(size(w,1),1),:);i3=i3(:);
     w=w(i2,:);w(:,2)=w(:,2)+i3;			% repeat in y
     w(1,3)=0;
     i2=ones(1,size(xi,1));
     na = [(1+w(:,1)*xi(1:4,1)').*(1+w(:,2)*xi(1:4,2)')/4 ...
      (1-w(:,1).^2).*(1-w(:,2))/2      (1-w(:,2).^2).*(1-w(:,1))/2 ...
      (1-w(:,1).^2).*(1+w(:,2))/2      (1-w(:,2).^2).*(1+w(:,1))/2];
     na(:,1:4)=na(:,1:4)-na(:,5:8)*[1 0 0 1;1 1 0 0;0 1 1 0;0 0 1 1]'/2;
     i3=na*FEnode(NNode(FEel0(jElt,iNode)),5:7);
     [node,i2] = femesh('AddNode',node,i3);
     i2=node(i2,1);

     i4=[1:8]; i4=i4(ones(opt(1)*opt(2),1),:);
     i5=[0:opt(1)*opt(2)-1]'*8;i5=i5(:,ones(1,8));
     i4 = i4+i5;
     i4=reshape(i2(i4),size(i4,1),size(i4,2));	% Renumber
     i4=[i4 i7(ones(size(i4,1),1),:)];		% Add property numbers
     opt=opt+1;

     % case for quad9 divide - - - - - - - - - - - - - - - - - - - -
     elseif strcmp(ElemP,'quad9')
     opt=opt-1;  i7 = FEel0(jElt,9:end);
     % basic cell
     xi = [-1 -1 0;1 -1 0;1 1 0;-1 1 0;0 -1 0;1 0 0;0 1 0;-1 0 0;0 0 0];
     coef=xi(:,1:2)+1;i1 = 2./opt;

     % w points in isoparametric coordinates
     % xi isoparametric coordinates of nodes
     % na shape functions are w points
     w=[coef(:,1)/opt(1) coef(:,2)/opt(2)]-1; % basic cell
     i2 = [1:size(w,1)]'; i2=i2(:,ones(opt(1),1)); i2=i2(:);
     i3 = [0:opt(1)-1]*i1(1);i3=i3(ones(size(w,1),1),:);i3=i3(:);
     w=w(i2,:);w(:,1)=w(:,1)+i3;				% repeat in x
     i2 = [1:size(w,1)]'; i2=i2(:,ones(opt(2),1)); i2=i2(:);
     i3 = [0:opt(2)-1]*i1(2);i3=i3(ones(size(w,1),1),:);i3=i3(:);
     w=w(i2,:);w(:,2)=w(:,2)+i3;			% repeat in y
     w(1,3)=0;
     i2=ones(1,size(xi,1));
     na = [(1+w(:,1)*xi(1:4,1)').*(1+w(:,2)*xi(1:4,2)')/4 ...
      (1-w(:,1).^2).*(1-w(:,2))/2      (1-w(:,2).^2).*(1-w(:,1))/2 ...
      (1-w(:,1).^2).*(1+w(:,2))/2      (1-w(:,2).^2).*(1+w(:,1))/2];
     na(:,1:4)=na(:,1:4)-na(:,5:8)*[1 0 0 1;1 1 0 0;0 1 1 0;0 0 1 1]'/2;
     i3=na*FEnode(NNode(FEel0(jElt,iNode(1:8))),5:7);
     [node,i2] = femesh('AddNode',node,i3);i2=node(i2,1);

     i4=[1:9]; i4=i4(ones(opt(1)*opt(2),1),:);
     i5=[0:opt(1)*opt(2)-1]'*9;i5=i5(:,ones(1,9));
     i4 = i4+i5;
     i4=reshape(i2(i4),size(i4,1),size(i4,2));	% Renumber
     i4=[i4 i7(ones(size(i4,1),1),:)];		% Add property numbers
     opt=opt+1;

     % case for hexa8 divide - 8/8/00 - - - - - - - - - - - - - - - - - - -
     elseif strcmp(ElemP,'hexa8')

     i2 = FEnode(NNode(FEel0(jElt,iNode)),5:7);    i7 = FEel0(jElt,9:end);
     i1=dv1(:)*2-1;i1= i1(:,ones(1,opt(3)*opt(2))); i3=i1(:);
     i1= ones(opt(1),1)*(dv2(:)*2-1)';i1=i1(:)*ones(1,opt(3)); i3(:,2)=i1(:);
     i1= ones(opt(1)*opt(2),1)*(dv3(:)*2-1)'; i3(:,3)=i1(:);
     % this is an isoparametric mapping with linear elements
     coef =[0 0 0;-1 -1 -1;1 -1 -1;1 1 -1;-1 1 -1;-1 -1 1;1 -1 1;1 1 1;-1 1 1];
     i4=zeros(size(i3,1),3);
     for j1=1:8
	n=i3.*coef(ones(size(i3,1),1)*(j1+1),1:3)+1;
	i4 = i4 + (n(:,1).*n(:,2).*n(:,3)/8)*i2(j1,:);
     end
     [node,i2] = femesh('AddNode',node,i4);i2=node(i2,1);    
     % 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));

     %Renumber
     i4=reshape(i2(i4),size(i4,1),size(i4,2));i4=[i4 i7(ones(size(i4,1),1),:)];

     % case for hexa20 divide - 8/8/00 - - - - - - - - - - - - - - - - - - -
     elseif strcmp(ElemP,'hexa20')
     opt=opt-1;    i7 = FEel0(jElt,21:end);

     % basic cell on -1 1 interval
     coef=[-1 -1 -1;1 -1 -1;1 1 -1;-1 1 -1;-1 -1 1;1 -1 1;1 1 1;-1 1 1
         0 -1 -1;1 0 -1;0 1 -1;-1 0 -1;-1 -1 0;1 -1 0;1 1 0;-1 1 0; %9:16
         0 -1 1;1 0 1;0 1 1;-1 0 1];

     i4=(coef+1)/2; 
     i2=i4(:,1)*dv1(2:end)'+(1-i4(:,1))*dv1(1:end-1)';
     i4=[i2(:) reshape(i4(:,2*ones(1,opt(1))),20*opt(1),1) ...
               reshape(i4(:,3*ones(1,opt(1))),20*opt(1),1)] ; % x repeat
     i2=i4(:,2)*dv2(2:end)'+(1-i4(:,2))*dv2(1:end-1)';
     i4=[reshape(i4(:,ones(1,opt(2))),size(i4,1)*opt(2),1) i2(:) ...
               reshape(i4(:,3*ones(1,opt(2))),size(i4,1)*opt(2),1)] ; %y repeat
     i2=i4(:,3)*dv3(2:end)'+(1-i4(:,3))*dv3(1:end-1)';
     i4=[reshape(i4(:,ones(1,opt(3))),size(i4,1)*opt(3),1) ...
        reshape(i4(:,2*ones(1,opt(3))),size(i4,1)*opt(3),1) i2(:)] ; % z repeat
     i4=i4*2-1;

     %initial i4; mapped i3; isoparametric mapping with linear elements
     i2 = FEnode(NNode(FEel0(jElt,1:8)),5:7);   i3=zeros(size(i4));
     coef=[0 0 0;coef];
     for j1=1:8
	n=i4.*coef(ones(1,size(i4,1))*(j1+1),1:3)+1;
	i3 = i3 + (n(:,1).*n(:,2).*n(:,3)/8)*i2(j1,:);
     end
     [node,i2] = femesh('AddNode',node,i3);i2=node(i2,1);

     i4=[1:20]; i4=i4(ones(opt(1)*opt(2)*opt(3),1),:);
     i5=[0:opt(1)*opt(2)*opt(3)-1]'*20;i5=i5(:,ones(1,20));
     i4 = i4+i5;
     i4=reshape(i2(i4),size(i4,1),size(i4,2));	% Renumber
     i4=[i4 i7(ones(size(i4,1),1),:)];		% Add property numbers
     opt=opt+1;

     % case for tria3 divide - - - - - - - - - - - - - - - - - - - -
     elseif strcmp(ElemP,'tria3')

      r3 = ones(opt(1),1)*linspace(1,2,opt(1));
      r3=[nonzeros(triu(r3)) nonzeros(triu(r3'))];
      r3=[2-r3(:,1) r3(:,2)-1];r3(:,3)=1-sum(r3,2);

      [node,i2]=femesh('addnode',node,r3*FEnode(NNode(FEel0(jElt,1:3)),5:7)); 
      i2=node(i2,1); 
      i3=cumsum(1:opt(1));i3=i3(:); i4=[];
      while length(i3)>1
       i4=[i4;i3(1:end-1) i3(2:end) i3(2:end)-1 ...
          ones(length(i3)-1,1)*FEel0(jElt,4:end)];
       if length(i3)>2;
         i4=[i4;i3(2:end-1)-1 i3(2:end-1) i3(3:end)-1 ...
            ones(length(i3)-2,1)*FEel0(jElt,4:end)];
       end
       i3=i3(2:end)-1;
      end
      i4(:,1:3)=i2(i4(:,1:3));
      %i4=[i2([1 4 6;4 5 6;4 2 5;6 5 3]) ones(4,1)*FEel0(jElt,4:end)];

     else warning(sprintf('%s division not supported',ElemP));i4=[];
     end % of case selection - - - - - - - - - - - - - - - - - - - -
     elt(size(elt,1)+[1:size(i4,1)],1:size(i4,2)) = i4;

   end % of jElt
  end % of jGroup

  FEnode = node;	FEel0 = elt;   yPlot=2;

end % subcommand selection - - - - - - - - - - - - - - - - - - - - - -

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

% extrude nRep tx ty tz
opt = comstr(CAM,[-1 1 0 0 1]);
femesh(sprintf('Rev %i 0 0 0 0 0 %.20g %.20g %.20g',opt(1:4)),varargin{2:end});

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

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
% FindElt (selectors)
if comstr(Cam,'el')

if     comstr(Cam,'elt')
  [out,out1]=feutil(['find' CAM],FEnode,FEelt,FEel0,varargin{2:end});
elseif comstr(Cam,'el0')
  [out,out1]=feutil(['find' CAM],FEnode,FEel0,FEelt,varargin{2:end});
else error('Not a valid element description matrix');end


% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
% FindSym[Plane Orig# nx ny nz] symmetry condition
elseif comstr(Cam,'sym') [CAM,Cam]=comstr(CAM,4);
if comstr(Cam,'plane') [CAM,Cam]=comstr(CAM,6);

  if comstr(Cam,'o')
    opt=comstr(CAM(2:length(CAM)),[-1 0 0 0 0 0 0]);
    i1=opt(1:3); opt=opt(3:6);
  else
    opt =  comstr(CAM,[-1 0 0 0 0]);i1=find(FEnode(:,1)==opt(1));
    if isempty(i1) i1=[0 0 0]';else i1=FEnode(i1(1),[5:7])'; end
  end
  if norm(opt(2:4))==0 error('FindSymPlane: you must define nx ny nz');end
  i1(:,2)=opt(2:4)'/norm(opt(2:4));
  [EGroup,nGroup]=getegroup(FEelt);
  NNode(FEnode(:,1))=[1:length(FEnode(:,1))]';

  %i2 contains the Nodes of interest, i3 the DOFs of interest
  i2 = femesh(sprintf('FindNode Group %i:%i',1,nGroup'));
  i6=[eye(3) i1(:,1);0 0 0 1]* ...
     [eye(3)-2*i1(:,2)*i1(:,2)' zeros(3,1);zeros(1,3) 1]* ...
     [eye(3) -i1(:,1);0 0 0 1];
  
  i6 = (i6*[FEnode(NNode(i2),5:7) ones(length(i2),1)]')';
  [FEnode,i3]=femesh('AddNode',FEnode,i6(:,1:3));i4=[];i4(i3)=[1:length(i3)]';
  if any(~i4(i2)) error('Some nodes have no symmetric equivalent');end

  i5 = femesh(sprintf('FindDof Group %i:%i',1,nGroup'));
  out=[]; % contains the symmetry condition matrix

  for j1=1:length(i2)
     if i2(j1)==i3(j1)
     elseif i2(i4(i2(j1)))>i2(j1)
        out=[out;fe_c(i5,[i2(j1)])-fe_c(i5,i2(i4(i2(j1))))];
     end
  end
  out=out';
end % of plane

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
% FindDof determines DOFs in a set of element groups
elseif comstr(Cam,'dof') [CAM,Cam]=comstr(CAM,4);

 if carg<=nargin&isa(varargin{carg},'double')&~finite(varargin{carg}(1))
   elt=varargin{carg};carg=carg+1;
 else elt=FEelt; end

 if isempty(elt) error('FindDof : Element definition matrix is empty');end
 sdtw('femesh FindDof is obsolete use feutil(''GetDof'',model) instead');
 out=feutil(['getdof ' CAM],FEnode,elt,FEel0,varargin{carg:end});

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

if comstr(Cam,'get')
 [out1,out]=feutil(['findnode' CAM(4:end)],FEnode,FEelt,FEel0,varargin{2:end});
elseif nargout<2
  out=feutil(['findnode' CAM],FEnode,FEelt,FEel0,varargin{2:end});
else
  [out,out1]=feutil(['findnode' CAM],FEnode,FEelt,FEel0,varargin{2:end});
end

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

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

if comstr(Cam,'node') 

  [out1,out]=feutil(['find' CAM],FEnode,FEelt,FEel0,varargin{2:end});

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

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

help femesh

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


if comstr(Cam,'node')
  feutil(['info' Cam],FEnode,FEelt);
elseif comstr(Cam,'el0') feutil('infoeltFEel0',FEel0);
elseif comstr(Cam,'el') [opt,CAM,Cam]=comstr(CAM,'eltid','%i');

  if carg<=nargin opt=varargin{carg};carg=carg+1;
  else opt=[]; end
  if isempty(opt) feutil('infoeltFEelt',FEelt)
  else feutil('eltidinfo',FEelt,opt); end

else

if comstr(Cam,'feel') % info on a particular element definition matrix

  if     Cam(5)=='0' feutil('infoeltFEel0',FEel0);
  elseif Cam(5)=='1' feutil('infoeltFEel1',FEel1);
  elseif Cam(5)=='t' feutil('infoeltFEelt',FEelt);
  else feutil('infoeltFEelt',FEelt);feutil('infoeltFEel0',FEel0);end
else feutil('infoeltFEelt',FEelt);feutil('infoeltFEel0',FEel0);
end

if DispFlag
 disp(sprintf('\nFEnode contains %i nodes, FEn0 contains %i nodes', ...
     size(FEnode,1),size(FEn0,1)));
end

end

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

if comstr(Cam,'el0') isel0=1; elt=FEel0; [CAM,Cam]=comstr(CAM,4);
else elt=FEelt; isel0=0; end

elt=feutil(['join' CAM],elt);

if isempty(elt)
else
  if isel0 FEel0=elt; if DispFlag femesh('infoFEel0'); end
  else     FEelt=elt; if DispFlag femesh('infoFEelt'); end;end
end

% ------------------------------------------------------------------------
% return model structure
elseif comstr(Cam,'model')  [CAM,Cam] = comstr(CAM,6);

out=femesh;
if comstr(Cam,'0') out.Elt=FEel0;end
out=feutil('rmfield',out,'El0');

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

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
% 'ObjectHoleInPlate'
% CenterNode# Edge1Node# Edge2Node# r1 r2 nDiv1 nDiv2 nQuadrant 
if comstr(Cam,'holeinplate') [CAM,Cam]=comstr(CAM,12);

NNode(FEnode(:,1))=[1:length(FEnode(:,1))]';
if comstr(Cam,'o')
  opt = comstr(CAM(2:end),-1);
  node =[0 0 0 0 opt(1:3); FEnode(NNode(opt(4:5)),:)];
  opt=[0 opt(4:end)]; 

else  opt = comstr(CAM,-1); node = FEnode(NNode(opt(1:3)),:);
end
if length(opt)~=8 error('wrong number of arguments'); end

% i1 basis for the hole, i3 combinations of basis vectors
i1 = [node(2,5:7)-node(1,5:7);node(3,5:7)-node(1,5:7)];
i1 = [i1(1,:)/norm(i1(1,:))*opt(4);i1(2,:)/norm(i1(2,:))*opt(5);i1];

i3 = [ones(opt(6)+1,1)   linspace(0,1,opt(6)+1)';
      linspace(1-1/opt(7),-1,2*opt(7))'  ones(2*opt(7),1);
      -ones(2*opt(6),1) linspace(1-1/opt(6),-1,2*opt(6))';
      linspace(-1+1/opt(7),1,2*opt(7))' -ones(2*opt(7),1);
      ones(opt(6),1)   linspace(-1+1/opt(6),0,opt(6))'];

%i3 = [atan2(i3(:,2),i3(:,1)) i3];
%i3 = [cos(i3(:,1)) sin(i3(:,1)) i3(:,[2 3])];
% Modified for regular angular spacing by M. Corus

ang_i3_1=linspace(0,pi/4,opt(6)+1);
ang_i3_2=linspace(0,pi/4,opt(7)+1);

ang_i3=[ang_i3_1 ...
        ang_i3_2(2:end)+pi/4 ang_i3_2(2:end)+pi/2 ...
        ang_i3_1(2:end)+3*pi/4 ang_i3_1(2:end)+pi ...
        ang_i3_2(2:end)+5*pi/4 ang_i3_2(2:end)+3*pi/2 ...
        ang_i3_1(2:end)+7*pi/4]';
i3 = [cos(ang_i3) sin(ang_i3) i3];



if opt(8)>=4
elseif opt(8)>=3.5 i3 = i3(1:size(i3,1)-  opt(6),:);
elseif opt(8)>=3   i3 = i3(1:size(i3,1)-  opt(6)-  opt(7),:);
elseif opt(8)>=2.5 i3 = i3(1:size(i3,1)-  opt(6)-2*opt(7),:);
elseif opt(8)>=2   i3 = i3(1:size(i3,1)-2*opt(6)-2*opt(7),:);
elseif opt(8)>=1.5 i3 = i3(1:size(i3,1)-3*opt(6)-2*opt(7),:);
elseif opt(8)>=1   i3 = i3(1:size(i3,1)-3*opt(6)-3*opt(7),:);
elseif opt(8)>=0.5 i3 = i3(1:size(i3,1)-3*opt(6)-4*opt(7),:); end

i2 = [i3(:,1)*i1(1,:)+i3(:,2)*i1(2,:);i3(:,3)*i1(3,:)+i3(:,4)*i1(4,:)];
i2 = i2([1:size(i2,1)/2;size(i2,1)/2+1:size(i2,1)],:);
i2 = [[1:size(i2,1)]' zeros(size(i2,1),3) i2+node(ones(size(i2,1),1),5:7)];


[FEnode,i3]=femesh('AddNode',FEnode,i2);
i1 = size(i2,1);
if max(opt(4:5))<sqrt(sum((node(2,5:7)-node(1,5:7)).^2 ...
       +(node(3,5:7)-node(1,5:7)).^2))
    FEel0 = [2:2:i1-2; 1:2:i1-2; 3:2:i1 ; 4:2:i1]';
else
    FEel0 = [1:2:i1-2; 2:2:i1-2; 4:2:i1; 3:2:i1]';
end;


FEel0 = reshape(FEnode(i3(FEel0),1),size(FEel0,1),size(FEel0,2));

FEel0 = [Inf abs('quad4'); FEel0 ones(size(FEel0,1),1)*[1 1]];
yPlot=2;

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
% 'ObjectSphere'
% CenterNode# r iter
elseif comstr(Cam,'sphere') [CAM,Cam]=comstr(CAM,7);

error(1)


else % objects in feutil

 out=feutil(varargin{:});
 if isnumeric(out)&~finite(out(1)) FEel0=out;out='done';end

end % subcommand selection - - - - - - - - - - - - - - -

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

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
% OptimModel: removes unused nodes from FEnode
if comstr(Cam,'model')   [FEnode,FEelt]=feutil('optimmodel',FEnode,FEelt);

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
% OptimModel: removes unused nodes from FEnode
elseif comstr(Cam,'nodenum')
  [FEnode,FEelt]=feutil('optimnodenum',FEnode,FEelt);

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


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

if comstr(Cam,'el0')
  FEel0=feutil(['orient' CAM(4:end)],FEnode,FEel0,[],varargin{carg:end});
else
  FEelt=feutil(['orient' CAM],FEnode,FEelt,FEel0,varargin{carg:end});
end


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

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

 disp('Use the fecom(''nodetext'') command')

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
% PlotElt
elseif comstr(Cam,'el')

if comstr(Cam,'elt')
  st=sprintf('Plotting %s','FEelt');disp(st);  
  if comgui<105 feplot(FEnode,FEelt,[],[],2);
  else 
     if ~isfield(FEInfo,'Feplot')|ishandle(FEInfo.Feplot)
      cf = comgui('guife'); FEInfo.Feplot=cf.opt(1);
     else cf=feplot(FEInfo.Feplot); end
     cf.model=femesh; 
  end
elseif comstr(Cam,'el0')
  st=sprintf('Plotting %s','FEel0');disp(st);  
  if comgui<105 feplot(FEnode,FEel0,[],[],2);
  else
     if ~isfield(FEInfo,'Feplot')|ishandle(FEInfo.Feplot)
       cf = comgui('guife'); FEInfo.Feplot=cf.opt(1);
     else cf=feplot(FEInfo.Feplot); end
     cf.model={FEnode,FEel0};end
else error('Not a valid call to femesh PlotEl'); end

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

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

 r1=struct('mdl',femesh,'name','femesh');
 eval('feutilg(''initmodel'',r1);');

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

[FEnode,FEel0]=feutil('quad42quad9',FEnode,FEel0);

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

[FEnode,FEel0]=feutil('quad42q5p',FEnode,FEel0);

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

[FEnode,FEel0]=feutil('quad42quadb',FEnode,FEel0);

% ------------------------------------------------------------------------
% Tria32Tria6
elseif comstr(Cam,'tria32tria6') [CAM,Cam] = comstr(CAM,12);

  [EGroup,nGroup]=getegroup(FEel0);

  i2 =[]; % i2 first group nodes, i3 second group nodes
  i3=[1 2;2 3;3 1];
   for jGroup = 1:nGroup %loop on element groups
   ElemF= feutil('getelemf',FEel0(EGroup(jGroup),:),jGroup);
   if strcmp(ElemF,'tria3')
     FEel0(EGroup(jGroup),1:6)=[Inf abs('tria6')];
     for j1=EGroup(jGroup)+1:EGroup(jGroup+1)-1
       NNode(FEnode(:,1))=[1:length(FEnode(:,1))]';
       i2 = FEnode(NNode(FEel0(j1,1:3)),5:7);
       i2=(i2(i3(:,1),:)+i2(i3(:,2),:))/2;
       [FEnode,i2]=femesh('AddNode',FEnode,i2);
       FEel0(j1,4:8)=[FEnode(i2(:),1)' FEel0(j1,4:5)];
     end
   end
  end % of jGroup loop

% ------------------------------------------------------------------------
% hexa82hexa20
elseif comstr(Cam,'hexa82hexa20') [CAM,Cam] = comstr(CAM,13);

  [EGroup,nGroup]=getegroup(FEel0);

  i2 =[]; % i2 first group nodes, i3 second group nodes
  i3=[1 2;2 3;3 4;4 1;1 5;2 6;3 7;4 8;5 6;6 7;7 8;8 5];
  if size(FEel0,2)<11; FEel0(1,11)=0;end
  for jGroup = 1:nGroup %loop on element groups
   ElemF= feutil('getelemf',FEel0(EGroup(jGroup),:),jGroup);
   if strcmp(ElemF,'hexa8')
     FEel0(EGroup(jGroup),1:7)=[Inf abs('hexa20')];
     cEGI=EGroup(jGroup)+1:EGroup(jGroup+1)-1;
     r2=zeros(length(cEGI)*size(i3,1),3);
     NNode=sparse(FEnode(:,1),1,1:size(FEnode,1));
     for jElt=1:length(cEGI);
       r3 = FEnode(NNode(FEel0(cEGI(jElt),1:8)),5:7);
       r3=(r3(i3(:,1),:)+r3(i3(:,2),:))/2;
       r2(size(i3,1)*jElt+[-size(i3,1)+1:0],:)=r3;
     end
     [FEnode,i2]=femesh('AddNode',FEnode,r2);
     for jElt=1:length(cEGI);
       i3=i2(size(i3,1)*jElt+[-size(i3,1)+1:0]);
       FEel0(cEGI(jElt),9:23)=[FEnode(i3,1)' FEel0(cEGI(jElt),9:11)];
     end
   end
  end % of jGroup loop
% ------------------------------------------------------------------------
% hexa82hexa27
elseif comstr(Cam,'hexa82hexa27') [CAM,Cam] = comstr(CAM,13);

  [EGroup,nGroup]=getegroup(FEel0);

  i2 =[]; % i2 first group nodes, i3 second group nodes
  i3=[1 2;2 3;3 4;4 1;1 5;2 6;3 7;4 8;5 6;6 7;7 8;8 5];
  i4=[1 2 3 4;1 4 8 5;1 2 6 5;5 6 7 8;2 6 7 3;3 7 8 4];
  for jGroup = 1:nGroup %loop on element groups
   ElemF= feutil('getelemf',FEel0(EGroup(jGroup),:),jGroup);
   if strcmp(ElemF,'hexa8')
     if comstr(Cam,'b') FEel0(EGroup(jGroup),1:8)=[Inf abs('hexa27b')];
     else               FEel0(EGroup(jGroup),1:7)=[Inf abs('hexa27')];
     end
     for j1=EGroup(jGroup)+1:EGroup(jGroup+1)-1
       NNode(FEnode(:,1))=[1:length(FEnode(:,1))]';
       i2 = FEnode(NNode(FEel0(j1,1:8)),5:7);
       i2=[(i2(i3(:,1),:)+i2(i3(:,2),:))/2;
           (i2(i4(:,1),:)+i2(i4(:,2),:)+i2(i4(:,3),:)+i2(i4(:,4),:))/4;
           mean(i2(1:8,:))];
       [FEnode,i2]=femesh('AddNode',FEnode,i2);
       FEel0(j1,1:29)=[FEel0(j1,1:8) FEnode(i2(:),1)' FEel0(j1,9:10)];
     end
   end
  end % of jGroup loop
% ------------------------------------------------------------------------
% penta62penta15 
elseif comstr(Cam,'penta62penta15')  [CAM,Cam] = comstr(CAM,15);

  [EGroup,nGroup]=getegroup(FEel0);

  i2 =[]; % i2 first group nodes, i3 second group nodes
  i3=[1 2;2 3;3 1;1 4;2 5;3 6;4 5;5 6;6 4];

  for jGroup = 1:nGroup %loop on element groups
   ElemF= feutil('getelemf',FEel0(EGroup(jGroup),:),jGroup);
   if strcmp(ElemF,'penta6')
     FEel0(EGroup(jGroup),1:8)=[Inf abs('penta15')];
     for j1=EGroup(jGroup)+1:EGroup(jGroup+1)-1
       NNode(FEnode(:,1))=[1:length(FEnode(:,1))]';
       i2 = FEnode(NNode(FEel0(j1,1:6)),5:7);
       i2=(i2(i3(:,1),:)+i2(i3(:,2),:))/2;
       [FEnode,i2]=femesh('AddNode',FEnode,i2);
       FEel0(j1,1:17)=[FEel0(j1,1:6) FEnode(i2(:),1)' FEel0(j1,7:8)];
     end
   end
  end % of jGroup loop
 
% ------------------------------------------------------------------------
% tetra42tetra10 
elseif comstr(Cam,'tetra42tetra10')  [CAM,Cam] = comstr(CAM,15);

  [EGroup,nGroup]=getegroup(FEel0);

  i2 =[]; % i2 first group nodes, i3 second group nodes
  i3=[1 2;2 3;3 1;1 4;2 4;3 4];

  for jGroup = 1:nGroup %loop on element groups
   ElemF= feutil('getelemf',FEel0(EGroup(jGroup),:),jGroup);
   if strcmp(ElemF,'tetra4')
     FEel0(EGroup(jGroup),1:8)=[Inf abs('tetra10')];
     for j1=EGroup(jGroup)+1:EGroup(jGroup+1)-1
       NNode(FEnode(:,1))=[1:length(FEnode(:,1))]';
       i2 = FEnode(NNode(FEel0(j1,1:4)),5:7);
       i2=(i2(i3(:,1),:)+i2(i3(:,2),:))/2;
       [FEnode,i2]=femesh('AddNode',FEnode,i2);
       FEel0(j1,1:12)=[FEel0(j1,1:4) FEnode(i2(:),1)' FEel0(j1,5:6)];
     end
   end
  end % of jGroup loop
 

% ------------------------------------------------------------------------
% hexa2tetra
elseif comstr(Cam,'hexa2tetra')  [CAM,Cam] = comstr(CAM,7);
  [EGroup,nGroup]=getegroup(FEel0);

  i2 =[]; % i2 first group nodes, i3 second group nodes
  i4=[];

  for jGroup = 1:nGroup %loop on element groups
   ElemF= feutil('getelemf',FEel0(EGroup(jGroup),:),jGroup);
   if strcmp(ElemF,'hexa8')
     for j1=EGroup(jGroup)+1:EGroup(jGroup+1)-1
       NNode(FEnode(:,1))=[1:length(FEnode(:,1))]';
       i2 = FEnode(NNode(FEel0(j1,1:8)),5:7);
       i3=[1 2 3 4;5 6 7 8;1 2 6 5;2 3 7 6;3 4 8 7;1 4 8 5];
       i2=(i2(i3(:,1),:)+i2(i3(:,2),:)+i2(i3(:,3),:)+i2(i3(:,4),:))/4;%xxx
       i2(end+1,:)=sum(i2,1)/6;
       [FEnode,i2]=femesh('AddNode',FEnode,i2);
       i3=[FEel0(j1,1:8) i2'];
       i5=[i3(:,[1 2 end-6 end]);i3(:,[2 3 end-6 end]);
           i3(:,[3 4 end-6 end]);i3(:,[4 1 end-6 end]);
           i3(:,[6 5 end-5 end]);i3(:,[7 6 end-5 end]);
           i3(:,[8 7 end-5 end]);i3(:,[5 8 end-5 end]);
           i3(:,[2 1 end-4 end]);i3(:,[6 2 end-4 end]);
           i3(:,[5 6 end-4 end]);i3(:,[1 5 end-4 end]);
           i3(:,[3 2 end-3 end]);i3(:,[7 3 end-3 end]);
           i3(:,[6 7 end-3 end]);i3(:,[2 6 end-3 end]);
           i3(:,[4 3 end-2 end]);i3(:,[8 4 end-2 end]);
           i3(:,[7 8 end-2 end]);i3(:,[3 7  end-2 end]);
           i3(:,[1 4 end-1 end]);i3(:,[4 8 end-1 end]);
           i3(:,[8 5 end-1 end]);i3(:,[5 1 end-1 end])];
       i5(:,5:6)=ones(24,1)*FEel0(j1,9:10);
       i4=[i4;i5];
     end
     FEel0 = [FEel0(1:EGroup(jGroup)-1,:);
      Inf abs('tetra4') zeros(1,size(FEel0,2)-7)
      i4 zeros(size(i4,1),size(FEel0,2)-6)
      FEel0(EGroup(jGroup+1):size(FEel0,1),:)];%xxx
      EGroup=getegroup(FEel0);
   end
  end % of jGroup loop

if 1==2
  [EGroup,nGroup]=getegroup(FEel0);

  i2 =[]; % i2 first group nodes, i3 second group nodes
  for jGroup = 1:nGroup %loop on element groups
   ElemF= feutil('getelemf',FEel0(EGroup(jGroup),:),jGroup);
   if strcmp(ElemF,'hexa8')
     i3 = FEel0(EGroup(jGroup)+1:EGroup(jGroup+1)-1,1:10);
     i3 =reshape([i3(:,[1 2 3 6 9 10]);i3(:,[1 3 4 8 9 10]);
                  i3(:,[1 3 6 8 9 10]);i3(:,[6 8 5 1 9 10]);
                  i3(:,[6 8 7 3 9 10])]',6,5*size(i3,1))';
     FEel0 = [FEel0(1:EGroup(jGroup)-1,:);
      Inf abs('tetra4') zeros(1,size(FEel0,2)-7)
      i3 zeros(size(i3,1),size(FEel0,2)-6)
      FEel0(EGroup(jGroup+1):size(FEel0,1),:)];
      EGroup=getegroup(FEel0);
   end
  end % of jGroup loop
end
% ------------------------------------------------------------------------
% hexa2penta
elseif comstr(Cam,'hexa2penta');  [CAM,Cam] = comstr(CAM,7);
  [EGroup,nGroup]=getegroup(FEel0);

  i2 =[]; % i2 first group nodes, i3 second group nodes
  i4=[];

  for jGroup = 1:nGroup %loop on element groups
   ElemF= feutil('getelemf',FEel0(EGroup(jGroup),:),jGroup);
   if strcmp(ElemF,'hexa8')
     for j1=EGroup(jGroup)+1:EGroup(jGroup+1)-1
       NNode(FEnode(:,1))=[1:length(FEnode(:,1))]';
       i2 = FEnode(NNode(FEel0(j1,1:8)),5:7);
       i3=[1 2 3 4;5 6 7 8];
       i2=(i2(i3(:,1),:)+i2(i3(:,2),:)+i2(i3(:,3),:)+i2(i3(:,4),:))/4;
       [FEnode,i2]=femesh('AddNode',FEnode,i2);
       i3=[FEel0(j1,1:8) i2'];
       i5=[i3(:,[1 2 end-1 5 6 end]);
           i3(:,[2 3 end-1 6 7 end]);
           i3(:,[3 4 end-1 7 8 end]);
           i3(:,[4 1 end-1 8 5 end])];
       i5(:,7:8)=ones(4,1)*FEel0(j1,9:10);
       i4=[i4;i5];
     end

     FEel0 = [FEel0(1:EGroup(jGroup)-1,:);
      Inf abs('penta6') zeros(1,size(FEel0,2)-7)
      i4 zeros(size(i4,1),size(FEel0,2)-8)
      FEel0(EGroup(jGroup+1):size(FEel0,1),:)];
      EGroup=getegroup(FEel0);
   end
  end % of jGroup loop
  if ~any(Cam==';');
   sdtw('_nb','Be carefull using hexa2penta, always check mesh');
  end

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

  [EGroup,nGroup]=getegroup(FEel0);

  i2 =[]; % i2 first group nodes, i3 second group nodes
  for jGroup = 1:nGroup %loop on element groups
   ElemF= feutil('getelemf',FEel0(EGroup(jGroup),:),jGroup);
   if strcmp(ElemF,'quad4')
     i3 = FEel0(EGroup(jGroup)+1:EGroup(jGroup+1)-1,1:6);
     i3 = reshape([i3(:,[1 2 3 5 6]) i3(:,[1 3 4 5 6])]',5,2*size(i3,1))';
     FEel0 = [FEel0(1:EGroup(jGroup)-1,:);
      Inf abs('tria3') zeros(1,size(FEel0,2)-6)
      i3 zeros(size(i3,1),size(FEel0,2)-5)
      FEel0(EGroup(jGroup+1):size(FEel0,1),:)];
      EGroup=getegroup(FEel0);
   end
  end % of jGroup loop

% ------------------------------------------------------------------------
% REFINE functions
elseif comstr(Cam,'refine') [CAM,Cam] = comstr(CAM,7);

% ('refine beam MaxLength')  - - - - - - - - - - - - - -
if comstr(Cam,'beam')  opt=comstr(comstr(CAM,5),[-1 1]);

[EGroup,nGroup]=getegroup(FEel0);
NNode(FEnode(:,1))=[1:length(FEnode(:,1))]'; elt=[];

for jGroup = 1:nGroup %loop on element groups
   [ElemF,i1,ElemP]= feutil('getelemf',FEel0(EGroup(jGroup),:),jGroup);
   if comstr(ElemP,'beam1')
    elt=[elt;FEel0(EGroup(jGroup),:)];
    cEGI = EGroup(jGroup)+1:EGroup(jGroup+1)-1;
    for jElt=cEGI
     l = norm([1 -1]*FEnode(NNode(FEel0(jElt,1:2)),5:7));
     i1=fix((l-epsl)/opt)+2; if i1>1
      node=FEnode(NNode(FEel0(jElt,1))*ones(i1,1),5:7) ...
           +linspace(0,1,i1)'*diff(FEnode(NNode(FEel0(jElt,1:2)),5:7));
      [FEnode,i1]=femesh('addnode',FEnode,node);i1=FEnode(i1,1);
      elt=[elt;i1(1:end-1) i1(2:end) ones(size(i1,1)-1,1)*FEel0(jElt,3:end)];
     else elt=[elt;FEel0(jElt,:)];
     end
    end
   else elt=[elt;FEel0(EGroup(jGroup):EGroup(jGroup+1)-1,:)];
   end % type of element
end % of loop on groups
FEel0=elt;

else out='unknown'; end

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

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
% RemoveElt (findelt selectors)
if comstr(Cam,'elt') 
 FEelt = feutil(['remove' CAM],FEnode,FEelt,[],varargin{2:end});
elseif comstr(Cam,'el0') 
 FEel0 = feutil(['removeelt' CAM(4:end)],FEnode,FEel0,[],varargin{2:end});


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

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

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
% RepeatSel       nITE tx ty tz Orig# AnleDeg nx ny nz
if comstr(Cam,'sel')  

model=feutil(['repeat', CAM],FEnode,FEel0,[],varargin{carg:end});

FEnode=model.Node; FEel0=model.Elt; yPlot=2;


else out='unknown'; end % subcommand selection - - - - - - - - - - - - - -
% ------------------------------------------------------------------------
elseif comstr(Cam,'reset'); [CAM,Cam] = comstr(CAM,7);
 FEnode=[]; FEelt=[];FEel0=[];
 FE=struct('Node',[],'Elt',[],'pl',[],'il',[]);

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

% TransSel dx dy dz - - - - - - - - - - - - - - - - - - - - - - - - - - -
if comstr(Cam,'sel')|comstr(Cam,'el')

  if Cam(1)=='s' [CAM,Cam] = comstr(CAM,4);else [CAM,Cam] = comstr(CAM,3);end
  opt = comstr(CAM,[-1]); if isempty(opt)
   opt=varargin{carg}; carg=carg+1;
  end
  if length(opt)<3 opt(3)=0; end
  opt=opt(1:3); i6 = [eye(3) opt(:);0 0 0 1]; %translations
  femesh('repeatsel',i6);yPlot=2;

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
% translation/rotation of selected group of nodes
elseif comstr(Cam,'node') opt = comstr(CAM(5:length(CAM)),[-1 0 0 0]);

if DispFlag
 disp(['Translating FEn0 ' sprintf('tx %15.3f ty %15.3f tz %15.3f',opt)])
end
FEn0(:,5:7) = FEn0(:,5:7)+opt(ones(size(FEn0,1),1),1:3);

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

% ------------------------------------------------------------------------
% REVOLUTION functions REV
elseif comstr(Cam,'rev');  [CAM,Cam] = comstr(CAM,4);

if ~isempty(strfind(Cam,'-trans'))
 i1=strfind(Cam,'-trans');
 RunOpt=struct('type','trans');CAM(i1+[0:5])='';[CAM,Cam] = comstr(CAM,1);
elseif ~isempty(strfind(Cam,'-tria'))
 i1=strfind(Cam,'-tria');
 RunOpt=struct('type','tria');CAM(i1+[0:4])='';[CAM,Cam] = comstr(CAM,1);
else;RunOpt=struct('type','');
end
% Rev               nRep Orig# AngleDeg nx ny nz tx ty tz
if ~isempty(strfind(Cam,'o'))
   i1 = strfind(Cam,'o');
   opt = comstr(CAM([1:i1-1 i1+1:length(CAM)]), ...
             [-1 10 0 0 0  360    0  0  1  0  0  0]);
   i4 = [1 0 0 0 opt(2:4)];opt=[opt(1) 0 opt(5:length(opt))];
else
  opt = comstr(CAM,[-1 10    0     360    0  0  1  0  0  0]);
  if opt(1)==1&opt(3)==360 error('Rev: nRep=1 and Angle=360 is not valid');end
  i4=find(FEnode(:,1)==opt(2));
  if isempty(i4) i4 = [1  0 0 0  0 0 0]; else i4 = FEnode(i4,:); end
end


if carg<=nargin&isa(varargin{carg},'double')
  ind=varargin{carg};carg=carg+1;ind=ind(:)';
  i7 =ind;  opt(1)=length(ind)-1;
else ind=linspace(0,1,opt(1)+1); i7=0:opt(1); end

% revolution/extrusion of the groups in FEel0
[EGroup,nGroup]=getegroup(FEel0);
elt=[];for jGroup=1:nGroup
  [ElemF,opt1,ElemP]= feutil('getelemf',FEel0(EGroup(jGroup),:),jGroup);
  cEGI = EGroup(jGroup)+1:EGroup(jGroup+1)-1;

  if     strcmp(ElemP,'beam1') 
     iNode=[1 2];elt(size(elt,1)+1,1:6)=[Inf abs('quad4')];
     r1=FEel0(cEGI,[3:4]);
  % transition ring using tria6
  elseif     strcmp(ElemP,'beam3') &strcmp(RunOpt.type,'trans');
     iNode=[1 2 3];elt(size(elt,1)+1,1:6)=[Inf abs('tria6')];
     ind=sort([ind ind(1:end-1)+diff(ind)/2]);
     i7=[i7(:) i7(:)+[diff(i7(:))/2;0]]';i7=i7(:);i7=i7(1:end-1)';
     opt(6:9)=opt(6:9)/2;opt(1)=opt(1)*2;
     if size(FEel0,2)<10; FEel0(1,10)=0;end
     r1=FEel0(cEGI,[4:5]);
  elseif     strcmp(ElemP,'beam3') &strcmp(RunOpt.type,'tria');
     iNode=[1 2 3];elt(size(elt,1)+1,1:6)=[Inf abs('tria6')];
     ind=sort([ind ind(1:end-1)+diff(ind)/2]);
     i7=[i7(:) i7(:)+[diff(i7(:))/2;0]]';i7=i7(:);i7=i7(1:end-1)';
     opt(6:9)=opt(6:9)/2;opt(1)=opt(1)*2;
     if size(FEel0,2)<10; FEel0(1,10)=0;end
     r1=FEel0(cEGI,[4:5]);
  elseif     strcmp(ElemP,'beam3') 
     iNode=[1 2 3];elt(size(elt,1)+1,1:6)=[Inf abs('quadb')];
     ind=sort([ind ind(1:end-1)+diff(ind)/2]);
     i7=[i7(:) i7(:)+[diff(i7(:))/2;0]]';i7=i7(:);i7=i7(1:end-1)';
     opt(6:9)=opt(6:9)/2;opt(1)=opt(1)*2;
     if size(FEel0,2)<10; FEel0(1,10)=0;end
     r1=FEel0(cEGI,[4:5]);
  elseif     strcmp(ElemP,'mass1') 
     iNode=[1];elt(size(elt,1)+1,1:6)=[Inf abs('beam1')];
     r1=ones(length(cEGI),1)*[1 1];
  elseif strcmp(ElemP,'tria3') 
     iNode=[1:3];elt(size(elt,1)+1,1:7)=[Inf abs('penta6')];
     r1=FEel0(cEGI,4:5);
  elseif strcmp(ElemP,'tria6') 

     iNode=[1:6];elt(size(elt,1)+1,1:8)=[Inf abs('penta15')];
     ind=sort([ind ind(1:end-1)+diff(ind)/2]);
     i7=[i7(:) i7(:)+[diff(i7(:))/2;0]]';i7=i7(:);i7=i7(1:end-1)';
     opt(6:9)=opt(6:9)/2;opt(1)=opt(1)*2;
     if size(FEel0,2)<17; FEel0(1,17)=0;end
     r1=FEel0(cEGI,9:10);

  elseif strcmp(ElemP,'quad4')
     iNode=[1:4];elt(size(elt,1)+1,1:6)=[Inf abs('hexa8')];
     r1=FEel0(cEGI,5:6);
  elseif strcmp(ElemP,'quadb')
     iNode=[1:8];elt(size(elt,1)+1,1:7)=[Inf abs('hexa20')];
     ind=sort([ind ind(1:end-1)+diff(ind)/2]);
     i7=[i7(:) i7(:)+[diff(i7(:))/2;0]]';i7=i7(:);i7=i7(1:end-1)';
     opt(6:9)=opt(6:9)/2;opt(1)=opt(1)*2;
     if size(FEel0,2)<10; FEel0(1,10)=0;end
     r1=FEel0(cEGI,9:10);
  elseif strcmp(ElemP,'quad9') iNode=[1:4];
     sdtw('_nb','quad9 extruded into hexa8');
     iNode=[1:4];elt(size(elt,1)+1,1:6)=[Inf abs('hexa8')];
     r1=FEel0(cEGI,10:11);
  else                     error([ElemP ' not supported by Rev']); end

  NNode(FEnode(:,1))=[1:length(FEnode(:,1))]';

  % i1 angles for revolution meridians i3 revolution basis i4 origin
  i1 = -ind*opt(3)*pi/180;

  % i2 initial nodes
  i2=NNode(FEel0(cEGI,iNode)');i2 = FEnode(i2,5:7);
  i6 = [1:size(i2,1)]'*ones(1,length(ind));
  i7 = i7(ones(1,size(i2,1)),:);  i1=i1(ones(1,size(i2,1)),:);

  if opt(3)~=0 %rotate the nodes
    i2 = i2 - i4(ones(length(iNode)*length(cEGI),1),5:7);
    
    i3 = basis(opt(4:6),[1 0 0],1);    i3 = i3(:,[2 3 1]);
    i2 = i2*i3;
    i2(:,1:2) = [sqrt(i2(:,1).^2+i2(:,2).^2) atan2(i2(:,1),i2(:,2))];

    % i2 rotated nodes
    i2 = i4(ones(size(i2,1)*(opt(1)+1),1),5:7) +  ...  % origin
         i2(i6,3)*i3(:,3)' + ... %vertical pos
         (i2(i6,1)*ones(1,3)) .* ...
                           (sin(i1(:)+i2(i6,2)) *i3(:,1)' + ...
                            cos(i1(:)+i2(i6,2)) *i3(:,2)');
  else 
    i2=i2(i6,:);
  end

  if size(r1,1)==length(cEGI) % propagate properties
   r1=[r1(:,ones(1,length(ind))) r1(:,2*ones(1,length(ind)))];
   r1=reshape(r1,size(r1,1)*length(ind),2);
  end


  i2=i2+i7(:)*opt(1,7:9); %translate nodes

  [FEnode,i2]=femesh('AddNode',FEnode,i2);
  i2=reshape(FEnode(i2,1),length(iNode),length(i2)/length(iNode))';
  i3=1:size(i2,1)-length(cEGI);  i6=size(elt,1);

    if strcmp(ElemP,'beam1')
      i5=find(i2(i3,1)-i2(i3,2));
      if length(i5)~=length(i3)
       warning('Removing degenerate beam from extrusion');
       i3=i3(i5); % xxx deal with MatId propagation
      end
      elt(i6+[1:length(i3)],1:6)= ...
              [i2(i3,1:2) i2(i3+length(cEGI),[2 1]) r1(i3,:)];
    elseif strcmp(ElemP,'beam3') &strcmp(RunOpt.type,'trans');
     i5=1:4:size(i2,1);
     for j1=i5(1:end-1)-1
      elt(end+[1:3],1:6)=[
      i2(j1+1,[1 2]) i2(j1+3,1) i2(j1+1,3) i2(j1+2,[3 1]) 
      i2(j1+1,2) i2(j1+5,2) i2(j1+3,1) i2(j1+3,2)  i2(j1+4,3)  i2(j1+2,3)   
     i2(j1+3,1) i2(j1+5,[2 1]) i2(j1+4,3) i2(j1+5,3) i2(j1+4,1)     ];
    end

    elseif strcmp(ElemP,'beam3') &strcmp(RunOpt.type,'tria');

      i3=[length(cEGI):2*length(cEGI):size(i2,1)-1]';i5=[-length(cEGI)+1:0];
      i3=i3(:,ones(size(i5)))+i5(ones(size(i3)),:);i3=i3(:);i5=length(cEGI);

      elt(i6+[1:2*length(i3)],1:8)= ...
       [i2(i3,1:2) i2(i3+2*i5,1)    i2(i3,3) i2(i3+i5,3)  ...
        i2(i3+i5,1) r1(i3,:) ;
        i2(i3,2) i2(i3+2*i5,[2 1])  i2(i3+i5,2) i2(i3+2*i5,3) ...
        i2(i3+i5,3) r1(i3,:)];

     opt(1)=opt(1)/2;opt(6:9)=opt(6:9)*2;

    elseif strcmp(ElemP,'beam3')
      i3=[length(cEGI):2*length(cEGI):size(i2,1)-1]';i5=[-length(cEGI)+1:0];
      i3=i3(:,ones(size(i5)))+i5(ones(size(i3)),:);i3=i3(:);i5=length(cEGI);

      elt(i6+[1:length(i3)],1:10)= [i2(i3,1:2) i2(i3+2*i5,[2 1]) ...
                              i2(i3,3) i2(i3+i5,2) i2(i3+2*i5,3) i2(i3+i5,1) ...
                               r1(i3,:)];

     opt(1)=opt(1)/2;opt(6:9)=opt(6:9)*2;

    elseif strcmp(ElemP,'tria6')

      i3=[length(cEGI):2*length(cEGI):size(i2,1)-1]';i5=[-length(cEGI)+1:0];
      i3=i3(:,ones(size(i5)))+i5(ones(size(i3)),:);i3=i3(:);i5=length(cEGI);
      elt(i6+[1:length(i3)],1:17)= [i2(i3,1:3) i2(i3+2*i5,1:3) ...
                                    i2(i3,4:6) i2(i3+i5,1:3) ...
                                    i2(i3+2*i5,4:6)  r1(i3,:)];
      opt(1)=opt(1)/2;opt(6:9)=opt(6:9)*2;

    elseif strcmp(ElemP,'mass1')
      elt(i6+[1:length(i3)],1:4)= ...
              [i2(i3,1) i2(i3+length(cEGI),1) r1(i3,:)];
    elseif strcmp(ElemP,'tria3')
      elt(i6+[1:length(i3)],1:8)= ...
              [i2(i3,:) i2(i3+length(cEGI),:) r1(i3,:)];

    elseif strcmp(ElemP,'quad4')|strcmp(ElemP,'quad9')
      elt(i6+[1:length(i3)],1:10)= ...
              [i2(i3,:) i2(i3+length(cEGI),:) r1(i3,:)];

    elseif strcmp(ElemP,'quadb')
      i3=[length(cEGI):2*length(cEGI):size(i2,1)-1]';i5=[-length(cEGI)+1:0];
      i3=i3(:,ones(size(i5)))+i5(ones(size(i3)),:);i3=i3(:);i5=length(cEGI);
      elt(i6+[1:length(i3)],1:22)= ...
       [i2(i3,1:4) i2(i3+2*i5,1:4) ... % corner nodes
              i2(i3,5:8) i2(i3+i5,1:4) i2(i3+2*i5,5:8) r1(i3,:)];
      %              r1(1:length(i3),:)]; was modified everywhere
     opt(1)=opt(1)/2;opt(6:9)=opt(6:9)*2;
    end
end % of loop on jGroup

FEel0=feutil('orient',FEnode,elt);
%yPlot=2;

% ------------------------------------------------------------------------
% RotateSel Orig# AnleDeg nx ny nz
elseif comstr(Cam,'rotatesel') [CAM,Cam] = comstr(CAM,10);

if ~isempty(strfind(Cam,'o'))
    i1=strfind(Cam,'o');
    opt=comstr(CAM([1:i1-1 i1+1:length(CAM)]),[-1 0 0 0  0   1 0 0]);
    i5=opt(1:3); opt=[0 opt(4:7)];
else
    opt=comstr(CAM,[-1 0  0  1  0  0]);
    i1=find(FEnode(:,1)==opt(1));
    if isempty(i1) 
     if DispFlag;disp('femesh Rotate : origin set to 0,0,0'); end
     i5=[0 0 0];
    else i5=FEnode(i1(1),5:7);end
end
if opt(2)==0 error('Need a non-zero angle for Rotate');end
opt(2)=opt(2)/180*pi;
if norm(opt(3:5))==0 error('Need a non-zero direction for Rotate');end
i5(3,:)=opt(3:5)/norm(opt(3:5));
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];

% i5 [origin;tx ty tz;rx ry rz];
i6=[eye(3) i5(1,:)';0 0 0 1]* ...
   [cos(opt(2))*eye(3)+i5(3,:)'*i5(3,:)*(1-cos(opt(2)))+ ...
	    sin(opt(2))*i5(4:6,1:3) zeros(3,1);zeros(1,3) 1]* ...
   [eye(3) -i5(1,:)';0 0 0 1];
%i6 = i6*[eye(3) i5(2,:)'*j1;0 0 0 1]; %translations

femesh('repeatsel',i6); %yPlot=2;

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

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
% 'selelt'
if comstr(Cam,'el')

 [i4,FEel0]=femesh(['Find' CAM],varargin{2:end});

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

 if comstr(Cam,'a')
   [EGroup,nGroup]=getegroup(FEel0);
   [CAM,Cam] = comstr(CAM,2);
   i1 = comstr(CAM,-1); 
   if isempty(i1)&nargin>=carg i1=varargin{carg};carg=carg+1; end
   if isempty(i1) error('No group selection specified'); end
   if any(i1<0|i1>nGroup) error('Invalid group selection'); end
   i2=[];for j1 = i1; i2 = [i2 EGroup(j1):EGroup(j1+1)-1]; end
   FEel0 = FEel0(i2,:);
   st=['FEel0 contains FEel0 group(s)' sprintf('  %i',i1)]; 
   if DispFlag; disp(st); end
 else
   [EGroup,nGroup]=getegroup(FEelt);
   i1 = comstr(CAM,-1);
   if isempty(i1)&nargin>=carg i1=varargin{carg};carg=carg+1; end
   if isempty(i1) error('No group selection specified'); end
   if any(i1<0|i1>nGroup) error('Invalid group selection'); end
   i2=[];for j1 = i1; i2 = [i2 EGroup(j1):EGroup(j1+1)-1]; end
   FEel0 = FEelt(i2,:);
   st=['FEel0 contains FEelt group(s)' sprintf('  %i',i1)]; 
   if DispFlag; disp(st); end
 end

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
% SelNode finds the node closest to the current point
elseif comstr(Cam,'node')

   NNode(FEnode(:,1))=[1:length(FEnode(:,1))]';
   if ~isempty(Cam) % selection of nodes given by number 
     i1 = comstr(CAM(5:length(CAM)),[-1]); % indices in node numbers
     i1 = i1(find(i1<=length(NNode)&i1>0));
     i1 = NNode(i1);	  % indices in position number
     i1 = i1(find(i1));
   else error('not available');
   end
   if ~isempty(i1)
    FEn0 = FEnode(i1,:);
     st=['Selected nodes (stored in FEn0) ' 10  ...
         sprintf('%6i%6i%6i%6i%6i%6i%6i%6i\n',FEnode(i1,1))];
   else st='No node selected'; end
   if DispFlag; disp(st); end


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

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

% 'SetGroup [i] Mat[j] Prop [k] Name [s] EGID (i) - - - - - - - - - - - - - -
if comstr(Cam,'group') [CAM,Cam] = comstr(CAM,6);

if comstr(Cam,'a')&~comstr(Cam,'all') elt=FEel0;st1='FEel0'; [CAM,Cam] = comstr(CAM,2);
else elt = FEelt;st1='FEelt'; end

[elt,st]=feutil(['setgroup' CAM],FEnode,elt,[],st1);
if DispFlag ; disp(st); end;

if comstr(st1,'FEel0') FEel0=elt; else FEelt=elt; end

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

% ------------------------------------------------------------------------
% stick pushes the nodes of FEn1 onto a surface of quad selected in FEel0
elseif comstr(Cam,'stick')

[EGroup,nGroup]=getegroup(FEel0);
[i1,r1] = feutil(sprintf('findnode group1:%i',nGroup),FEnode,FEel0,[]);
NNode=[];NNode(i1)=1:length(i1);
r1=r1(:,5:7);
if carg<=nargin&isa(varargin{carg},'double')
  node=varargin{carg};carg=carg+1;
else node=FEn1; end
ind = find(finite(FEel0(:,1)));

for j1=1:size(node,1)
  r2 = ((r1-node(j1*ones(size(r1,1),1),5:7)).^2)*[1;1;1];
  [r3,i2]=min(r2);
  [i3,i4] = find(FEel0(ind,1:4)==i1(i2));
  [r3,i4]=min(sum(r2(NNode(FEel0(ind(i3),1:4)))'));i3=ind(i3(i4));
  r3 = r1(NNode(FEel0(i3,1:4)),:);
  r3=r3-ones(size(r3,1),1)*mean(r3);
  r4 = basis(r3(2,:)/norm(r3(2,:))+r3(3,:)/norm(r3(3,:)),r3(3,:));
  out1(j1,1:3)=r4(:,3)'*((node(j1,5:7)-r1(i2,:))*r4(:,3));
  node(j1,5:7) = node(j1,5:7)-r4(:,3)'*((node(j1,5:7)-r1(i2,:))*r4(:,3));
end
if nargout==0 FEn1=node; else out=node; end

% ------------------------------------------------------------------------
% SymSel Orig# nx ny nz
elseif comstr(Cam,'symsel') [CAM,Cam] = comstr(CAM,7);

if comstr(Cam,'o')
    opt=comstr(CAM(2:length(CAM)),[-1 0 0 0 0 0 0]);
    i5=opt(1:3); opt=opt(3:6);
else
    opt =  comstr(CAM,[-1 0 0 0 0]);i1=find(FEnode(:,1)==opt(1));
    if isempty(i1) i5=[0 0 0];else i5=FEnode(i1(1),[5:7]); end
end

if norm(opt(2:4))==0 error('Need a non-zero direction for Sym');end
i5(3,:)=opt(2:4)/norm(opt(2:4));
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];

% i5 [origin;0 0 0;tx ty tz;rx ry rz];
i6=[eye(3) i5(1,:)';0 0 0 1]* ...
   [eye(3)-2*i5(3,:)'*i5(3,:) zeros(3,1);zeros(1,3) 1]* ...
   [eye(3) -i5(1,:)';0 0 0 1];
%i6 = i6*[eye(3) i5(2,:)'*j1;0 0 0 1]; %translations

femesh('repeatsel',i6);%yPlot=2;

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

out.Node=[1 0 0 0   0.  10. 0.
          2 0 0 0   0.  0.  0.
          3 0 0 0   10. 0.  0.
          4 0 0 0   10. 10. 0.
          5 0 0 0   2.  2.  0.
          6 0 0 0   8.  3.  0.
          7 0 0 0   8.  6.  0.
          8 0 0 0   4.  7.  0.   ];

out.Elt=[Inf abs('quad4');
       1 2 5 8 1 1
       2 3 6 5 1 1
       3 4 7 6 1 1
       4 1 8 7 1 1
       5 6 7 8 1 1 ];

% ------------------------------------------------------------------------
% Test a series of test objects
elseif comstr(Cam,'test') [CAM,Cam] = comstr(CAM,5); 

RunOpt=struct('Struct',0,'Div',[],'Back',0,'Plot',0);
i1=strfind(Cam,'struct');
if isempty(i1); RunOpt.Struct=0;
else RunOpt.Struct=1;CAM(i1+[0:5])='';[CAM,Cam]=comstr(CAM,1);
end

i1=strfind(Cam,'-back');
if ~isempty(i1);RunOpt.Back=1;CAM(i1+[0:4])='';[CAM,Cam]=comstr(CAM,1);end
i1=strfind(Cam,'back');
if ~isempty(i1);RunOpt.Back=1;CAM(i1+[0:3])='';[CAM,Cam]=comstr(CAM,1);end
i1=strfind(Cam,'plot');
if ~isempty(i1);RunOpt.Plot=1;CAM(i1+[0:3])='';[CAM,Cam]=comstr(CAM,1);end

% 100 steel
% 101 water
% 102 air
FE.pl=m_elastic('dbval 100 Steel','dbval 101 Water','dbval 102 Air');
r1=FE.pl(find(FE.pl(:,1)==100),:);r1(1)=112;r1([3 5])=2*r1([3 5]);
FE.pl=m_elastic(FE.pl,'dbval 112',r1);

% 110 Shell 1cm
% 111 full solid
% 112 1x1 beam section
% 113 full 2D solid
FE.il=[110 fe_mat('p_shell',1,1) 0 0 0 .01 0]; % default shell thickness
FE.il= p_solid(FE.il,'dbval 113 Plane stress',[111 fe_mat('p_solid','SI',1) 0 0 0]);
FE.il= p_beam(FE.il,'dbval 112 rectangle .02 .02');

st={'2bay','ubeam','1dbar'};ElemF='';
for j1=1:length(st); if comstr(Cam,st{j1}) ElemF=st{j1}; ElemP=st{j1}; end; end

if isempty(ElemF) % match a standard element test
 st={'q4pb',[100 113],'q8pb',[100 113],'t3pb',[100 113],'t6pb',[100 113],...
     'q4p',[100 113],'q8p',[100 113],'t3p',[100 113], ...
     't6p',[100 113],'q5p',[100 113],'q9a',[100 113],...
     'hexa8b',[100 111],'hexa20b',[100 111],'hexa27b',[100 111], ...
     'penta6b',[100 111],'penta15b',[100 111],'tetra4b',[100 111], ...
     'tetra10b',[100 111],...
     'hexa8',[100 111],'hexa20',[100 111],'penta6',[100 111], ...
     'penta15',[100 111],'tetra4',[100 111],'tetra10',[100 111],...
     'tria3',[100 110],'tria6',[100 110],'quad4',[100 110], ...
     'quadb',[100 110],'quadc',[100 110], ...
     'quad9',[100 110],'mitc4',[100 110], ...
     'dktp',[100 110], ...
     'bar1',[100 112],'beam1',[100 112], ...
     'flui4',[102 111],'flui6',[102 111],'flui8',[102 111],...
     'celas',[1 1],'cbush',[1 1],'mck3',[1 1]};
 st=reshape(st,2,length(st)/2)';
 ElemF=''; i5=[]; j1=1;
 for j1=1:size(st,1); 
   if comstr(Cam,st{j1,1}) ElemF=st{j1,1};i5=st{j1,2}; break; end; 
 end
 if isempty(ElemF) error('Not a femesh teststruct'); end
 ElemP=feval(ElemF,'parent'); FEel0=[Inf abs(ElemF)];
end

%
div=strfind(Cam,'divide'); % Allow for divide definition in command
if ~isempty(div);
 i1=div;[div,i2,i3,i4]=sscanf(Cam(i1+6:end),'%i');
 CAM(i1:i1+4+i4)='';[CAM,Cam]=comstr(CAM,1);
 if length(div)<3;div(end+1:3)=3;end
end

if ~isempty(div);RunOpt.Div=1;
elseif ~isempty(strmatch(ElemP,strvcat('quad4','quadb','tria3','tria6')))
 if RunOpt.Struct; div=[6 6];else; div=[1 1];end % division default for 2D
elseif RunOpt.Struct; div=[3 3 3];
else div=[1 1 1];
end

% - - - - - - - - - - - - - - - - - - - - - - - create elements
% - - - - - - - - - - - - - - -  - - - - - - - - - - 2D 
if comstr(ElemP,'quad4') | comstr(ElemP,'quadb')

  FEel0(1,1:6)=[Inf abs('quad4')];
  FEnode=[1 0 0 0  0 0 0;2 0 0 0  0 1 0;3 0 0 0 1 1 0;4 0 0 0 1 0 0];
  FEel0(end+1,1:6)=[1 4 3 2 i5];
  femesh(sprintf('divide%i %i',div(1:2)),';');
  if comstr(ElemP,'quadb')  femesh('quad42quadb',';');  end
  femesh(sprintf('set group a1 name %s matid %i proid %i',ElemF,i5));
elseif comstr(ElemP,'tria3') | comstr(ElemP,'tria6')
  FEnode=[1 0 0 0  0 0 0;2 0 0 0  0 1 0;3 0 0 0 1 1 0;4 0 0 0 1 0 0];
  if RunOpt.Struct;  
   femesh('teststructquad4 divide 1 1'); femesh('quad2tria');
   femesh(sprintf('divide%i %i',div(1:2)),';');
  else      FEel0=[Inf abs('tria3') ;1 2 3 i5 0 ];
  end
  if comstr(ElemP,'tria6')  femesh('tria32tria6'); end
  femesh(sprintf('set group a1 name %s matid %i proid %i',ElemF,i5));
elseif comstr(ElemP,'q5p')
  FEnode=[1 0 0 0  0 0 0;2 0 0 0  0 1 0;3 0 0 0 1 1 0;
          4 0 0 0  1 0 0;5 0 0 0 .5 .5 0];
  FEel0(end+1,1:7)=[1 4 3 2 5 i5];
elseif comstr(ElemP,'quad9')
  FEnode=[1 0 0 0  0 0 0;4 0 0 0  0 1 0;3 0 0 0 1 1 0;2 0 0 0 1 0 0;
  5 0 0 0  .5 .0 .0; 6 0 0 0 1. .5 .0;7 0 0 0 .5 1. 0.;8 0 0 0 .0 .5 0
  9 0 0 0 .5 .5 0];
  FEel0(end+1,1:11)=[1:9 i5];
% - - - - - - - - - - - - - - -  - - - - - - - - - - 1D 
elseif comstr(ElemP,'beam1')
  FEnode=[1 0 0 0  0 0 0;2 0 0 0  1 0 0;3 0 0 0  0 0 1];
  FEel0(end+1,1:5)=[1:2 i5 3];
  femesh(sprintf('divide%i',div(1)),';');
elseif comstr(ElemP,'mck3')
  FEnode=[1 0 0 0  0 0 0;2 0 0 0  1 0 0;3 0 0 0  0.5 0 0];
  FEel0(end+1,1:10)=[1:3 i5 3 10 10 1000 1000 .2];
  femesh(sprintf('divide%i',div(1)),';');
elseif comstr(ElemP,'celas')
  FEnode=[1 0 0 0  0 0 0;2 0 0 0  1 0 0];
  FEel0(end+1,1:4)=[1:2 i5 3];
  femesh(sprintf('divide%i',div(1)),';');
elseif comstr(ElemP,'cbush')
  FEnode=[1 0 0 0  0 0 0;2 0 0 0  1 0 0];
  FEel0(end+1,1:4)=[1:2 i5 3];
  femesh(sprintf('divide%i',div(1)),';');
% - - - - - - - - - - - - - - -  - - - - - - - - - - 3D 
elseif comstr(ElemP,'tetra4')
  if ~RunOpt.Struct
   FEnode=[1 0 0 0  0 0 0;2 0 0 0  1 0 0;3 0 0 0 1 1 0;4 0 0 0 0 0 1];
   FEel0(end+1,1:6)=[1 2 3 4 i5];
  else
   if ~isempty(div);femesh(sprintf('testhexa8 divide %i %i %i',div(1:3)));
   else;femesh('teststructhexa8');
   end
   femesh hexa2tetra; 
  end
  femesh(strcat('set group a1 name ',ElemF,sprintf(' mat %i pro %i',i5(1),i5(2))));
elseif comstr(ElemP,'tetra10')
  if ~RunOpt.Struct; femesh('test tetra4');femesh('tetra42tetra10');
  else
   femesh('testhexa8 divide 1 1 1');
   femesh(sprintf('divide%i %i %i',div(1:3)),';');
   femesh hexa2tetra;  femesh('tetra42tetra10');
   FEel0(1,1:length(ElemF)+1)=[Inf abs(ElemF)];
  end
elseif comstr(ElemP,'penta6')
  if ~RunOpt.Struct
    FEnode=[1 0 0 0  0 0 0;2 0 0 0  1 0 0;3 0 0 0  1 1 0;
            4 0 0 0  0 0 1;5 0 0 0  1 0 1;6 0 0 0  1 1 1];
    FEel0(end+1,1:8)=[1:6 i5];
  else
   femesh('testhexa8 divide 1 1 1');
   femesh(sprintf('divide%i %i %i',div(1:3)),';');
   femesh('hexa2penta;');
  end
  femesh(strcat('set group a1 name ',ElemF,sprintf(' mat %i pro %i',i5(1),i5(2))));
elseif comstr(ElemP,'penta15')
  if ~RunOpt.Struct;   femesh('test penta6');femesh('penta62penta15');
  else
   femesh('testhexa8 divide 1 1 1');
   femesh(sprintf('divide%i %i %i',div(1:3)),';');
   femesh('hexa2penta;');femesh penta62penta15; 
  end
  femesh(strcat('set group a1 name ',ElemF,sprintf(' mat %i pro %i',i5(1),i5(2))));
elseif comstr(ElemP,'hexa8')
  FEnode=[1 0 0 0  0 0 0;2 0 0 0  1 0 0;3 0 0 0  1 1 0;4 0 0 0  0 1 0;
          5 0 0 0  0 0 1; 6 0 0 0  1 0 1;7 0 0 0  1 1 1;8 0 0 0  0 1 1];
  FEel0(end+1,1:10)=[1:8 i5];
  if RunOpt.Struct|RunOpt.Div; 
   femesh(sprintf('divide%i %i %i',div(1:3)),';');
  end
  femesh(strcat('set group a1 name ',ElemF,sprintf(' mat %i pro %i',i5(1),i5(2))));
elseif comstr(ElemP,'hexa20')
  femesh('testhexa8 divide 1 1 1');
  femesh(sprintf('divide%i %i %i',div(1:3)),';');
  femesh('hexa82hexa20');
  femesh(strcat('set group a1 name ',ElemF,sprintf(' mat %i pro %i',i5(1),i5(2))));
elseif comstr(ElemP,'hexa27b')
  if ~RunOpt.Struct; femesh('test hexa8');femesh('hexa82hexa27b');
  else
   femesh('testhexa8 divide 1 1 1');
   femesh(sprintf('divide%i %i %i',div(1:3)),';');
   femesh('hexa82hexa27b');
   femesh(strcat('set group a1 name ',ElemF,sprintf(' mat %i pro %i',i5(1),i5(2))));
  end
% - - - - - - - - - - - - - - -  - - - - - - - - - - OTHERS
elseif comstr(ElemP,'2bay')
  FEnode=[1 0 0 0  0 0 0;2 0 0 0    0 1 0;
          3 0 0 0  1 0 0;4 0 0 0    1 1 0;
          5 0 0 0  2 0 0;6 0 0 0    2 1 0];
  FEel0=[Inf abs('beam1'); 1 3 1 1 0 0;2 4 1 1 0 0;3 4 1 2 0 0;1 4 1 3 0 0;
         Inf abs('beam1'); 3 5 1 1 0 0;4 6 1 1 0 0;5 6 1 2 0 0;3 6 1 3 0 0];
  FEelt=FEel0; 
  FE.pl= m_elastic(FE.pl,'dbval 1 Steel');
  if nargout==1
    FE.il = [1 1 5.0000e-09 5.0000e-09   5.0000e-09   2.0000e-05   % longerons
     2 1     5.0000e-09   5.0000e-09   5.0000e-09   2.0000e-05	  % diagonals
     3 1     5.0000e-09   5.0000e-09   5.0000e-09   2.0000e-05]; % battens
   out=femesh; out.DOF=feutil('getdof',out);
   if RunOpt.Plot; feplot(out);end; return;
  end

elseif comstr(ElemP,'ubeam')

 FEnode=[1 0 0 0  -.5 -.5 0;2  0 0 0  -.5+1/6 -.5 0;3 0 0 0  -.5 .5-1/6 0
         4 0 0 0  -.5+1/6 .5-1/6 0;5 0 0 0  -.5 .5 0;6 0 0 0 -.5+1/6 .5 0
         7 0 0 0 .5-1/6 .5 0;8 0 0 0 .5 .5 0;9 0 0 0 .5-1/6 .5-1/6 0
         10 0 0 0 .5 .5-1/6 0;11 0 0 0 .5-1/6 -.5 0;12 0 0 0 .5 -.5 0];
 FEelt=[Inf abs('quad4');4 6 5 3 1 1;9 10 8 7 1 1];
 FEel0=[Inf abs('quad4');1 2 4 3 1 1;11 12 10 9 1 1];
 femesh(';divide 5 1;addsel;');
 FEel0=[Inf abs('quad4');4 6 7 9 1 1];
 femesh(';divide 4 1;addsel;');
 femesh('join group 1:3');
 femesh(';selgroup1;extrude 10 0 0 .25;');
 FEelt=FEel0; 
 FEelt(2:end,1:8)=FEelt(2:end,[5:8 1:4]);
 FEelt=feutil('orient;',FEnode,FEelt);FEel0=FEelt;
 FE.pl=m_elastic('dbval 1 Steel');
 FE.il=p_solid('dbval 1 full');

 if nargout==1;  out=femesh;  if RunOpt.Plot;feplot(out);end;return; end

elseif comstr(ElemP,'bar')  % 'test bar'
error('Change name xxx');
  FEnode=[1 0 0 0    0   0  0;2 0 0 0    50  0  0];
  FEel0=[Inf abs('bar1')];
  FEel0(2,:)=[1 2 1 1 1];
  femesh(sprintf('divide%i',div(1)),';');
  FEelt=FEel0; %FEel0=[];
  FE.pl=m_elastic('dbval 1 Steel');
  FE.il= p_beam('dbval 1 rectangle 1 1');

else warning('No element created')
end

% - - - - - - - - - - - - - - - - - - - - - - - - - end of create elements

if RunOpt.Back;
  out=struct('Node',FEnode,'Elt',FEel0,'pl',FE.pl,'il',FE.il);out1=[];
  return;
end

if nargout==1 & ischar(out)
 out=struct('Node',FEnode,'Elt',FEel0,'pl',FE.pl,'il',FE.il);
end

Rand=strfind(Cam,'rand'); % Allow internal node randomization
if ~isempty(Rand);
 i1=Rand;[Rand,i2,i3,i4]=sscanf(Cam(i1+4:end),'%g');
 CAM(i1:i1+2+i4)='';[CAM,Cam]=comstr(CAM,1);
else; Rand=0;end

if RunOpt.Struct %  - - - - - - - - - - - - - - - - STRUCT TESTS

model=struct('Node',FEnode,'Elt',FEel0,'pl',FE.pl,'il',FE.il);

switch Cam(1:3)
case {'bea','bar'} % - - - - - - - - - - - - - - - - 1D

  model=femesh('test2bay'); model.DOF=[];
  if ~exist('CAM'); Cam='eig'; end
  if comstr(Cam,'bar1')
   model.Elt=feutil('set group 1:2 name bar1',model);
   model=fe_case(model,'fixdof','2D',.03);
  end
  % reassign new material to half the structure
  model.Elt(feutil('findelt withnode {x>1}',model),3)=112;

  model=fe_case(model,'fixdof','Edge','x==0');
  [m,k,mdof]=fe_mk(model);
  [Case,model.DOF]=fe_mknl('init',model);
  k1=fe_mknl('assemble',model,Case,1);  
  if ~isequal(mdof,Case.DOF) i1=fe_c(Case.DOF,mdof,'ind');k1=k1(i1,i1);end  
  if norm(diag(k)-diag(k1))/norm(diag(k))>sqrt(eps) 
   error('Problem in NL assembly');
  end

  if ~isempty(strfind(Cam,'eig'))
    def=fe_eig(model,[105 6 1e3 11]);
  else % load
    data=struct('sel','groupall','dir',[0 1 0]);
    model=fe_case(model,'FVol','Gravity',data);
    def = fe_load(model);
  end

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

  if Rand; % Randomize node positions
   i1=find(FEnode(:,5)<1);
   FEnode(i1,5)=FEnode(i1,5).*(1+rand(size(i1))/div(1)*Rand);
   i1=find(FEnode(:,6)<1);
   FEnode(i1,6)=FEnode(i1,6).*(1+rand(size(i1))/div(2)*Rand);
  end
  model.Node=FEnode;model.Elt=feutil('orient;',model);
 
  if ~isempty(strfind(Cam,'_0'));     model.il(find(model.il(:,1)==113),3)=0; 
  elseif ~isempty(strfind(Cam,'_1')); model.il(find(model.il(:,1)==113),3)=1; 
  elseif ~isempty(strfind(Cam,'_2')); model.il(find(model.il(:,1)==113),3)=2; 
  end

  if nargin==2 model.pl=varargin{2}; end % deals with specified pl (eg. anisotropic)

  % attribute material 112 to 1/2 structure
  i2=feval(ElemF,'prop'); ind=feutil('findelt withnode {x>.5}',model);
  model.Elt(ind,i2(1))=112;

  model=fe_case(model,'fixdof','Edge','x==0');

  if comstr(Cam(3),'p') % 2D
   data=struct('sel','groupall','dir',[1 0 0]);
   model=fe_case(model,'FVol','Gravity',data);
   data=struct('sel','y==1','eltsel','withnode {x>0.51}', ...
               'def',1,'DOF',.02);
   model=fe_case(model,'Fsurf','Surface load',data);
   data=struct('sel','y==1','eltsel','withnode {x>0.51}', ...
               'def',1,'DOF',.19);
   model=fe_case(model,'Fsurf','Pressure load',data);
  else % SHELL
   data=struct('sel','groupall','dir',[0 0 9.81]);
   model=fe_case(model,'FVol','Gravity',data);
  end
 
  if ~isempty(strfind(Cam,'eig'))  def=fe_eig(model,[105 10 1e3 11]);
  else                             def = fe_load(model);
  end

case {'tet','pen','hex'} % - - - - - - - - - - - - - - - - - - - - - 3D 

  model.Elt=feutil('orient;',model);

  if Rand; % Randomize node positions
   i1=find(FEnode(:,5)<1);
   FEnode(i1,5)=FEnode(i1,5).*(1+rand(size(i1))/div(1)*Rand);
   i1=find(FEnode(:,6)<1);
   FEnode(i1,6)=FEnode(i1,6).*(1+rand(size(i1))/div(2)*Rand);
   i1=find(FEnode(:,7)<1);
   FEnode(i1,7)=FEnode(i1,7).*(1+rand(size(i1))/div(3)*Rand);
  end
  FEnode(:,7)=FEnode(:,7)/3/max(FEnode(:,7));model.Node=FEnode;

  % attribute material 112 to 1/2 structure
  ind=feutil('findelt innode {x>.5}',model);
  i2=feval(ElemF,'prop');
  model.Elt(ind,i2(1))=112;
  model=fe_case(model,'fixdof','Base','x==0');

  if ~isempty(strfind(Cam,'eig')) def=fe_eig(model,[105 10 1e3 11]);
  else % load
    data=struct('sel','groupall','dir',[0 0 9.81]);
    model=fe_case(model,'FVol','Gravity',data);
    def = fe_load(model);
    data=struct('sel','z==1/3','eltsel','withnode {x>0.51}', ...
               'def',1,'DOF',.01);
    model=fe_case(model,'Fsurf','Surface load',data);
    data=struct('sel','z==1/3','eltsel','withnode {x>0.51}', ...
                'def',1,'DOF',.19);
    model=fe_case(model,'Fsurf','Pressure load',data);
  end
 
case 'flu'   % - - - - - - - - - - - - - - - - - - - fluids

   model=femesh('model0');
   FEnode(:,7)=FEnode(:,7)/3/max(FEnode(:,7));model.Node=FEnode;

   [m,k,mdof]=fe_mk(model);
   [Case,model.DOF]=fe_mknl('init',model);
   k1=fe_mknl('assemble',model,Case,1);  
   if ~isequal(mdof,Case.DOF) i1=fe_c(Case.DOF,mdof,'ind');k1=k1(i1,i1);end  
   if norm(diag(k)-diag(k1))/norm(diag(k))>sqrt(eps) 
    error('Problem in NL assembly');
   end
   if isempty(strfind(Cam,'eig')) & isempty(strfind(Cam,'load'))
     Cam='eig';
   end
   if ~isempty(strfind(Cam,'eig')) % eig
    model = fe_case(model,'AddToCase 1','FixDof','clamped dofs','z==0');
    def=fe_eig(model,[105 10 1e3 11]);
   end


otherwise % - - - - - - - - - - - - - - - - - - - Others
end


 % the same for all elements
 if ~RunOpt.Back; feplot(model,def);end
 if ~isempty(strfind(Cam,'str')) % try display stress
  if isempty(strfind(Cam,'eig'))
   [m,k,mdof]=fe_mknl(model);def.def=ofact(k,def.def);
   feplot('initdef',def);
  end
  fecom(';colordatastress'); 
 elseif ~isempty(strfind(Cam,'ener')) % try display stress
  if isempty(strfind(Cam,'eig'))
   [m,k,mdof]=fe_mknl(model);def.def=ofact(k,def.def);
   feplot('initdef',def);
  end
  fecom(';colordataenerk'); 
 end
 out=model;out1=def; % return def

else % Not a basic elt test

 %if ~isempty(div); femesh(strcat('divide',sprintf(' %i',div)));end
 if isstruct(out); 
  out.Node=FEnode;out.Elt=FEel0;
  out.Elt=feutil('orient;',out);
 end
end % basic test of elements

if RunOpt.Plot
 if isfield(out,'Elt') mdl=out; 
 else mdl=femesh; if isempty(FEelt) mdl.Elt=FEel0; end
 end
 feplot(mdl);  
 if sp_util('issdt') 
  cf=feplot;assignin('base','cf',cf);out=cf.mdl;
  disp('Base workspace variable ''cf'' points to the feplot figure')
 end
end

% ------------------------------------------------------------------------
% UnJoin command disconnects nodes that are common to two groups
elseif comstr(Cam,'unjoin')
 [CAM,Cam] = comstr(CAM,7);  opt = fix(comstr(CAM,-1));
 if isempty(opt)
  opt=varargin{carg};carg=carg+1;opt1=varargin{carg};carg=carg+1;
 else opt1=opt(2); opt=opt(1); end

  NNode(FEnode(:,1))=[1:length(FEnode(:,1))]';
  [EGroup,nGroup]=getegroup(FEelt);

 % i2 first group nodes, i3 second group nodes
 if ischar(opt)
   i2=feutil(sprintf('findnode inelt {%s}',opt),FEnode,FEelt);
 else
   i2=feutil(['findnode group' sprintf('%i ',opt)],FEnode,FEelt);
 end
 if ~ischar(opt1) opt1=sprintf('group %i',opt1);end
 [ind,elt]=feutil(sprintf('findelt %s',opt1),FEnode,FEelt);
 i3=feutil('findnode groupall',FEnode,elt);
  
 [i4,i5]=intersect(i3,i2); % Nodes to duplicate
 [FEnode,i8] = femesh('AddNodeNew',FEnode,FEnode(NNode(i4),:));
 out=[i4 FEnode(i8,1)];
 NNode=sparse(FEnode(:,1),1,FEnode(:,1));NNode(i4)=FEnode(i8,1); % New nodes
 
 for jGroup=1:nGroup
    ElemF= feutil('getelemf',FEelt(EGroup(jGroup),:),jGroup);
    i7=fe_super('node',ElemF);
    cEGI=EGroup(jGroup)+1:EGroup(jGroup+1)-1;
    i6=intersect(cEGI,ind);
    if ~isempty(i6)
      FEelt(i6,i7)=reshape(full(NNode(FEelt(i6,i7))),length(i6),length(i7));
    end
 end
 disp('Interface separated FEnode and FEelt modified');

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

if nargout==0&ischar(out)&comstr(out,'done') clear out;
elseif ischar(out)&comstr(out,'unk')
 sdtw('_nb','%s is not a known femesh command',CAM0);clear out
end

if isfield(FE,'Node') % new format with persistent variable
     FE.Node=FEnode;FE.Elt=FEelt;
     if ~isempty(FEel0) FE.El0=FEel0; end
     if ~isempty(FEn0) FE.N0=FEn0; end
     if ~sdtdef('isdeployed'); r1=dbstack;if length(r1)>1& ...
          ~isempty(strfind(r1(1).name,'femesh'))
        assignin('caller','FEnode',FEnode);       
        assignin('caller','FEelt',FEelt);       
        assignin('caller','FEel0',FEel0);       
     end;end

end
