function [o1,o2,o3]=fe_load(varargin)

%FE_LOAD Construction of loads (input shape matrices)
%
%	Syntax : load = fe_load(model,Case)
%
%       This function evaluates the CASE (see doc('sdt/case')) using the
%       defined MODEL (see doc('sdt/fem'))
%       
%       CASE is a structure with fields CASE.Stack and Case.DOF
%
%	Supported loads are DofLoad, DofSet, FVol and FSurf see more details
%         in doc('fe_load')
%
%       See also help fe_case, fe_c
%                doc  fe_load, sdt/case, femk, fe_case

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

if comstr(varargin{1},'cvs')
 o1='$Revision: 1.63 $  $Date: 2006/06/09 09:56:48 $'; return;
end

epsl=sdtdef('epsl');
if nargin==0; help('fe_load'); return; end
% string commands - - - - - - - - - - - - - - - - - - - - - - - - - - - -
if ischar(varargin{1}) 
 [CAM,Cam]=comstr(varargin{1},1);carg=2;
% building the load - - - - - - - - - - - - - - - - - - - - - - - - - - - -
if comstr(Cam,'buildu'); [CAM,Cam]=comstr(CAM,7);

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

 o1=fe_load(model,Case);
 o2=[];
 if comstr(Cam,'time');
   
  if carg<=nargin; t=varargin{carg}; carg=carg+1; end
  if carg<=nargin; opt=varargin{carg}; carg=carg+1; end

  if ~isstruct(opt); error(1); end

  if isstruct(o1) & isfield(o1,'def'); opt.nc=size(o1.def,2);
  else;                                opt.nc=size(o1,2); 
  end

  ft=[]; % each column for each fc.def
  if isstruct(o1) & isfield(o1,'curve') &  ~isempty(o1.curve) 
   ft=o1.curve;
   if isempty(ft); error('Curve not found'); end
   ft=fe_curve('returny -extrapby0',ft,t,model);
  elseif isfield(opt,'TimeVector');  ft=zeros(opt.nt,opt.nc);ft(1:opt.nt,:)=1;
  elseif length(opt.Opt)>=6 ;
         ft=zeros(opt.nt,opt.nc);ft(1:opt.Opt(6),:)=1;  % Step of fixed length
  else;  ft=ones(opt.nt,opt.nc);
  end
  o2=ft;
 end % comstr(Cam,'time')
else
 sdtw('_nb','Not a fe_load command');
end

return;
else % string commands or traditionnal model,Case
 carg=1;Cam=''; CAM='';
 model=varargin{carg};carg=carg+1;
end

%-----------------------------------------------------------------------------------
%-----------------------------------------------------------------------------------
if ~isstruct(model) % - - - - - - - - - - - - - - - - - - - - - - - - -
  error('you must define a MODEL and a CASE to call FE_LOAD');
else
 NoT=0; Case=[]; % Give load on active DOFs
 while carg<=nargin
  r1=varargin{carg}; carg=carg+1;
  if ischar(r1)&  strcmpi(r1,'not') NoT=1;
  elseif ischar(r1) 
    Case=fe_case(model,['getcase' r1]);
  elseif isfield(r1,'Stack')
    Case=r1;
  end
 end % Loop on arguments
 if isempty(Case) Case=fe_case(model,'getcase');end
 % really only case was given
 if isempty(Case.Stack)&~isfield(model,'Elt') 
         Case=model;model=struct('DOF',Case.DOF);
 elseif isa(Case,'cell')&size(Case,2)==3; Case=struct('Stack',{Case});
 elseif isstruct(Case)&isfield(Case,'Stack')&size(Case.Stack,2)==3
   if ~NoT; Case=fe_case(model,'gett',Case);end
 else; error('Not a valid case specification');
 end
end % ~isstruct

% building the load - - - - - - - - - - - - - - - - - - - - - - - - - - - -
if isfield(model,'Node')
 NNode(model.Node(:,1))=[1:length(model.Node(:,1))]';
 node=model.Node;
