function [out,out1]=elem0(CAM,varargin);

% Default element function used for elements that conform strictly to an
% OpenFEM standard
%
% See also : doc elem0


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

if comstr(CAM,'cvs')
 out='$Revision: 1.67 $  $Date: 2006/05/17 15:47:43 $'; return;
end

if nargin==0;doc('elem0');return;
else;  [CAM,Cam]=comstr(CAM,1);end


if comstr(Cam,'callmat_og');

 out='[k1,m1]=elem0(''mat_og'',jElt,NodePos,Case.Node,pointers,integ,constit,gstate,elmap,InfoAtNode,EltConst,def);';

% Elements of the Mat_Og family use a [1.01 ... N.01 1.02 ... N.02 ...]
% internal DOF sorting because this is appropriate for BLAS calls
% The standard numbering returned by elements is [1.01 1.02 ... 2.01 2.02 ...]
elseif comstr(Cam,'elmapmat_og');

  ind=varargin{1};
  i2=prod(ind); i1=reshape(1:i2,ind(2),ind(1))';
  out=reshape(1:i2^2,i2,i2); out(i1,i1)=out;

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
% standard GroupInit call for og elements
elseif comstr(Cam,'groupinitog');[CAM,Cam]=comstr(CAM,12);

 if isempty(Cam);
  out=sprintf('[Case.GroupInfo{jGroup,8},pointers]=%s(''Constants'',pointers,integ,constit);',varargin{1});
 elseif comstr(Cam,'shell')
 out=sprintf('[Case.GroupInfo{jGroup,8},pointers,Case.GroupInfo{jGroup,7}]=%s(''Constants'',pointers,integ,constit,model,Case,cEGI);',varargin{1});
end

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
% standard DOF building for variable field elements
elseif comstr(Cam,'dofcall');
 out='[i2,i3]=p_solid(''BuildDof'',model,cEGI,nd,ElemF);';
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
% standard Matcall that allows selection of mat_og and mat_of strategies
elseif comstr(Cam,'matcall');

 if varargin{1}(end,1)==-9999 % Elements issued from Modulef in of_mk_sub.c
  out='mat_of'; out1=1; % call SymFlag
 else % Newer elements handled by of_mk.c MatrixIntegration
   out='mat_og'; out1=0; % mat_og and non symmetric
 end
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
elseif comstr(Cam,'rhsmat_og'); out='rhs_og';

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
% Matlab implementation of single element matrix integration call. 
% This is only meant for development purposes
elseif comstr(Cam,'rhs_og');

 [DofPos,NodePos,node, ...
      pointers,integ,constit,gstate, ...
      elmap,InfoAtNode,EltConst,def,out]=deal(varargin{:});

for jElt=1:size(DofPos,2);

 point=pointers(:,jElt);
 nodeE=node(NodePos(:,jElt),5:7);
 of_mk('buildndn',point(4),EltConst,nodeE);

 Nw=EltConst.Nw;Nnode=EltConst.Nnode; Ndof=Nnode*length(EltConst.DofLabels);
 if point(4)==3; Ndim=3; else; Ndim=2;end

 if ~isfield(EltConst,'VectMap'); 
  i1=reshape(1:Ndof,Ndof/Nnode,Nnode)'; EltConst.VectMap=i1(:);
 end
 inde=double(DofPos(:,jElt))+1; in2=find(inde>0);
 inde=inde(EltConst.VectMap);

switch double(point(5))
case 100  % Volume load (not proportional to density)

 defe=def(1:Ndim,NodePos(:,jElt));
 for jW=0:Nw-1
  r1=EltConst.NDN(:,ones(Ndim,1)*(jW+1))*diag((defe*EltConst.NDN(:,jW+1))* ...
   EltConst.jdet(jW+1)* EltConst.w(jW+1,4));
  out(inde(in2))=out(inde(in2))+r1(in2);
 end % loop on jW

case 101  % Pressure load

 keyboard

case 102 % inertia load (not proportional to density)

 if isfield(def,'def'); 
  defe=reshape(def.def(double(DofPos(:,jElt))+1), ...
   length(EltConst.DofLabels),EltConst.Nnode);
 else; defe=def(:,NodePos(:,jElt));
 end
 for jW=0:Nw-1
  r1=EltConst.NDN(:,[jW jW jW]+1)*diag((defe*EltConst.NDN(:,jW+1))* ...
   EltConst.jdet(jW+1)* EltConst.w(jW+1,4)*constit(1+point(7)));
  % it is assumed (NOT CHECKED) that rho(1) contains density
  out(inde(in2))=out(inde(in2))+r1(in2);
 end % loop on jW