else node=[1];
end
if ~isfield(model,'DOF')|isempty(model.DOF); 
  model.DOF=feutil('getdof',model);
end

mdof=model.DOF; i1=find(mdof>0);
%nd=sparse(round(rem(mdof(i1),1)*100),fix(mdof(i1)),1:length(i1),100, ...
%    max([fix(max(mdof(i1)))+1 max(node(:,1))]));
N=length(mdof);
%if max(max(nd))>N error('Something wrong with MDOF'); end
b=[]; bset=struct('DOF',[],'def',[]);bset.lab={};
lab={}; curve={}; bi=[];dofi=[]; r1=[];
pl=[]; ID=[];

for j1=1:size(Case.Stack,1) % loop on loads

r1=Case.Stack{j1,3};

switch comstr(Case.Stack{j1,1},-27)
case 'fvol'

 if ischar(r1.sel) [i1,elt]=feutil(['findelt' r1.sel],model);
 elseif ~finite(r1.sel(1)) elt=r1.sel; 
 else error('Data.Sel must contain elements or a selection'); end
 if isempty(elt); error('Element selection for volume load is empty');end
 b1=volume_load(model,Case,elt,r1,NNode);

case 'fsurf'

 if isempty(pl) pl=fe_mat('getpl',model); il=fe_mat('getil',model);end
 b1=surface_load(model,pl,il,r1,NNode);

case 'dofload' 

 [r2,r3,b1]=fe_c(model.DOF,r1.DOF);
 if size(r1.def,1)==length(r3); b1=b1'*r1.def;
 elseif isempty(r2)
  sdtw('_nb','Load ''%s'' affects no DOF',Case.Stack{j1,2});
 else % actually some DOFs are missing
  b1=b1'*r1.def(fe_c(r1.DOF,r2,'ind'),:);
 end

case 'dofset' 
  ind=find(rem(round(r1.DOF*100),100)==0); % deal with wild card spec
  if ~isempty(ind);
   i1=fe_c(model.DOF,r1.DOF(ind),'ind');
   r1.DOF(ind)=[];r1.def(ind,:)=[];r1.DOF=[bset.DOF;model.DOF(i1)];
   r1.def(end+[1:length(i1)],1:length(i1))=eye(length(i1));
  end
 if nargout==2 % return a DOFsets as a second ouput

  if isempty(bset.DOF); bset.DOF=r1.DOF;
  else; bset.DOF=unique(round([bset.DOF;r1.DOF]*100))/100;
  end
  ind=fe_c(bset.DOF,r1.DOF,'ind');
  bset.def(ind,end+[1:size(r1.def,2)])=r1.def;
  if size(r1.def,2)==1;
   st1=fe_curve('datatypecell','displacement');
   bset.lab(end+1,1:3)={Case.Stack{j1,2} st1{2:3}};
  elseif size(r1.def,2)>1
   st1=fe_curve('datatypecell','displacement');
   for j2=1:size(r1.def,2) 
    bset.lab(end+1,1:3)={sprintf('%s-u%i',Case.Stack{j1,2},j2) st1{2:3}}; 
   end
  end
  if isfield(r1,'curve');bset.curve=r1.curve;end
  if isfield(r1,'ID');bset.ID=r1.ID;end
  b1=[];

 elseif ~isfield(model,'K')|isempty(model.K)| ...
   size(model.K{1},1)~=size(model.DOF,1)
    sdtw('_nb',['%s ''%s'' entry ignored because model.K was' ...
      ' not defined'], Case.Stack{j1,1:2}); b1=[];
 elseif isfield(Case,'TIn') & isfield(Case,'TIn')
  b1=spalloc(length(model.DOF),0,1); for j2=1:length(model.K)
   try; r2=model.K{j2}*Case.TIn; 
   catch;eval('r2=feutilb(''a*b'',model.K{j2},Case.TIn);')
   end
   b1(:,end+[1:size(r2,2)])=r2;
  end
 else
  ind=fe_c(model.DOF,r1.DOF,'ind');cind=1:length(model.DOF);cind(ind)=0;
  cind=find(cind);
  b1=zeros(length(model.DOF),0); for j2=1:length(model.K)
   r2=model.K{j2}(cind,ind)*r1.def; b1(cind,end+[1:size(r2,2)])=r2;
  end
 end

case {'fixdof','keepdof','par','rigid','sensdof','mpc','nastran', ...
   'rbe3','sensstrain'}
 b1=[];
otherwise
  b1=[]; fprintf('%s entry not supported by fe_load',Case.Stack{j1,1})
end

 if     isfield(r1,'Curve'); st1=r1.Curve;
 elseif isfield(r1,'curve'); st1=r1.curve;
 else st1='';
 end

 % labels are ported from load definition
 if isfield(r1,'lab')&size(r1.lab,1)==size(b1,2); 
   lab(size(b,2)+[1:size(r1.lab,1)],1:size(r1.lab,2))=r1.lab;
   ind=size(b,2); for j2=1:size(b1,2) 
    if iscell(st1); curve{ind+j2}=st1{j2};else; curve{ind+j2}=st1;end
   end
 elseif size(b1,2)==1; 
  ind=size(b,2); lab{ind+1,1}=Case.Stack{j1,2};  
  if iscell(st1); curve(ind+[1:length(st1)])=st1;
   else; curve{ind+1}=st1;end
 elseif size(b1,2)>1
  ind=size(b,2); for j2=1:size(b1,2) 
   lab{ind+j2,1}=sprintf('%s-u%i',Case.Stack{j1,2},j2); 
   if iscell(st1); curve{ind+j2}=st1{j2};else; curve{ind+j2}=st1;end
  end
 end

 %if isfield(Case,'T')&~isempty(Case.T)&size(b1,1)==size(Case.T,1)&~NoT 
 %  b1=Case.T'*b1;
 %end % done below and should thus not be done twice (see clean up below)
 if isfield(r1,'ID'); % Fix the ID if needed
  if length(r1.ID)<size(b1,2); 
    r1.ID(end+1:size(b1,2))=max([r1.ID(:);max(ID)])+ ...
     [1:size(b1,2)-length(r1.ID)];
  end
  r1.ID=r1.ID(:); ID(end+[1:size(b1,2)],1)=r1.ID(1:size(b1,2));
 end

 if isempty(b) b=b1; else b=[b b1];end
 b1=[];

end % loop on stack

o1=struct('DOF',model.DOF,'def',b,'lab',{lab},'fun',[0 13], ...
  'ID',ID,'curve',{curve});
% fun(2)=13 is for UFF DataType excitation force

o2=bset;
if isfield(Case,'T')&size(Case.T,1)==size(o1.def,1)&~NoT 
  o1.def=Case.T'*o1.def; o1.DOF=Case.DOF;
elseif isfield(Case,'T')&size(Case.T,2)==size(o1.def,1)
  o1.DOF=Case.DOF;
end


% ---------------------------------------------------------------------------
% ---------------------------------------------------------------------------
% ---------------------------------------------------------------------------
function b=volume_load(model,Case,elt,r1,NNode)

b=zeros(size(model.DOF,1),1);
[EGroup,nGroup]=getegroup(elt);
node=model.Node;

if ~isequal(model.Elt,elt);
 model.Elt=elt;[Case,mdof]=fe_mknl('initnoconkeep',model);
elseif ~isfield(Case,'GroupInfo')
 [Case,mdof]=fe_mknl('initnoconkeep',model);
else; mdof=model.DOF; end
if isfield(r1,'dir')
 r2=elem0('VectFromDir',model,r1);
 jGroup=1; 
 % Some adjustments if there is a problem with the field at node
 [ElemF,i1,ElemP]= feutil('getelemf',elt(EGroup(jGroup),:),jGroup);
 switch ElemF
 case {'q4p','t3p','q5p','q8p','q8a','q9a','t3a','t3p','t6a','t6p'}
   r2=r2(1:2,:);
 case 'mitc4'
   r2(5,:)=0; 
 end