case 103 % Initial stress load (F = \int \sigma_0 \epsilon dv)

 % row, NDN, DDL, NwStart NwRule
 r2=EltConst.StrainDefinition{1}; Nstrain=length(EltConst.StrainLabels{1});
 Nnode=EltConst.Nnode;
 for jW=0:Nw-1
  for ji=0:size(r2,1)-1 % Strain definition rows
   jj=(r2(ji+1,2)-1)*Nw+jW; % NDN column
   jk=(r2(ji+1,3)-1)*Nnode; % DOF block
   ib=double(DofPos(EltConst.VectMap(jk+[1:Nnode]),jElt))+1;
   r1=EltConst.NDN(:,jj+1)*gstate(jW*Nstrain+r2(ji+1,1),jElt);
   out(ib)=out(ib)+r1*EltConst.jdet(jW+1)* EltConst.w(jW+1,4);
  end

 end % loop on jW

otherwise error('Not a supported load type');
end
end % loop on elements of the group

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
elseif comstr(Cam,'mat_og');

  varg=varargin; 
  if ~isa(varg{4},'int32'); 
   varg{4}=int32(varg{4});evalin('caller','pointers=int32(pointers);');
  end
  if ~isa(varg{2},'int32'); 
   varg{2}=int32(full(varg{2}));evalin('caller','pointers=int32(pointers);');
  end
  [out,out1]=of_mk('matrixintegration',varg{:});
  if ~isempty(varg{8});
   out=out(varg{8});if ~isempty(out1); out1=out1(varg{8});end
  end

% -----------------------------------------------------------------------
elseif comstr(Cam,'matrixintegration');
% -----------------------------------------------------------------------
% -----------------------------------------------------------------------
% -----------------------------------------------------------------------
% Matlab implementation of single element matrix integration call. 
% THIS IS ONLY MEANT FOR DEVELOPMENT PURPOSES
% -----------------------------------------------------------------------
% -----------------------------------------------------------------------
% -----------------------------------------------------------------------


[DofPos,NodePos,node, ...
      pointers,integ,constit,gstate, ...
      elmap,InfoAtNode,EltConst,def,k,AsmPoint]=deal(varargin{:});

if isempty(k); 
 jElt=DofPos; DofPos=evalin('caller','DofPos');
elseif isa(k,'struct');
 cEGI=k.cEGI; out1=k.Ener;
end

Nw=EltConst.Nw;Nnode=EltConst.Nnode; Ndof=Nnode*length(EltConst.DofLabels);
NdofPerField=size(EltConst.NDN,1);
if ~isfield(EltConst,'material'); EltConst.material='';end
try; Ncondense=EltConst.Ncondense;catch; Ncondense=0; end

for jElt=1:size(DofPos,2) % loop on elements

point=pointers(:,jElt);
nodeE=node(NodePos(:,jElt),5:7);
of_mk('buildndn',point(4),EltConst,nodeE);

if any(point(5)==[5 6])  % 3D Geometric nonlinear stiffness matrix

 Debug=0; try; Debug=evalin('base','Debug');end

 defe=def(double(DofPos(:,jElt))+1,1); %1xyz 2xyz, ...
 ke=zeros(Ndof);

 % Step 1 : build the local gradient vector at integration points
 % lestate (9 * ninteg) : here u_{i,j}