else; r2=elem0('VectFromDir',model,r1,model.DOF);
end

for j1=1:size(Case.GroupInfo,1);
  if isempty(Case.GroupInfo{j1,5}); Case.GroupInfo{j1,5}=zeros(1,100);end
end
if isfield(r1,'MatDes'); b=fe_mknl('assemble',model,Case,r1.MatDes,r2,NNode);
else; b=fe_mknl('assemble',model,Case,100,r2,NNode);
end
b=fe_c(mdof,model.DOF,b')';

return 

% ---------------------------------------------------------------------------
% ---------------------------------------------------------------------------
% ---------------------------------------------------------------------------
function b=surface_load(model,pl,il,r1,NNode)

elt=model.Elt; node=model.Node;

if isfield(r1,'eltsel') 
 [i1,elt]=feutil('findelt',node,model.Elt,[],r1.eltsel);
 if isempty(elt)
  error('Element selection for surface load selects no elements');
 end
end
model.Elt=elt; 
% Generate surface properties - - - - - - - - - - - - - - - - -
r3=unique(feutil('mpid',elt),'rows');r3(find(~any(r3,2)),:)=[];
for j1=1:size(r3,1);
 try;
  i2=[]; if ~isempty(model.il);i2=find(model.il(:,1)==r3(j1,2));end
  if isempty(i2);i2=size(model.il,1)+1;end
  model.il(i2,1:4)=[r3(j1,2) fe_mat('type','p_solid','SI',3) 0 3];
 end
end % define surface properties
[Case,mdof]=fe_mknl('initnocon',model);
b=zeros(size(mdof,1),1);

[EGroup,nGroup]=getegroup(elt);

% Obsolete selfacei call
if isfield(r1,'set')&~isempty(r1.set)
 ind=[];nind=[];
 FaceSet=stack_get(model,'set',r1.set,'GiveData');
 if ~isfield(FaceSet,'data')
  error(sprintf('Set ''%s'' not found',r1.set));
 end
 
elseif isa(r1.sel,'cell')&size(r1.sel,1)==1& ...
                           strcmp(comstr(r1.sel{1,2},-27),'selfacei')
 ind=r1.sel{1,4};
% Node numbers are given
elseif isnumeric(r1.sel); ind=r1.sel;
% Node selection command is given
elseif ischar(r1.sel) & ~isempty(r1.sel) 
  ind=feutil(['findnode' r1.sel],node,elt);
  if isempty(ind); error('No node selected for surface load'); end
% Node selection given as a Stack
elseif isa(r1.sel,'cell')&size(r1.sel,2)==4
  ind=feutil('findnode',node,elt,[],r1.sel);
else
 error('Currently surface load must be defined with selfacei selection');
end

eltid=feutil('eltidfix;',elt); 
if isempty(ind)
else %if ~isempty(ind) % Edge/face set defined by nodes
 nind(ind)=ind; if max(node(:,1))>length(nind);nind(max(node(:,1)))=0;end
 if isempty(ind) disp(r1); error('No node selected for surface load'); end
 for jGroup=1:nGroup
  EltConst=Case.GroupInfo{jGroup,8};
  % find the associated parent element
  [ElemF,i1,ElemP]= feutil('getelemf',elt(EGroup(jGroup),:),jGroup);
  if isfield(EltConst,'MatrixIntegrationRule') % *b family - - - - - - - - - -
   r2=feutil('selelt selface & innode',node,elt,[],ind);
   i2=find(r2(:,3)==r2(:,4)); if ~isempty(i2);
    i2=intersect(i2,feutil('findelt eltname quad4',node,r2));
    r2(end+1,1:6)=[Inf abs('tria3')];
    r2(end+[1:length(i2)],1:5)=r2(i2,[1:3 5 6]);r2(i2,:)=[];
   end
   if ~isfield(r1,'Elt'); r1.Elt=r2;
   else;r2.Elt(end+[1:size(r2,1)],1:size(r2,2))=r2;
   end
  else % MODULEF strategy with face numbers - - - - - - - - - - - - - - - -
   if i1<0 cEGI=[];  % Ignore groups for display only 
   else; cEGI = EGroup(jGroup)+1:EGroup(jGroup+1)-1;
   end

   if any(strcmp(ElemP, ...  % 2-D elements
    {'q4p','t3p','q5p','q8p','q8a','q9a','t3a','t3p','t6a','t6p'})) | ...
    any(strcmp(ElemF, ...  % 2-D elements
    {'q4p','t3p','q5p','q8p','q8a','q9a','t3a','t3p','t6a','t6p'}))
   i3=feval(ElemP,'edge')';dims=2;
   else; % 3D elements
     i3=feval(ElemP,'face')'; dims=3;
   end
   if ~isempty(cEGI);for j1=1:size(i3,2)
     i1=reshape(nind(elt(cEGI,i3(:,j1))),length(cEGI),size(i3,1));
     ind=find(all(i1,2)); % indices of included faces
     if ~isempty(ind) eltid(cEGI(ind),j1+1)=1; end
   end;end % Loop on faces
   [i3,i4]=find(eltid(:,2:end));
   FaceSet=struct('ID',0,'data',[eltid(i3,1) i4]);
   eltid=eltid(:,1);
  end % *b or MODULEF
 end % Loop on groups
end % Define face/edge set from nodes

% *b family strategy - - - - - - - - - - - - - - - - - - - - - - - - - - -
% WARNING THIS IS REALLY EXPERIMENTAL
if isfield(r1,'Elt'); 

% really dumb for large node number
spb=sparse(round(rem(r1.DOF,1)*100),fix(r1.DOF)+1,r1.def, ...
  100,max(ceil(max(r1.DOF)),max(node(:,1)))); % sparse reindexing of load 
i1=find(spb(:,1)); for j1=i1(:)' spb(j1,:)=spb(j1,1); end

i1=find(mdof>0);
nd=sparse(round(rem(mdof(i1),1)*100),fix(mdof(i1)),1:length(i1),100, ...
    max([fix(max(mdof(i1)))+1 max(model.Node(:,1))]));

 st1={'quad4','sload4';'quad8','sload8';'tria3','sload3';
      'tria6','sload6'};
 % need to have 

[EGroup,nGroup]=getegroup(r1.Elt);
% loop on load elements - - - - - - - - - - - - - - - - - - - - - - - - - 
for jGroup=1:nGroup
 cEGI = EGroup(jGroup)+1:EGroup(jGroup+1)-1;
 [ElemF,i1,ElemP]= feutil('getelemf',r1.Elt(EGroup(jGroup),:),jGroup);
 switch ElemF
 case 'quad4'; rule=integrules(ElemF,2); 
 case 'quadb'; rule=integrules(ElemF,9);
 case 'tria3'; rule=integrules(ElemF,-2);
 otherwise; rule=integrules(ElemF);
 end % case
 rule.VectMap=reshape(1:3*rule.Nnode,3,rule.Nnode)';
 inode=feval(ElemF,'node');
 NodePos=reshape(NNode(r1.Elt(cEGI,inode)'),length(inode),length(cEGI));
 rule.Nw=size(rule.w,1);
 if ~isfield(rule,'bas'); rule.bas=zeros(9,size(rule.N,1));end

 for jElt=1:length(cEGI)
  i1=NodePos(:,jElt); nodeE=Case.Node(i1(find(i1)),[5:7]);
  of_mk('buildndn',23,rule,nodeE);
  for jw=1:rule.Nw
     i1=r1.Elt(cEGI(jElt),inode);
     for ji=1:3; % surface loads and pressure
      r2=[spb(ji,i1)+spb(19,i1)*rule.bas(6+ji,jw)]*rule.NDN(:,jw); 
      if r2~=0; in1=full(nd(ji,i1));
       b(in1)=b(in1)+rule.NDN(:,jw)*r2*rule.jdet(jw)*rule.w(jw,4);
      end
     end
  end % jW
 end% jElt

end % group

else % MODULEF strategy - - - - - - - - - - - - - - - - - - - - - -

% build sparse matrix with loading
spb=sparse(round(rem(r1.DOF,1)*100),fix(r1.DOF)+1,r1.def, ...
  100,max(ceil(max(r1.DOF)),max(node(:,1)))); % sparse reindexing of load 
i1=find(spb(:,1)); for j1=i1(:)' spb(j1,:)=spb(j1,1); end
faces=sparse(FaceSet.data(:,1),FaceSet.data(:,2),1);

% loop on elements - - - - - - - - - - - - - - - - - - - - - - - - - 
for jGroup=1:nGroup

% find the associated parent element
[ElemF,i1,ElemP]= feutil('getelemf',elt(EGroup(jGroup),:),jGroup);
if i1<0 cEGI=[];  % Ignore groups for display only 
else cEGI = EGroup(jGroup)+1:EGroup(jGroup+1)-1;
 [i2,i3]=intersect(eltid(cEGI),FaceSet.data(:,1));
 cEGI=cEGI(i3);
end

if ~isempty(cEGI) % Some faces are loaded

  if any(strcmp(ElemF, ...  % 2-D elements
    {'q4p','t3p','q5p','q8p','q8a','q9a','t3a','t3p','t6a','t6p'}))
   FaceIndex=feval(ElemP,'edge')';
  else FaceIndex=feval(ElemP,'face')'; end

i4=feval(ElemP,'node');
i2=max(max(elt(cEGI,i4))); if i2>length(nind) nind(i2)=0;end
iopt=[0 0 length(feval(ElemF,'dof')) length(i4) 100 0 2 zeros(1,2*size(i3,2))];
iopt=int32(iopt);

vol=zeros(3,double(iopt(5)));
pointers=Case.GroupInfo{jGroup,2};
integ=Case.GroupInfo{jGroup,3};
constit=Case.GroupInfo{jGroup,4};

% Switch on element type - - - - - - - - - - - - - - - - - - - - - - - -

switch ElemF
case {'tetra4','hexa8','penta6','tetra10','hexa20','penta15'}
  % idof : help of_mk 
  % [coor nodeid], volf(3*noe), surf(3 ddl x nnof=4 x nface), 
  %                           pres(nnof=4 x nface)
if length(constit)<4; error('constit initialization failed');end

for jElt = 1:length(cEGI)

  nodeE=node(NNode(elt(cEGI(jElt),i4)),[5:7 1]);
  % define loaded faces
  point=pointers(:,jElt);point(5)=100; point(9:20)=0;

  % flags for loaded surface and pressure faces
  ind=full(faces(eltid(cEGI(jElt)),:));
  if length(ind)<size(FaceIndex,2); ind(size(FaceIndex,2))=0;end
  if length(ind)>size(FaceIndex,2); ind=ind(1:size(FaceIndex,2));end
  point(9:8+size(FaceIndex,2))=ind;
  point(9+size(FaceIndex,2):8+2*size(FaceIndex,2))=ind;

  % build the sur and pre values
  sur = zeros(3*size(FaceIndex,1),size(FaceIndex,2)); i5=size(sur,1);
  pre = zeros(size(FaceIndex,1),size(FaceIndex,2)); 
  i1=reshape(elt(cEGI(jElt),FaceIndex),size(FaceIndex,1),size(FaceIndex,2));
  ind=find(ind);
  sur(1:3:i5,ind)=reshape(spb(1,i1(:,ind)),size(i1,1),length(ind));
  sur(2:3:i5,ind)=reshape(spb(2,i1(:,ind)),size(i1,1),length(ind));
  sur(3:3:i5,ind)=reshape(spb(3,i1(:,ind)),size(i1,1),length(ind));

  pre(:,ind)=reshape(spb(19,i1(:,ind)),size(i1,1),length(ind));

  bi=of_mk(ElemF,int32(point),integ,constit,nodeE,[sur(:);pre(:)],vol);
  if ~isempty(bi)
    in1=double(Case.GroupInfo{jGroup,1}(:,cEGI(jElt)-EGroup(jGroup)))+1;
    in2=find(in1);in1=in1(in2);    b(in1,1)=b(in1,1)+bi(in2);
  end
end % loop on elements

case {'q4p','t3p','q5p','q8p','q8a','q9a','t3a','t3p','t6a','t6p'}
% 2-D elements - - - - - - - - - - - - - - - - - - - - - - - - - -
% (coor[noe,2], fomega[2,noe], fgamma[2,2*nb_edge] (fx,fy), ...
%  press[2*nb_edg], ...
%  no_ref[n_edge,2](int) (i,1) = 0 fgamma is 0, (i,2)=0 press is 0
%  alpha[3] (xx,xy,yy), theta[noe], car, iopt[1](int), be
%
% xxx [idof,bi]=of_mk(ElemF,[ndof noe 100 0 2 iopt no_ref], ...
%            node(i1,[5:7 1]),[fgamma(:),press(:);T],vol);
% The material properties are used for thermal problems.
if length(constit)<4; error('constit initialization failed');end

for jElt = 1:length(cEGI)

  nodeE=node(NNode(elt(cEGI(jElt),i4)),[5:7 1]);

  % define loaded faces
  point=pointers(:,jElt);point(5)=100; point(9:20)=0;

  % flags for loaded surface and pressure faces
  ind=full(faces(eltid(cEGI(jElt)),:));
  if length(ind)<size(FaceIndex,2); ind(size(FaceIndex,2))=0;end
  if length(ind)>size(FaceIndex,2); ind=ind(1:size(FaceIndex,2));end
  point(9:8+size(FaceIndex,2))=ind;
  point(9+size(FaceIndex,2):8+2*size(FaceIndex,2))=ind;

  % build the sur and pre values
  sur = zeros(2*size(FaceIndex,1),size(FaceIndex,2)); i5=size(sur,1);
  pre = zeros(size(FaceIndex,1),size(FaceIndex,2)); temp=pre;
  i1=reshape(elt(cEGI(jElt),FaceIndex),size(FaceIndex,1),size(FaceIndex,2));
  ind=find(ind);
  sur(1:2:i5,ind)=reshape(spb(1,i1(:,ind)),size(i1,1),length(ind));
  sur(2:2:i5,ind)=reshape(spb(2,i1(:,ind)),size(i1,1),length(ind));
  pre(:,ind)=reshape(spb(19,i1(:,ind)),size(i1,1),length(ind));
  temp(:,ind)=reshape(spb(20,i1(:,ind)),size(i1,1),length(ind));
  bi=of_mk(ElemF,int32(point),integ,constit,nodeE, ...
    [sur(:);pre(:);temp(:)],vol);
 
  if any(~finite(bi)) error(1);end

  if ~isempty(bi)
    in1=double(Case.GroupInfo{jGroup,1}(:,cEGI(jElt)-EGroup(jGroup)))+1;
    in2=find(in1);in1=in1(in2);    b(in1,1)=b(in1,1)+bi(in2);
  end

end % loop on elements

case {'beam1','bar1','beam3','mass1','celas','rigid','mass2','cbush'}
% 1-D elements - - - - - - - - - - - - - - - - - - - - - - - - - -
 sdtw('Surface Load not defined for : %s (group %i) \n',ElemF,jGroup);
 idof=[];bi=[];

otherwise 

 try;
  % xxx call for in element surface load here
 end

 
 sdtw('Surface Load not yet supported for : %s \n',ElemF);
 idof=[];bi=[];

end;end % of element type and cEGI is not empty

end % of loop on jGroup
end % strategy

b=fe_c(model.DOF,mdof,b')';