% constants 
 ind_ts_eg=[1 6 5 6 2 4 5 4 3]; % e_11 22 33 23 31 12
 ind_ts_egt=[1 6 5 6 2 4 5 4 3]'; % e_11 22 33 23 31 12
 if isempty(EltConst.material);
   if isequal(EltConst.DofLabels(:)',{'u','v','w'});
    EltConst.material='Elastic3DNL';
   end
 end
switch EltConst.material
case 'hyperelastic'
 d2I2dcdc=[0 1 1 0 0 0;1 0 1  0 0 0;1 1 0 0 0 0;0 0 0 -.5 0 0;
           0 0 0 0 -.5 0;0 0 0 0 0 -.5];
 ci_ts_eg=[1 5 9 8 7 4];ci_ts_egt=ci_ts_eg';
 dI1dc=eye(3);

 for jW=0:Nw-1; % HERE IS THE LOOP ON INTEGRATION POINTS 

   r1=EltConst.NDN(:,[Nw 2*Nw 3*Nw]+jW+1);
   U_ij=reshape(defe,3,Nnode)*r1; F_ij=U_ij+eye(3);

   C=F_ij'*F_ij;
   I=[trace(C) (trace(C)^2-sum(C(:).^2))/2 det(C)]; % invariants
   CI=inv(C);

   d2I3dcdc=[0 C(9) C(5) -C(8) 0    0;
             C(9) 0 C(1)    0 -C(7) 0;
             C(5) C(1) 0    0  0   -C(4);
             -C(8) 0 0     -C(1)/2  C(4)/2 C(7)/2
              0 -C(7) 0     C(2)/2 -C(5)/2 C(8)/2
              0  0   -C(4)  C(3)/2 C(6)/2 -C(9)/2];
   dI2dc=I(1)*eye(3)-C; dI3dc=I(3)*CI;

 [dWdI,d2WdI2]=EnHeart(integ,constit,I);
  % xxx call Marina xxx
   % Equation (15)
   Sigma=2*(dWdI(1)*dI1dc+dWdI(2)*dI2dc+dWdI(3)*dI3dc);

  if Debug==1 % Test of rivlin cube
   L=evalin('base','L');
   test(1)=norm(diag(L)+eye(3)-F_ij);
   test(2)=norm(C-(eye(3)+diag(L))^2);
   r1=0;for ji=1:3;for jj=1:ji-1; r1=r1+(1+L(ji))^2*(1+L(jj))^2;end;end
   test(3)=norm(I-[sum((1+L).^2) r1 prod((1+L).^2)]);
   test(4)=norm(dI2dc- (I(1)*eye(3)-diag((1+L).^2)));
   test(5)= norm(dI3dc- (I(3)*diag((1+L).^(-2))));
   for ji=1:3; 
    r1(ji)=2*(dWdI(1)+dWdI(2)*(I(1)-(1+L(ji))^2) + dWdI(3)*I(3)*(1+L(ji))^-2);
   end
   test(6)=norm(Sigma-diag(r1))/norm(r1);
   if norm(test,'inf')>1e-5; error('Error on state check');end
   assignin('base','Sigma',Sigma);

  end % Test of rivlin cube

   % Equation (16)
   d2wde2=dWdI(2)*d2I2dcdc+dWdI(3)*d2I3dcdc + ...
    d2WdI2(1)*dI1dc(ci_ts_egt)*dI1dc(ci_ts_eg) + ...
    d2WdI2(2)*dI2dc(ci_ts_egt)*dI1dc(ci_ts_eg) + ...
    d2WdI2(3)*dI3dc(ci_ts_egt)*dI1dc(ci_ts_eg) + ...
    d2WdI2(4)*dI1dc(ci_ts_egt)*dI2dc(ci_ts_eg) + ...
    d2WdI2(5)*dI2dc(ci_ts_egt)*dI2dc(ci_ts_eg) + ...
    d2WdI2(6)*dI3dc(ci_ts_egt)*dI2dc(ci_ts_eg) + ...
    d2WdI2(7)*dI1dc(ci_ts_egt)*dI3dc(ci_ts_eg) + ...
    d2WdI2(8)*dI2dc(ci_ts_egt)*dI3dc(ci_ts_eg) + ...
    d2WdI2(9)*dI3dc(ci_ts_egt)*dI3dc(ci_ts_eg); 
   d2wde2=d2wde2*4;
    %d2wde2(ind_ts_eg,ind_ts_eg)
   if Debug>0 & min(eig(d2wde2))<0 error('Negative constitutive law');end
   % dW/de = dW/di di/de

  %DD(ci_ts_eg,ci_ts_eg)
   Mecha3DInteg(ke,F_ij,d2wde2,Sigma, ...
    def,EltConst.w(:,4),EltConst.jdet,EltConst.NDN,Nnode,Ndof,Nw,jW)

  end % END THE LOOP ON INTEGRATION POINTS 

  % assignin('base','d2wde2_hyper',d2wde2);
case 'Elastic3DNL'

 for jW=0:Nw-1; % HERE IS THE LOOP ON INTEGRATION POINTS 

   r1=EltConst.NDN(:,[Nw 2*Nw 3*Nw]+jW+1);
   U_ij=reshape(defe,3,Nnode)*r1; F_ij=U_ij+eye(3);

  i1=EltConst.ConstitTopology{1}; %i1=i1(ind_ts_eg,ind_ts_eg);
  d2wde2=constit(double(i1)+double(point(6+1)));
  Sigma=d2wde2(ind_ts_eg,ind_ts_eg)*U_ij(:);
  if ~isempty(gstate)
   Sigma=Sigma+gstate(jW*6+ind_ts_eg,jElt);
  end
  % assignin('base','d2wde2_elas',d2wde2);
  %DD(ci_ts_eg,ci_ts_eg)
  Mecha3DInteg(ke,F_ij,d2wde2,reshape(Sigma,3,3), ...
   def,EltConst.w(:,4),EltConst.jdet,EltConst.NDN,Nnode,Ndof,Nw,jW)

 end % END THE LOOP ON INTEGRATION POINTS 
case 'fs_matrix'  % Follower pressure (fluid/structure coupling)

 i1=[-Nnode+1:0]';Axr=zeros(3);Axs=zeros(3);
 ke=zeros(Nnode*4);F=def(:,2);fe=zeros(Nnode*4,1);

 for jW=0:Nw-1; % HERE IS THE LOOP ON INTEGRATION POINTS 

   
   p=EltConst.N(jW+1,:)*defe(EltConst.VectMap(end+i1));
   eDofPos=EltConst.VectMap([ ...
     EltConst.Nnode*1+i1 EltConst.Nnode*2+i1 EltConst.Nnode*3+i1]);
%evalin('caller',['fe_c(model.DOF([' sprintf('%i ',double(DofPos(eDofPos(:),1))+1) ']))'])
   x=defe(eDofPos)+nodeE;
   axr=EltConst.Nr(jW+1,:)*x;
   axs=EltConst.Ns(jW+1,:)*x;
   Axr([6 7 2])=axr;Axr([8 3 4])=-axr; 
   dq=zeros(4*Nnode,3);
   dq(Nnode+i1,1)=EltConst.N(jW+1,:)';
   dq(Nnode*2+i1,2)=EltConst.N(jW+1,:)';
   dq(Nnode*3+i1,3)=EltConst.N(jW+1,:)';
   dqr=zeros(3,4*Nnode);
   dqr(1,Nnode+i1)=EltConst.Ns(jW+1,:);
   dqr(2,Nnode*2+i1)=EltConst.Ns(jW+1,:);
   dqr(3,Nnode*3+i1)=EltConst.Ns(jW+1,:);
   ke=ke+dq*(EltConst.w(jW+1,4)*Axr)*dqr;
   Axr([6 7 2])=axs;Axr([8 3 4])=-axs; 
   dqr(1,Nnode+i1)=EltConst.Nr(jW+1,:);
   dqr(2,Nnode*2+i1)=EltConst.Nr(jW+1,:);
   dqr(3,Nnode*3+i1)=EltConst.Nr(jW+1,:);
   ke=ke+p*dq*(-EltConst.w(jW+1,4)*Axr)*dqr;
  % RHS computation : P x,r ^ x,s \delta v
  fe=p*(EltConst.w(jW+1,4)*cross(axr,axs)'*EltConst.N(jW+1,:))';fe=fe(:);
  i2=double(DofPos(eDofPos,jElt))+1; F(i2)=F(i2)+fe;
 end % END THE LOOP ON INTEGRATION POINTS 
 sp_util('setinput',def,F,int32(1:length(def)),int32(2));

otherwise;
  error('Not a valid EltConst.material');
end 

% Linear fluid structure coupling (mass[2] and stiffness[1] matrices)
elseif isequal(EltConst.material,'fs_matrix') % - - - - - - - - - - - - -

 ke=zeros(Ndof);
 for jW=0:Nw-1; % HERE IS THE LOOP ON INTEGRATION POINTS 
   % ({u}.{n}) p ds
   for ji=0:2; 
     coef=EltConst.bas(7+ji,jW+1)*EltConst.jdet(jW+1)*EltConst.w(jW+1,4);
     if coef; jj=Nnode*ji+[1:Nnode]; jk=Nnode*3+[1:Nnode];
      if point(5)==1;
       ke(jj,jk)=ke(jj,jk)-coef*EltConst.N(jW+1,:)'*EltConst.N(jW+1,:);
      elseif point(5)==2
       ke(jk,jj)=ke(jk,jj)+coef*EltConst.N(jW+1,:)'*EltConst.N(jW+1,:);
      end
     end
   end
 end

else % standard methodology for MatrixIntegration rule - - - - - - - -

rule=double(EltConst.MatrixIntegrationRule{point(5)});
% check constit topology
%DD=double(EltConst.ConstitTopology{point(5)});DD(find(DD))=constit(DD(find(DD)));disp([jElt DD(1)])
ke=zeros(double(integ(3))+Ncondense);
for j2=1:size(rule,1) % terms of rule
  wstep=rule(j2,6);
  for jw=1:double(rule(j2,7));
    i1=rule(j2,5); if i1<0; coef=-1; else;coef=1;end
    i1=abs(i1);i2=jw+rule(j2,8);
    coef=coef*EltConst.w(i2,4)*EltConst.jdet(i2)* ...
     constit((jw-1)*wstep+i1+1+double(point(7)));
    if coef
       for j4=1:NdofPerField
         iblock=(rule(j2,2)+j4-1)*size(ke,1)+rule(j2,1)+1;
            for j3=1:NdofPerField
             ke(iblock)=ke(iblock)+coef* ...
             EltConst.NDN(j3+(jw+rule(j2,3)-1)*NdofPerField)* ...
             EltConst.NDN(j4+(jw+rule(j2,4)-1)*NdofPerField);
             iblock=iblock+1;
            end
       end
    end
  end % jw
end % rule

% Volume integration if defined
if size(def,2)>1&isfield(EltConst,'RhsDefinition')
  rrule=double(EltConst.RhsDefinition);
  F=def(:,2); 
  Ndof=size(InfoAtNode,1)/EltConst.Nnode;
  jdet=EltConst.jdet;w=EltConst.w(:,4);Nnode=EltConst.Nnode;
  for jRule=1:size(rrule,1)
   r3=rrule(jRule,:);
   switch r3(1)
   case 101 % Volume load defined in InfoAtNode
   %101(1) InfoAtNode1(2) InStep(3) NDNOff1(4) 
   % FDof1(5) NDNCol(6) NormalComp(7)
   % w1(8) nwStep(9)
   if isempty(InfoAtNode); break;end
    for jw=r3(8)+[0:r3(9)-1]
     r1=EltConst.NDN(:,jw+r3(4)+1)'* ...
       InfoAtNode(r3(2)+[0:Nnode-1]*r3(3)+1,jElt);
     r1=r1*jdet(jw+1)*w(jw+1);
     if r3(7)>=0; r1=r1*EltConst.bas(r3(7)+1,jw+1);end
     in1=double(DofPos(EltConst.VectMap(Nnode*r3(5)+[1:Nnode]),jElt))+1;
     F(in1)=F(in1)+  EltConst.NDN(:,jw+r3(6)+1)*r1;
    end % jw
   otherwise; sdtw('_nb','% Rhs type not supported');
   end
  end % jRule

  sp_util('setinput',def,F,int32(1:length(def)),int32(2));
end % volume integration if needed

end % matrix type

if Ncondense; % possibly condense matrix if needed
 i1=double(EltConst.CondenseInd)+1; ke=ke(i1,i1);
 ind=1:size(ke,1)-Ncondense;cind=size(ke,1)-Ncondense+1:size(ke,1);
 ke=ke(ind,ind)-ke(ind,cind)* (ke(cind,cind)\ke(cind,ind));
end

%       ke=of_mk('matrixintegration',jElt,NodePos,Case.Node, ...
%       pointers,integ,constit,gstate, ...
%       elmap,InfoAtNode,EltConst,def);
%       i1=[Case.DofPerElt(jGroup);SymFlag;0;DofPos(:,jElt)];
%       of_mk('asmsparse',k,int32(i1),ke,elmap);
 if isa(k,'struct')
  in1=double(DofPos(EltConst.VectMap,jElt))+1; in2=(find(in1)); in1=in1(in2);
  vect=zeros(size(DofPos,1),1); if all(EltConst.jdet<0); ke=-ke;end
  for j2=1:size(def,2)
   vect(in2)=def(in1,j2); out1(cEGI(jElt),j2)=real(vect'*ke*vect);
  end
  keyboard
 elseif ~isempty(k); % assemble element matrix if needed
  i1=[AsmPoint;DofPos(:,jElt)]; of_mk('asmsparse',k,int32(i1),ke,elmap);
 end
 out=ke;
end % loop on elements
if isa(k,'struct')
  sp_util('setinput',k.Ener,out1, ...
   int32(1:size(k.Ener,1)),int32(1:size(k.Ener,2)));
end

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
% Matlab implementation of single element stress computation call. 
elseif comstr(Cam,'stress_og');

[NodePos,jGroup,Case,def,pointers,integ,constit,EltConst]= ...
  deal(varargin{:});

node=Case.Node;

if comstr(version,'7'); rule=EltConst.StressRule{1};
else; rule=double(EltConst.StressRule{1});
end
out=[];

for jElt=1:size(NodePos,2)

 i2=NodePos(:,jElt);nodeE=node(i2,[5:7 1]);
 point=int32(pointers(:,jElt));Ns=length(EltConst.StrainLabels{1});
 if ~isempty(out)
 elseif point(4)==3
  out=zeros(Ns*EltConst.Nw,size(def.def,2)*size(NodePos,2));
 elseif point(4)==2
  out=zeros(Ns*EltConst.Nw,size(def.def,2)*size(NodePos,2));
 end

try;
 ke=of_mk('StressObserve',EltConst,int32(rule),constit,nodeE,point);
 Ndof=double(point(4)); nw=EltConst.Nw; Nnode=EltConst.Nnode;
catch
 point=double(point);
 nw=EltConst.Nw; Nnode=EltConst.Nnode;
 of_mk('buildndn',point(4),EltConst,nodeE);
 if point(4)==3;     Ndof=3;
 elseif point(4)==2; Ndof=2;
 else; error('Not a valid StressObserve case');
 end
 ke=zeros(Ns*nw,Nnode*Ndof);
 for jConst=1:size(rule,1) % loop on elements of rule
  coef=constit(rule(jConst,3)+1+point(7));
  for jN=1:Nnode;for jW=1:nw;
    i1=(jW-1)*rule(jConst,5)+rule(jConst,1)+1;
    i2=rule(jConst,4)+(jN-1)*Ndof+1;
    ke(i1,i2)= ke(i1,i2)+ ...
     EltConst.NDN(jN,rule(jConst,2)+jW)*coef;
  end;end
 end
end % m-file implementation

 state=zeros(Ndof*Nnode,1);
 i3=double(Case.GroupInfo{jGroup,1}(:,jElt))+1;i4=find(i3);i3=i3(i4);
 state(i4,1:size(def.def,2))=def.def(i3,:);

 out(:,jElt:size(NodePos,2):end)=ke*state;

end % loop on elements
out1=nw;

% -----------------------------------------------------------------------
elseif comstr(Cam,'mooney');

 [out,out1]=EnHeart(varargin{:});

% -----------------------------------------------------------------------
elseif comstr(Cam,'initgstate');

 model=varargin{1};Case=varargin{2};
 [EGroup,nGroup]=getegroup(model.Elt);
 out1=struct('DOF',model.DOF,'def',zeros(length(model.DOF),2), ...
   'lab',{{'CurrentState';'RHS'}});
 for j1=1:size(Case.Stack,1)
  switch comstr(Case.Stack{j1,1},-27)
  case 'fvol'  % Place the volume loads at the beginning of InfoAtNode
   data=Case.Stack{j1,3};
   if ~ischar(data.sel); 
    sdtw('_nb','Ignoring FVol(%s) defined by elements',Case.Stack{j1,2});
    break;
   end
   [i1,elt]=feutil(horzcat('findelt',data.sel),model);
   % field specified at nodes
   if isfield(data,'def')&isfield(data,'DOF')
     if max(fix(data.DOF))==0; % constant field
      r1=zeros(size(model.DOF));
      for j2=1:length(data.DOF); 
        r1(fe_c(model.DOF,data.DOF(j2),'ind'))=data.def(j2);
      end
     else; % field specified at DOFs
      r1=fe_c(model.DOF,data.DOF,data.def(:,1)')';
     end
     for jGroup=1:nGroup;if isa(Case.GroupInfo{jGroup,8},'struct')
       % fill in the volume forces
       cEGI = intersect([EGroup(jGroup)+1:EGroup(jGroup+1)-1]',i1);
       if ~isempty(cEGI)
        gstate=r1(double(Case.GroupInfo{jGroup,1})+1);
        cEGI=cEGI-EGroup(jGroup);
        if any(size(Case.GroupInfo{jGroup,7})<size(gstate))
         Case.GroupInfo{jGroup,7}(size(gstate,1),size(gstate,2))=0;
        end
        Case.GroupInfo{jGroup,7}(1:size(gstate,1),cEGI)= ...
        Case.GroupInfo{jGroup,7}(1:size(gstate,1),cEGI)+gstate(:,cEGI);
       end
     end;end  % fill in the volume forces
   elseif isfield(data,'dir');
    r2=elem0('VectFromDir',model,data);
    for jGroup=1:nGroup;if isa(Case.GroupInfo{jGroup,8},'struct')
       % fill in the volume forces
       cEGI = intersect([EGroup(jGroup)+1:EGroup(jGroup+1)-1]',i1);
       if ~isempty(cEGI)
        [ElemF,EGID]= feutil('getelemf',model.Elt(EGroup(jGroup),:),jGroup);
        inode=fe_super('node',ElemF);
        NodePos=model.Elt(cEGI,inode)'; 
        if any(Case.GroupInfo{jGroup,2}(4,1)==[2 12]); ind=1:2;
        else; ind=1:3
        end
        gstate=reshape(r2(ind,NodePos),length(ind)*size(NodePos,1), ...
          size(NodePos,2));
        cEGI=cEGI-EGroup(jGroup);
        if any(size(Case.GroupInfo{jGroup,7})<size(gstate))
         Case.GroupInfo{jGroup,7}(size(gstate,1),size(gstate,2))=0;
        end
        Case.GroupInfo{jGroup,7}(1:size(gstate,1),cEGI)= ...
        Case.GroupInfo{jGroup,7}(1:size(gstate,1),cEGI)+gstate(:,cEGI);
       end
    end;end % jGroup

   end
  end
 end
 out=Case;

% -----------------------------------------------------------------------
% Builds a matrix whose column give components of a field at node
elseif comstr(Cam,'vectfromdir')

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

node=model.Node;
% Build the field at nodes (r2) based in .dir field
if isnumeric(r1.dir)&prod(size(r1.dir))==3 
 r2=r1.dir(:)'; r2=r2(ones(size(model.Node,1),1),:)';
elseif isa(r1.dir,'cell')&length(r1.dir)==3
 r2=zeros(size(model.Node,1),3);
 x=model.Node(:,5); y=model.Node(:,6); z=model.Node(:,7); 
 for j1=1:length(r1.dir) 
  if isnumeric(r1.dir{j1}); r2(:,j1)=r1.dir{j1};
  elseif any(r1.dir{j1}=='x'|r1.dir{j1}=='y'|r1.dir{j1}=='z')
    r3=eval(r1.dir{j1}); r3=r3(:);
    if size(r3,1)~=size(r2,1)&size(r3,1)~=1 
     error('String giving an inconsistent evaluation of volume load at nodes');
    end
    r2(:,j1)=r3(:);
  else
   try; r2(:,j1)=str2num(r1.dir{j1});
   catch
    try; r2(:,j1)=eval(r1.dir{j1}); 
    catch; r2(:,j1)=feval(r1.dir{j1},node);
    end
   end
  end
 end
 r2=r2';
else
  error('Not a supported volume load direction format');
end
out=r2;

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

carg=1;
i1=varargin{carg};carg=carg+1;
ke=varargin{carg};carg=carg+1;
if carg<nargin; me=varargin{carg};carg=carg+1;else; me=[];end

for j1=size(ke,1):-1:i1
 ke(j1,j1)=1/ke(j1,j1); r1=ke(j1,j1); % set pointer kcc_i
 ind=1:j1-1;
 ke(ind,j1)=-ke(ind,j1)*r1; vect=ke(ind,j1); % set pointer k_ic
 ke(ind,ind)=ke(ind,ind)-1/r1*vect*vect'; % do a dger call
 if ~isempty(me); 
  me(ind,ind)=me(ind,ind)-vect*me(j1,ind)-me(ind,j1)*vect'+ vect*me(j1,j1)*vect';
 end
end
out=ke(1:i1-1,1:i1-1); out1=me(1:i1-1,1:i1-1);

% -----------------------------------------------------------------------
else; sdtw('_nb','''%s'' not supported',CAM);
end


% -----------------------------------------------------------------------
function [dWdI,d2WdI2]=EnHeart(integ,constit,I)

%C1=0.3MPa, C2=0.2MPa, K=0.3MPa

% Cenerg : C1 C2 K
%if length(constit)<3;constit=[25. 0.25 1274.];end
%if length(constit)<3;constit=[0 1 3 2 3];end
if length(constit)<3;constit=[0 .3 .2 .3];end


if (constit(1) == 0)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%hperelastic type 0: C1*(J1-3)+C2*(J2-3)+K*(J3-1)^2%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

dWdI(1) = constit(2)*I(3)^(-1./3.);
dWdI(2) = constit(3)*I(3)^(-2./3.);
dWdI(3) = - 1./3.* constit(2)*I(1)*I(3)^(-4./3.) ...
         - 2./3.* constit(3)*I(2)*I(3)^(-5./3.) ...
             +   constit(4)*(1-I(3)^(-1/2)) ;

d2WdI2=[0 0 -1./3.*constit(2)*I(3)^(-4./3.) ;
       0 0  -2./3.*constit(3)*I(3)^(-5./3.);
      -1./3.*constit(2)*I(3)^(-4./3.)  -2./3.*constit(3)*I(3)^(-5./3.) ...
       4/9*constit(2)*I(1)*I(3)^(-7./3.)...
             + 10./9.* constit(3)*I(2)*I(3)^(-8./3.)...
             + 1/2 * constit(4)*I(3)^(-3/2)];

else
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%hyperelastic type 1: C1*(J1-3)+C2*(J2-3)+K*(J3-1)-K*ln(J3)%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

dWdI(1) = constit(2)*I(3)^(-1./3.);
dWdI(2) = constit(3)*I(3)^(-2./3.);
dWdI(3) = - 1./3.* constit(2)*I(1)*I(3)^(-4./3.) ...
          - 2./3.* constit(3)*I(2)*I(3)^(-5./3.) ...
          + 0.5 * constit(4)*I(3)^(-1./2.) ...
          - 0.5*constit(4)*I(3)^(-1);

d2WdI2=[0 0 -1./3.*constit(2)*I(3)^(-4./3.) ;
        0 0  -2./3.*constit(3)*I(3)^(-5./3.);
       -1./3.*constit(2)*I(3)^(-4./3.)  -2./3.*constit(3)*I(3)^(-5./3.) ...
        4./9.*constit(2)*I(1)*I(3)^(-7./3.)...
              + 10./9.* constit(3)*I(2)*I(3)^(-8./3.)...
              - 1./4.*constit(4)*I(3)^(-3./2.) ...
              + 0.5*constit(4)*I(3)^(-2)];
end;

% -------------------------------------------------------------------------
function Mecha3DInteg(ke,F_ij,d2wde2,Sigma,def,w,jdet,NDN,Nnode,Ndof,Nw,jW)




try
   %ind_ts_eg=[1 6 5 6 2 4 5 4 3];
   %d2wde2=d2wde2*0;Sigma=evalin('base','Sigma0');Sigma=Sigma(ind_ts_eg);
  of_mk('Mecha3DInteg',ke,F_ij,d2wde2,Sigma,w,jdet,NDN,int32([Nnode Ndof Nw jW]));
  %k1=ke;ke=ke*0; error(1);
catch

   ind_ts_eg=[1 6 5 6 2 4 5 4 3];
   % energy : (dW_n,j F_ni) (d2wde)_(ij kl) (dW_m,l F_mk)
   % strain indices n,j m,l, DD indices  i k
   % DD= [F_ni (d2wde)_(ij kl) F_mk]_(nj ml)
   DD=zeros(9); 
   for jn=0:2; for jj=0:2; for jm=0:2; for jl=0:2; 
     for ji=0:2; for jk=0:2
      DD(jn+3*jj+9*jm+27*jl+1)=DD(jn+3*jj+9*jm+27*jl+1)+ ...
       F_ij(jn+ji*3+1)* ...
       d2wde2(ind_ts_eg(ji+3*jj+1)+6*(ind_ts_eg(jk+3*jl+1)-1))* ...
       F_ij(jm+jk*3+1);
     end;end
     if jn==jm; 
       DD(jn+3*jj+9*jm+27*jl+1)=DD(jn+3*jj+9*jm+27*jl+1)+ ...
        Sigma(jj+3*jl+1);
     end
   end;end;end;end

try;

    for ji=0:2;for jj=0:2;for jk=0:2;for jl=0:2;
    coef=DD(ji+jj*3+jk*9+jl*27+1)*jdet(jW+1)*w(jW+1);
    of_mk('k<-k+a*x*y',ke,NDN,NDN,coef, ...
     int32([Nnode*ji+Ndof*Nnode*jk ... % block in stiffness matrix
            Nnode*[Nw*(1+jj)+jW Nw*(1+jl)+jW]   ... % columns in NDN
            Ndof Nnode Nnode 1 1]));
   end;end;end;end
catch
   for ji=0:2;for jj=0:2;for jk=0:2;for jl=0:2;
    coef=DD(ji+jj*3+jk*9+jl*27+1)*jdet(jW+1)*w(jW+1);
    in1=Nnode*ji+[1:Nnode];in2=Nnode*jk+[1:Nnode];
    x=NDN(:,Nw*(jj+1)+jW+1);y=NDN(:,Nw*(jl+1)+jW+1);
    ke(in1,in2)=ke(in1,in2)+coef*x*y';
   end;end;end;end
   assignin('caller','ke',ke);
end

   %if norm(k1-k)/norm(k)>eps; error(1);end
end

% Computation of \int_{\Omega} \Sigma : \delta e = \int \Sigma : (F^T.Dpde)
% Sigma_ij F_kj D_ki = r1_ik : (D_ki = N_,i(jW) u_k)
if size(def,2)>1
  F=def(:,2); 
  DofPos=evalin('caller','DofPos(EltConst.VectMap,jElt)');
  DofPos=reshape(double(DofPos)+1,Nnode,Ndof/Nnode);
  r1=Sigma*F_ij';
  for ji=0:2; for jk=0:2;
   coef=r1(ji+3*jk+1)*jdet(jW+1)*w(jW+1);
   fe=NDN(:,Nw*(ji+1)+jW+1);
   F(DofPos(:,jk+1))=F(DofPos(:,jk+1))+coef*fe; % In C direct add to def(:,2)
  end;end
  sp_util('setinput',def,F,int32(1:length(def)),int32(2));

end
