#include <stdio.h>
#include <string.h>
#include "mex.h"
#include "math.h"
#include "../mex/of_def.h"

extern char* pre_cvs();
#define size_3 3	  /* dimension of matrix */

#ifdef FORMATLAB
#include "matrix.h"
#include "../mex/of_def.h"
#include "../mex/of_mk_pre.c"
#include "../mex/hyper.h"
#define of_inc_heart "../mex/of_mk_heart.c"
#ifdef of_inc_heart
#   include "../mex/cardio.h"
#   include "../mex/visco_1.h"
#endif


/* --------------------------------------------------------------*/


#else /* Scilab compile */
#include "stack-c.h"
#include "of_def.h"
#include "of_mk_pre.c"
#endif



void mexFunction (int nlhs, mxArray *plhs[],
                  int nrhs, const mxArray *prhs[]) 
/*-----------------------------------------------------------------------*/
{
  double *out, *constit, *node, *out1, *out2, *estate, *defe, *eltconst, *InfoAtNode;
  char CAM[15], ierr[100]={'\0'};
  int i1, typ;
  int *point, *integ;

  /*-----------------------------------------------------------------------*/
  /* INIT used by all                                                      */

  if (nrhs==0)
    {
      plhs[0] =  mxCreateDoubleMatrix(1,1,mxREAL);
      out = mxGetPr(plhs[0]);
      out[0]=5.0003;
      return;
    }

  i1 = mxGetM(prhs[0])*mxGetN(prhs[0]);
  mxGetString(prhs[0],CAM,15);
  /* CAM = mxCalloc(i1+1,sizeof(char)); mxGetString(prhs[0],CAM,i1+1);*/

  if ((nrhs==1) && (!strcmp("cvs",CAM))) {
    mxArray *st, *st1, *st2;
    char    *ch;
    int     *dims;

    st=mxCreateString("$Revision: 1.86 $  $Date: 2006/05/22 10:53:59 $");
    dims=mxCalloc(2,sizeof(int)); dims[0]=3; dims[1]=2;
    plhs[0]=mxCreateCellArray(2,dims);
    mxSetCell(plhs[0],0,      mxCreateString("of_mk"));
    mxSetCell(plhs[0],1,      mxCreateString("of_mk_subs"));
    mxSetCell(plhs[0],2,      mxCreateString("of_mk_pre"));

    mxSetCell(plhs[0],3,st);
    ch=pre_cvs(); st1=mxCreateString(ch);
    st2=pre_cvs2(); 
    mxSetCell(plhs[0],4,st1);
    mxSetCell(plhs[0],5,st2);
    mxFree(dims); 

    return;
  }

  /* INIT and other non standard calls. TEST IS STRING WITH MORE THAN 8 char
   -------------------------------------------------------------- */
  if (i1>=8){

    if (!strcmp(CAM,"mitcinit")) {
	plhs[0] = mxCreateDoubleMatrix(96,1,mxREAL);
	out = mxGetPr(plhs[0]);
	edemit4_(out);
    /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    } else if (!strcmp("xkx_trans",CAM)) {
    /* k=of_mk('xkx_trans',x,k);  */
    double *offset;
    int    type,ym,yn;

    if(nlhs!=1) mexErrMsgTxt ("One output required.");
    if (mxGetM (prhs[1]) != 3 || mxGetN (prhs[1]) != 3) {
      mexErrMsgTxt ("First argument X must be 3 by 3 in c = xkx(X,K)");
    }
    if (mxIsSparse(prhs[2])) mexErrMsgTxt ("k must be full");
    /* create matrix for the return argument */
    ym=mxGetM(prhs[2]);yn=mxGetN(prhs[2]);
    if (mxIsComplex(prhs[2])) { 
      plhs[0] = mxCreateDoubleMatrix (ym, yn, mxCOMPLEX);
    } else  plhs[0] = mxCreateDoubleMatrix (ym, yn, mxREAL);
    if (nrhs>3) {
     if   (mxGetM(prhs[3])==0) offset=NULL; 
     else if (mxGetM(prhs[3])!=3) { mexErrMsgTxt ("offset must be 3xN");}
     else offset=mxGetPr(prhs[3]); 
    } else offset=NULL;
    if (nrhs==5) type=(int)*mxGetPr(prhs[4]); else type=0;
    if (mxIsComplex (prhs[1])) 
          mexErrMsgTxt ("Matrix of basis change must not be complex.");
    if (mxIsComplex (prhs[2])) { 
      x_k_x(mxGetPr(prhs[1]),mxGetPi(prhs[2]),mxGetPi(plhs[0]),
                                                   offset,type,ym,yn);
    }
    /* call the subroutine xkx */
    x_k_x (mxGetPr (prhs[1]),mxGetPr(prhs[2]),mxGetPr(plhs[0]),
                                                   offset,type,ym,yn);

    /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    }  else if (!strcmp(CAM,"princstress"))  {
        double    *A, *b, c[size_3], DUMMY[1][1], WORK[3*size_3];
        int       j1, ok, c1=size_3, c2=3*size_3, c3=1;
        char      c4='N';
	/*
              of_mk('princstress',eye(3))
 	*/
        if (mxIsSparse(prhs[1])) mexErrMsgTxt("Full matrix needed.");
        A=mxGetPr(prhs[1]);
        if ((mxGetM(prhs[1])!=3) || (mxGetN(prhs[1])!=3)) 
                                 mexErrMsgTxt("Only 3x3 matrix");

        plhs[0]= mxCreateDoubleMatrix(size_3,1,mxREAL); b=mxGetPr(plhs[0]);

	/* of_dgeev(&c4,&c4,&c1,A,&c1,b,c,DUMMY,&c3,DUMMY,&c3, WORK, &c2, &ok);  */
        /*   'N','N',3,input, 3,eigen,VL,  1,   VR,    1  , out, 2*3, out,info)
              SUBROUTINE DGEEV( JOBVL, JOBVR, N, A, LDA, WR, WI, VL, LDVL, VR,
                                $                  LDVR, WORK, LWORK, INFO ) */
	/*
ALTERNATIVE IS:

      SUBROUTINE DSYEV( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, INFO )
*
*  -- LAPACK driver routine (version 3.0) --
*     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
*     Courant Institute, Argonne National Lab, and Rice University
*     June 30, 1999
*
*     .. Scalar Arguments ..
      CHARACTER          JOBZ, UPLO
      INTEGER            INFO, LDA, LWORK, N
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION   A( LDA, * ), W( * ), WORK( * )
*     ..
*
*  Purpose
*  =======
*
*  DSYEV computes all eigenvalues and, optionally, eigenvectors of a
*  real symmetric matrix A.
*
*  Arguments
*  =========
*
*  JOBZ    (input) CHARACTER*1
*          = 'N':  Compute eigenvalues only;
*          = 'V':  Compute eigenvalues and eigenvectors.
*
*  UPLO    (input) CHARACTER*1
*          = 'U':  Upper triangle of A is stored;
*          = 'L':  Lower triangle of A is stored.
*
*  N       (input) INTEGER
*          The order of the matrix A.  N >= 0.
*
*  A       (input/output) DOUBLE PRECISION array, dimension (LDA, N)
*          On entry, the symmetric matrix A.  If UPLO = 'U', the
*          leading N-by-N upper triangular part of A contains the
*          upper triangular part of the matrix A.  If UPLO = 'L',
*          the leading N-by-N lower triangular part of A contains
*          the lower triangular part of the matrix A.
*          On exit, if JOBZ = 'V', then if INFO = 0, A contains the
*          orthonormal eigenvectors of the matrix A.
*          If JOBZ = 'N', then on exit the lower triangle (if UPLO='L')
*          or the upper triangle (if UPLO='U') of A, including the
*          diagonal, is destroyed.
*
*  LDA     (input) INTEGER
*          The leading dimension of the array A.  LDA >= max(1,N).
*
*  W       (output) DOUBLE PRECISION array, dimension (N)
*          If INFO = 0, the eigenvalues in ascending order.
*
*  WORK    (workspace/output) DOUBLE PRECISION array, dimension (LWORK)
*          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
*
*  LWORK   (input) INTEGER
*          The length of the array WORK.  LWORK >= max(1,3*N-1).
*          For optimal efficiency, LWORK >= (NB+2)*N,
*          where NB is the blocksize for DSYTRD returned by ILAENV.
*
*          If LWORK = -1, then a workspace query is assumed; the routine
*          only calculates the optimal size of the WORK array, returns
*          this value as the first entry of the WORK array, and no error
*          message related to LWORK is issued by XERBLA.
*
*  INFO    (output) INTEGER
*          = 0:  successful exit
*          < 0:  if INFO = -i, the i-th argument had an illegal value
*          > 0:  if INFO = i, the algorithm failed to converge; i
*                off-diagonal elements of an intermediate tridiagonal
*                form did not converge to zero.
*
*  =====================================================================
*/

        for (j1=0; j1<size_3; j1++) { 
          if (c[j1]!=0 ) mexErrMsgTxt("Complex eigenvalue exist");
        }
        if (ok!=0) mexErrMsgTxt("Problem to compute eigenvalues");
        /* for (j1=0; j1<size_3; j1++) mexPrintf("%f  %f \n", b[j1],c[j1]);*/


/* -----------------------------------------------------------------------*/
    }  else if (!strcmp(CAM,"buildndn"))  {

        /* rules=of_mk('buildndn',13,integrules,nodeE);
           rules=of_mk('buildndn',2 or 3 or 23,integrules,nodeE);
           rules=of_mk('buildndn',2 or 3 or 23,out,x);               */
        mxArray       *field; 
        double        *nodeE, *NDN, *Nr, *Ns, *Nt, *N, *w, *jdet, *bas, *J;
        int           Nw, Nnode, Nshape, t1, *ti;

        if (mxIsInt32(prhs[1])) { ti = mxGetData(prhs[1]);t1=ti[0];}
        else                    { t1 = (int)(*mxGetPr(prhs[1])); }

        if (!mxIsStruct(prhs[2])) mexErrMsgTxt("EltConst must be a structure");
        if (mxIsEmpty(prhs[2])) mexErrMsgTxt("EltConst must not be empty");
   
        field = mxGetField(prhs[2], 0,"Nr");
        /* number of weight points, shapes*/
        Nw = mxGetM(field); Nshape= mxGetN(field);
        Nr =  mxGetPr (field);
        if ((t1==3) || (t1==2) || (t1==23)) Ns =  mxGetPr (mxGetField(prhs[2], 0,"Ns"));   
        if (t1==3)  Nt =  mxGetPr (mxGetField(prhs[2], 0,"Nt"));   
        N =  mxGetPr (mxGetField(prhs[2], 0,"N"));
        w =  mxGetPr (mxGetField(prhs[2], 0,"w"));
        nodeE = mxGetPr (prhs[3]); Nnode = mxGetM(prhs[3]);
        if ( mxGetM(mxGetField(prhs[2], 0,"jdet"))!=Nw )
          mexErrMsgTxt("Problem of row allocation in jdet");
        if (    mxGetM(mxGetField(prhs[2], 0,"NDN"))
             != mxGetN(mxGetField(prhs[2], 0,"N"))  )
          mexErrMsgTxt("Problem of row allocation in NDN ");

	if (t1==3) {
          if (Nw*4!= mxGetN(mxGetField(prhs[2], 0,"NDN"))) 
            mexErrMsgTxt("Problem of column allocation in NDN");
        } else if (t1==2 || t1==23) {
          if (Nw*3!= mxGetN(mxGetField(prhs[2], 0,"NDN"))) 
            mexErrMsgTxt("Problem of column allocation in NDN");
        } else if (t1==13) {
          if (Nw*3!= mxGetN(mxGetField(prhs[2], 0,"NDN"))) 
            mexErrMsgTxt("Problem of column allocation in NDN");
	} else mexErrMsgTxt("Not a valid NDN mode");
        jdet        = mxGetPr (mxGetField(prhs[2], 0,"jdet"));
        NDN         = mxGetPr (mxGetField(prhs[2], 0,"NDN"));

        if (t1==2)       BuildNDN2d(w,N,Nr,Ns,nodeE,jdet,NDN,Nw,Nnode,Nshape);
        else if (t1==3)  BuildNDN3d(w,N,Nr,Ns,Nt,nodeE,jdet,NDN,Nw,Nnode,Nshape);
        else if (t1==13) BuildNDN1d(w,N,Nr,nodeE,jdet,NDN,Nw,Nnode,Nshape);
        else if (t1==23) {
         field=mxGetField(prhs[2], 0,"bas"); 
         if (field==NULL||mxGetM(field)!=9||mxGetN(field)!=Nw) bas=NULL;
         else  bas = mxGetPr (field);
         field=mxGetField(prhs[2], 0,"J"); 
         if (field==NULL||mxGetM(field)!=4||mxGetN(field)!=Nw) J=NULL; 
         else  J = mxGetPr (field);
         BuildNDNS3d(w,N,Nr,Ns,nodeE,jdet,NDN,Nw,Nnode,Nshape,bas,J);
	}
        else if (t1==13) {

	}
/* ------------------------------------------------------------------------*/
	/* CONDENSE :

 ke=t3p('testmat'); me=ke{2}; ke=ke{1}; [a,b]=elem0('condense',5,ke,me);
 [kr,mr]=of_mk('condense',5,ke,me);
 norm(kr-a)
 norm(mr-b)

*/
}  else if (!strcmp(CAM,"condense"))  {

  int i1, M;
  double *ke, *me, *kr, *mr;

  if (nrhs!=4) mexErrMsgTxt("4 input arguments needed");
  i1=(int)(*mxGetPr(prhs[1]))-1;
  M=mxGetM(prhs[2]);   /* xxx check sizes */
  ke=mxGetPr(prhs[2]);
  me=mxGetPr(prhs[3]); /* xxx add check */

  plhs[0]= mxCreateDoubleMatrix(i1,i1,mxREAL); kr=mxGetPr(plhs[0]);
  if (mxGetM(prhs[3])>0) {
    plhs[1]= mxCreateDoubleMatrix(i1,i1,mxREAL); mr=mxGetPr(plhs[1]);
  } else { mr=NULL;  }

 

  elt_condense(i1, ke, kr, me, mr, M);  


/* ------------------------------------------------------------------------*/
      /* applies the matrix integration rule defined in R1     
        ke=of_mk('matrixintegration',jElt[1],NodePos[2],Case.Node[3], ...
            pointers[4],integ[5],constit[6],gstate[7],elmap[8],
            InfoAtNode[9],EltConst[10],def.def)
           of_mk('matrixintegration',DofPos[1],NodePos[2],Case.Node[3], ...
            pointers[4],integ[5],constit[6],gstate[7],elmap[8],
            InfoAtNode[9],EltConst[10],def.def[11],k[12],opt[13])

      */
}  else if (!strcmp(CAM,"matrixintegrat"))  {

        /*  k=zeros(8); of_mk('matrixintegration',rules,r1,constit,k); */
        mxArray     *field;
        double      *NDN, *jdet, *w, *ke, *kblock, *constit, coef, 
                    *nodeE, *node, *N, *Nr, *Ns, *Nt,
                    *o_pr, *def, *defi, *gstate, *InfoAtNode, *bas, *J,
                    *Ener, *vect, vtmp, *ke_tmp, *Be, *defelem;
        int         jw, j1, j2, j3, j4, *rule_terms, Nterms, Nw, Nshape,Mk, 
                    NdofPerField, Nelt, wstep,i1,i2, jMat, *o_ir, *o_jc, *opt,
                    *elmap, *DofPos, DofPerElt, StrategyType, *VectMap;
        int         jElt, *NodePos, Nnode, Nmnode,*point,*integ,
                    ind_ts_eg[9]={0,5,4,5,1,3,4,3,2}, *topo, Ntopo, Nstate,
                    ci_ts_eg[6]={0,4,8,7,6,3}, Mdef, Ndef, opt1[5],
		    *cEGI, Mener, MInfoAtNode, Ncondense,
                    *rrule, Nrule, Mrule, *ind, Nk, NDNtype=0;
	/* constitutive energy function handle (Integ,Constit,I,dWdI,d2WdI2) */

#define of_inc_fs "../mex/of_mk_fsmatrix.c"
#define MatrixIntegrationStep 0
#include of_inc_fs
#include of_inc_heart
#ifdef of_inc_user
#   include of_inc_user
#endif

     /* inputs */
     NodePos=mxGetData(prhs[2]); 
     if (!mxIsInt32(prhs[2])) mexErrMsgTxt("NodePos must be int32");

     Nnode=mxGetM(prhs[2]); nodeE=mxCalloc(Nnode*4,sizeof(double));
     Nmnode=mxGetM(prhs[3]); node=mxGetPr(prhs[3]);

     if (!mxIsStruct(prhs[10])) mexErrMsgTxt("EltConst must be a structure");
     if (mxIsEmpty(prhs[10])) mexErrMsgTxt("EltConst must not be empty");

     field = mxGetField(prhs[10], 0,"Ncondense");
     if (field!=NULL) {
       Ncondense=(int)*mxGetPr (field);
       ke_tmp=mxGetPr(mxGetField(prhs[10], 0,"CondenseMat"));
       ind=(int*)mxGetData(mxGetField(prhs[10], 0,"CondenseInd"));
       Nk=mxGetM(mxGetField(prhs[10], 0,"CondenseMat"));
     } else {Ncondense=0;}

     field = mxGetField(prhs[10], 0,"NDN");
     if (field!=NULL) {
          NDN = mxGetPr (field);NdofPerField=mxGetM(field);
          jdet    = mxGetPr (mxGetField(prhs[10], 0,"jdet"));
          field   = mxGetField(prhs[10], 0,"w");
          w  = mxGetPr (field); Nw=mxGetM(field); w = w+3*Nw;
          Nshape=mxGetN(mxGetField(prhs[10], 0,"Nr"));
     }
     N=mxGetPr(mxGetField(prhs[10], 0,"N"));
     Nr=mxGetPr(mxGetField(prhs[10], 0,"Nr"));
     Ns=mxGetPr(mxGetField(prhs[10], 0,"Ns"));
     field = mxGetField(prhs[10], 0,"Nt");
     if (field!=NULL) Nt=mxGetPr(field); else Nt=NULL;
     field=mxGetField(prhs[10], 0,"bas"); 
     if (field==NULL||mxGetM(field)!=9||mxGetN(field)!=Nw) bas=NULL; 
     else  bas = mxGetPr (field);
     field=mxGetField(prhs[10], 0,"J"); 
     if (field==NULL||mxGetM(field)!=4||mxGetN(field)!=Nw) J=NULL; 
     else  J = mxGetPr (field);
     InfoAtNode=mxGetPr(prhs[9]); MInfoAtNode=mxGetM(prhs[9]);
     field=mxGetField(prhs[10], 0,"VectMap");
     if (field!=NULL) {VectMap=(int*)mxGetData(field);}
     field=mxGetField(prhs[10], 0,"RhsDefinition");
     if (field!=NULL) { 
       rrule=(int*)mxGetData(field);Mrule=mxGetM(field);Nrule=mxGetN(field);
     }

     if (nrhs>12) { /* direct assembly of a complete group  - - - - - - - - */

       Nelt=mxGetN(prhs[1]); jElt=0;
       o_pr=mxGetPr(prhs[12]); o_ir=mxGetIr(prhs[12]); o_jc=mxGetJc(prhs[12]); 
       opt=mxGetData(prhs[13]); DofPerElt=opt[0]; 
       elmap=mxGetData(prhs[8]);DofPos=mxGetData(prhs[1]);
       if (opt[1]==-1) { /* xxx ENER GetField(cEGI), GetField(Ener) */
         cEGI=(int*)mxGetData(mxGetField(prhs[12], 0,"cEGI"));
         Ener=mxGetPr(mxGetField(prhs[12], 0,"Ener"));
         Mener=mxGetM(mxGetField(prhs[12], 0,"Ener"));
         rrule=NULL;
       }
       def=mxGetPr(prhs[11]);defi=mxGetPi(prhs[11]);
       Mdef=mxGetM(prhs[11]);Ndef=mxGetN(prhs[11]);

     } else { /* assembly of a single element matrix */
       jElt=(int)*mxGetPr(prhs[1])-1; Nelt=jElt+1;
       opt=NULL; rrule=NULL;
     } /*  assembly of a single matrix or a full group - - - - - - - - - - - */

     /* initialize element output matrices - - - - - - - - - - - - - - */

     point   = (int*)mxGetData(prhs[4]);
     integ   = (int*)mxGetData(prhs[5])+point[5];
     Mk=integ[2]+Ncondense; 

     if (opt==NULL && point[4]==0) { /* two matrices M,K return as outputs */
       plhs[0] =  mxCreateDoubleMatrix(Mk,Mk,mxREAL);
       plhs[1] =  mxCreateDoubleMatrix(Mk,Mk,mxREAL);
       if (nlhs==1) mexErrMsgTxt("m and k should be returned for MatDes 0");
     } else if (opt==NULL) { /* one matrix returned or assembled */ 
       plhs[0] =  mxCreateDoubleMatrix(Mk,Mk,mxREAL);
       plhs[1] =  mxCreateDoubleMatrix(0,0,mxREAL);
     } else if (opt[1]==-1) { /* energy computations assembled */
       vect=(double*)mxCalloc(DofPerElt,sizeof(double));  
       plhs[0] =  mxCreateDoubleMatrix(Mk,Mk,mxREAL);
     } else { /* one matrix returned or assembled */ 
       plhs[0] =  mxCreateDoubleMatrix(Mk,Mk,mxREAL);
       plhs[1] =  mxCreateDoubleMatrix(0,0,mxREAL);
     }

 jMat=0; StrategyType=0;
 while (jMat<3) { /* matrices to generate allow for multiple */

   point   = (int*)mxGetData(prhs[4]);
   if      (point[4]==0 &&jMat==0) { /* mass */
       point[4]=2;jMat=2; ke=mxGetPr(plhs[1]); Mk=mxGetM(plhs[1]);
   } else if (point[4]==2 && jMat==2) { /* stiffness assembled in mass+stiff*/
       point[4]=1; ke=mxGetPr(plhs[0]); Mk=mxGetM(plhs[0]);
       jElt=(int)*mxGetPr(prhs[1])-1;
   } else if (point[4]==1 &&jMat==2) { /* exit when point[4]==0 and m and k done */ 
      point[4]=0;jMat=3;break; /*exit*/
   } else if (point[4]==3  &&jMat==0) {   /* viscous damping */
       jMat=3; ke=mxGetPr(plhs[0]); Mk=mxGetM(plhs[1]);
       for (j1=0;j1<Mk*Mk;j1++) ke[j1]=0.;
       break; /* xxx do something if defined */
   } else { /* standard case one generates one element matrix */
       ke=mxGetPr(plhs[0]); Mk=mxGetM(plhs[0]); jMat=3; 
   }
   /* prepare for different matrix assembly strategies - - - - - - - - - - */
   field = mxGetField(prhs[10],0,"material");
   if (field==NULL || point[4]<4) { /* generic multiphysic linear element */
    field = mxGetField(prhs[10],0,"MatrixIntegrationRule");
    if ((field==NULL)||(mxGetM(field)==0)) {
      field = mxGetField(prhs[10],0,"material");
      if (field==NULL) mexErrMsgTxt("Not a supported MatrixIntegrationRule 1");
    } else if (point[4]-1>mxGetN(field)) {
      field = mxGetField(prhs[10],0,"material");
      if (field==NULL) mexErrMsgTxt("Not a supported MatrixIntegrationRule 2");
    } else {
     field=mxGetCell(field,point[4]-1); /* rule for given matrix type */
     rule_terms   = mxGetData(field); Nterms=mxGetM(field);
     StrategyType=1;field=NULL;
    }
   }
   if (field!=NULL) { 
     mxGetString(field,CAM,15);
     if (!strcmp("Elastic3DNL",CAM)) { /* geometric non linear elastic 3D */
       StrategyType=2;
       field = mxGetCell(mxGetField(prhs[10],0,"ConstitTopology"),0);
       if (!mxIsInt32(field)) mexErrMsgTxt("ConstitTopology must be int32");
       topo=mxGetData(field);Ntopo=mxGetM(field);
     } else if (!strcmp("hyperelastic",CAM)) { /* hyperelastic */
       StrategyType=3;
       NDNtype=1;
       field = mxGetCell(mxGetField(prhs[10],0,"ConstitTopology"),0);
       if (!mxIsInt32(field)) mexErrMsgTxt("ConstitTopology must be int32");
       topo=mxGetData(field);Ntopo=mxGetM(field);
       Be=mxCalloc(3*Nnode,sizeof(double));
       defelem=mxCalloc(3*Nnode*Ndef,sizeof(double));
     } 
#undef MatrixIntegrationStep 
#define MatrixIntegrationStep 1
#include of_inc_fs
#include of_inc_heart
#ifdef of_inc_user
#   include of_inc_user
#endif
     else {
       mexPrintf("'%s' not supported",CAM);mexErrMsgTxt(" not supported");
     } 
     gstate=mxGetPr(prhs[7]);Nstate=mxGetM(prhs[7]);
   }
for (;jElt<Nelt;jElt++) { /* loop on elements to assemble */
   int *CurDofPos;

     CurDofPos=DofPos+jElt*DofPerElt;
     point   = (int*)mxGetData(prhs[4])+jElt*mxGetM(prhs[4]);
     point[4]= *((int*)mxGetData(prhs[4])+4);

     integ   = (int*)mxGetData(prhs[5])+point[5];
     constit =       mxGetPr(prhs[6])+point[6];
     for (j2=0;j2<Nnode;j2++) {/* nodeE nodes of current element */
         i1=NodePos[j2+Nnode*jElt]-1;
         nodeE[j2]=node[i1+Nmnode*4]; 
         nodeE[j2+Nnode]=node[i1+Nmnode*5]; 
         nodeE[j2+2*Nnode]=node[i1+Nmnode*6]; 
     }

     if (point[3]==3) { /* build the jacobian and [N N,x N,y N,z] */
       if (NDNtype == 1) {
	 BuildNDN3dW(w,N,Nr,Ns,Nt,nodeE,jdet,NDN,Nw,Nnode,Nshape);/*Nw,Nnode,4*/
       }
       else {
	 BuildNDN3d(w,N,Nr,Ns,Nt,nodeE,jdet,NDN,Nw,Nnode,Nshape);/*Nnode,Nw,4*/
       }
     } else if (point[3]==2) {
        BuildNDN2d(w,N,Nr,Ns,nodeE,jdet,NDN,Nw,Nnode,Nshape);
     } else if (point[3]==23) {
         BuildNDNS3d(w,N,Nr,Ns,nodeE,jdet,NDN,Nw,Nnode,Nshape,bas,J);
     } else {mexErrMsgTxt("Not a supported NDN method");}

#undef MatrixIntegrationStep 
#define MatrixIntegrationStep 2

     switch (StrategyType) {
     case 1: {
       /*  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
        GENERIC LINEAR MULTIPHYSIC ELEMENTS
        standard mat_og matrix assembly with MatrixIntegrationRule
      */
        double *RHS, r1;
        for (j2=0;j2<Nterms;j2++) { /* loop on terms of integration rule */ 
	 
         wstep=rule_terms[j2+Nterms*5]; 
         for (jw=0;jw<rule_terms[j2+Nterms*6];jw++) {/* loop on integ. points*/ 

	  i1=rule_terms[j2+4*Nterms];
          if (i1<0) {coef=-1;} else coef=1;
	  i1=abs(i1); i2=jw+rule_terms[j2+Nterms*7];
          coef = coef* w[i2]*jdet[i2] * constit[ jw*wstep+i1 ];
          if (coef!=0) { 
	  /*    mexPrintf("consit(%i)*%.2f(%i)*NDN(i,%i)\n", 
	       jw*wstep+abs(i1)+1,w[i2],i2,(1+jw+rule_terms[j2+2*Nterms]));  */
          for (j4=0;j4<NdofPerField;j4++) { 
	    kblock = ke + (rule_terms[j2+Nterms]+j4)*Mk + rule_terms[j2]; 
            for (j3=0;j3<NdofPerField;j3++) {
             *kblock += coef*NDN[j3+(jw+rule_terms[j2+2*Nterms])*NdofPerField] *
                             NDN[j4+(jw+rule_terms[j2+3*Nterms])*NdofPerField ];
             kblock++; 
	    } /* j3 */
	   } /* j4 */

          } /* coef!= 0*/
         } /* point */
        } /* terms of constitutive relation */

        if (rrule!=NULL && (Ndef>1) && (point[4]==1 | point[4]==5) ) { 
         int Ndof, jj, jk, jrule;
/*
if size(def,2)>1&isfield(EltConst,'RhsDefinition')
  rrule=double(EltConst.RhsDefinition);%[101 0 1 0 0 0 -1 0 8]
  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
*/

         RHS=def+Mdef; /* def.def(:,2) */
         Ndof=(int)(MInfoAtNode/Nnode); /* 1 = 8/8*/

         for (jrule=0;jrule<Mrule;jrule++) {
          switch (rrule[jrule]) {
            case 101 :  /* Volume load defined in InfoAtNode */
            {
	      if (!(mxGetN(prhs[9])*MInfoAtNode)) 
                       mexErrMsgTxt("InfoAtNode is empty");

             for (jw=rrule[7+Mrule*jrule];
                  jw< rrule[7+Mrule*jrule]+rrule[8+Mrule*jrule];jw++) { 
              r1=0.;
	      /* r1=EltConst.NDN(:,jw+r3(4)+1)'
                    *InfoAtNode(  r3(2)+[0:Nnode-1]*r3(3)+1  ,jElt);  */
              for (jk=0;jk<Nnode;jk++) {
                r1 += NDN[jk+Nnode*(jw+rrule[3+Mrule*jrule])] * 
                      InfoAtNode[ rrule[1+Mrule*jrule]+jk*rrule[2+Mrule*jrule]
                                 +MInfoAtNode*jElt];                
	      }
              r1*=jdet[jw]*w[jw]; /* r1=r1*jdet(jw+1)*w(jw+1); */
	      /* if r3(7)>=0; r1=r1*EltConst.bas(r3(7)+1,jw+1);end  */
              if (rrule[6+Mrule*jrule]>=0) mexErrMsgTxt(".bas not implemented"); 
	      /* in1=double(DofPos(EltConst.VectMap(Nnode*r3(5)+[1:Nnode]),jElt));
		 F(in1)=F(in1)+  EltConst.NDN(:,jw+r3(6)+1)*r1; */
              for (jk=0;jk<Nnode;jk++) { 
		RHS[CurDofPos[ VectMap[ Nnode*(rrule[4+Mrule*jrule])+jk ]-1 ]]
                   += NDN[ jk+Nnode*(jw+rrule[5+Mrule*jrule]) ]*r1;
	      }      

	     } /* jw */
	    }
	    default : 
            { /* mexErrMsgTxt("Not a supported RHS method"); */ }
	  } /* switch */
	 } /* jrule */
	}
        break;
     }
     case 3:  
     { /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     geometric non linear hyperelastic 3D
     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
       double e0[6], e[6], F_ij[9], Sigma[6], d2wde2[36],zero=0;
       double d2I3dcdc[36], dI2dc[6], dI3dc[6], I[3], dWdI[3], d2WdI2[9], *RHS; 
      
       int *CurDofPos,nb,nul=0,unit=1,jk,jj,j1,j2;

       CurDofPos=DofPos+jElt*DofPerElt;

       nb=3*Nnode;
       of_dcopy(&nb,&zero,&nul,Be,&unit);

       for (j1=0;j1<nb;j1++) { for (j2=0;j2<Ndef;j2++) {
	   defelem[j1+nb*j2]=def[CurDofPos[j1]+Mdef*j2];
	 }}

 
       for (jw=0;jw<Nw;jw++) {/* loop on integ. points --------------------- */

	 for (j1=0;j1<6;j1++) { 
	   for (j2=0;j2<6;j2++) {
	     d2wde2[j1+6*j2]=0;
	   }	   
	   Sigma[j1]=0;
	 }
	 
	
	 elemcalc(jw,Nw,Nnode,defelem,NDN,e0,e,F_ij,I,dI2dc,dI3dc,d2I3dcdc,nb);
	
        
	 EnPassiv(integ,constit,I,dWdI,d2WdI2);

	 matPassiv(opt1,Sigma,dWdI,d2WdI2,dI2dc,dI3dc,d2I3dcdc,d2wde2);

	 Mecha3DIntegH(ke,Be,F_ij,d2wde2,Sigma,w,jdet,NDN,Nnode,nb,Nw,jw);

	
       }

       RHS=def+Mdef; /* def.def(:,2) */
       for (jk=0;jk<3;jk++) {
	 for (jj=0;jj<Nnode;jj++) {
	   RHS[CurDofPos[jk+3*jj]]+=Be[Nnode*jk+jj];
	 }
       }
     }
     break;

     case 2: 
     { /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
       geometric non linear elastic 3D (see elem0.m for m file implementation)
      defe=def(double(DofPos(:,jElt))+1,1); 
      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(:);
      % assignin('base','d2wde2_elas',d2wde2);
     %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)
      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -  */

      double U_ij[9], F_ij[9], d2wde2[36], Sigma[9];
      int ji,jj,jk;
      for (ji=0;ji<36;ji++) { 
       if (topo[ji]==0) d2wde2[ji]=0; else; d2wde2[ji]=constit[topo[ji]-1];
      }
      /*	for (ji=0;ji<36;ji+=3)  { 
	  mexPrintf("(%i) %.2f (%i) %.2f (%i) %.2f\n",ji,d2wde2[ji],
           ji+1,d2wde2[ji+1],ji+2,d2wde2[ji+2]);
	   }*/

      for (jw=0;jw<Nw;jw++) {/* loop on integ. points*/
        /* build U and F */
        for (ji=0;ji<9;ji++)  U_ij[ji]=0;
        for (ji=0;ji<3;ji++){ for (jj=0;jj<3;jj++) {for (jk=0;jk<Nnode;jk++) {
         U_ij[ji+3*jj]+=def[CurDofPos[ji+3*jk]]*NDN[jk+Nnode*(Nw*(jj+1)+jw)];
        }}}
	for (ji=0;ji<9;ji++)  {   F_ij[ji]=U_ij[ji];}
        F_ij[0]+=1;F_ij[4]+=1;F_ij[8]+=1;
        /* compute Sigma */
        for (ji=0;ji<9;ji++)  Sigma[ji]=0;
        for (ji=0;ji<9;ji++) { for (jj=0;jj<9;jj++) {
	  Sigma[ji]+=d2wde2[ind_ts_eg[ji]+6*ind_ts_eg[jj]]*U_ij[jj];
        }}
        if (gstate!=NULL) { /* allow for non elastic stress */
          jj=Nstate*jElt+jw*6; for (ji=0;ji<9;ji++) {   
	    Sigma[ji]+=gstate[jj+ind_ts_eg[ji]];
          }
        } /* allow for non elastic stress */
        Mecha3DInteg(ke,F_ij,d2wde2,Sigma,w,jdet,NDN,Nnode,3*Nnode,Nw,jw);
      }
        break;
#include of_inc_fs
#include of_inc_heart
#ifdef   of_inc_user
#   include of_inc_user
#endif

     case 999:   
       /* When introducing a new element family you need to 
        use standard variables : integ,constit,nodeE,NDN
        to fill in
         ke : the element matrix (which will be assemble below)
        and possibly a current RHS as illustrated in the Hyperelastic
        strategy. You will then need to use the DofPositions of the current
        element
         CurDofPos=DofPos+jElt*DofPerElt;
        and add the elemental RHS to the current force in def.def(:,2) 
         F=def+Mdef; 
       */
       mexErrMsgTxt( "User element family place holder" );

      } default:
       mexPrintf("jMat %i (strategy %i)",jMat,StrategyType);
       mexErrMsgTxt( "Not a supported matrix strategy" );
    
     } /* standard mat_og matrix assembly */

     if (Ncondense) {
       double alpha;
       int j3;
       /* xxx JM renumeroter la matrice, condenser les DDLs
          voir elem0.m ligne 402
          pour finaliser les inits cherche Ncondense ci-dessus
          Ne faire aucune alloc.
       */
       opt1[2]=1; opt1[3]=1; /* no increment */
       for (j1=0; j1<Nk; j1++) { for (j2=0; j2<Nk; j2++) {
        ke_tmp[j1+Nk*j2]=ke[ind[j1]+Nk*ind[j2]];
       }} 
       /* condense step by step */
       for (j1=0; j1<Ncondense; j1++) {
         opt1[0]=Nk-j1-1; opt1[1]=Nk-j1-1; opt1[4]=Nk;
         alpha=-1./ke_tmp[(Nk-j1-1)*Nk+Nk-j1-1];
         of_dger(&opt1[0], &opt1[1],  &alpha,
                 &ke_tmp[(Nk-j1-1)*Nk], &opt1[2],
                 &ke_tmp[(Nk-j1-1)*Nk], &opt1[3], /* valid if  symmetric */
                 &ke_tmp[0],            &opt1[4]);
	 /*  mexPrintf("%i %i %i %i %i  \n",
	     opt1[0],opt1[1],(Nk-j1-1)*Nk+Nk-j1-1,(Nk-j1-1)*Nk,(Nk-j1-1)*Nk); */
       }
       /*  Par rapport a ton condense_matrix remettre 
            une phase de renumerotation de ke_tmp -> ke  */
       j3=0;
       for (j1=0; j1<Nk-Ncondense; j1++) { for (j2=0; j2<Nk-Ncondense; j2++) {
        ke[j3]=ke_tmp[j1+Nk*j2];
        j3++;
       }} 

     }
    
     if ((opt!=NULL) && (opt[1]!=-1)){ 
       /* assemble current element matrix if appropriate */
       AssembleSparse(o_ir,o_jc,o_pr,opt[0],opt[1],opt[2],
         CurDofPos,elmap,NULL,NULL,ke);
       for (j2=0;j2<Mk*Mk;j2++) ke[j2]=0;
     } else if (opt==NULL) {}
     else {/* energy output  - - - - - - - - - - - - - - - - - - -*/
        for (j2=0;j2<Ndef;j2++) {
          for (j1=0;j1<DofPerElt;j1++) {
            i1 = CurDofPos[VectMap[j1]-1];
            if (i1<0) vect[j1]=0.; else { vect[j1]=def[i1+j2*Mdef]; }
          }
          vtmp=0.;
          for (j1=0;j1<DofPerElt;j1++) { for(j3=0;j3<DofPerElt;j3++) {
            vtmp+=vect[j3]*ke[j3+j1*DofPerElt]*vect[j1];
          }}
          if (defi!=NULL) {
           for (j1=0;j1<DofPerElt;j1++) {
            i1 = CurDofPos[VectMap[j1]-1];
            if (i1<0) vect[j1]=0.; else { vect[j1]=defi[i1+j2*Mdef]; }
           }
           for (j1=0;j1<DofPerElt;j1++) { for(j3=0;j3<DofPerElt;j3++) {
            vtmp+=vect[j3]*ke[j3+j1*DofPerElt]*vect[j1];
           }}
          } /* if imaginary part consider it */
          Ener[cEGI[jElt]-1 + j2*Mener]=vtmp;
        } /* j2 */
       for (j2=0;j2<Mk*Mk;j2++) ke[j2]=0.;
     } /* (opt!=NULL) - - - - - - - - - - - - - - - - - - - - - - - */



   } /* end of jElt loop on elements */ 
 } /* end multiple matrix case */

/* strategy dependent clean up */
#undef MatrixIntegrationStep 
#define MatrixIntegrationStep 3
switch (StrategyType) {
 case 3: {
   if (Be!=NULL) mxFree(Be); 
   if (defelem!=NULL) mxFree(defelem);
   break;
 }
#include of_inc_fs
#include of_inc_heart
#ifdef   of_inc_user
#   include of_inc_user
#endif
}

if (opt==NULL) {} else if (opt[1]==-1) { mxFree(vect); }

mxFree(nodeE);


/* ------------------------------------------------------------------------*/
} else if (!strcmp(CAM,"StressObserve"))  {

	/*
hexa8b('testeigstress');

 nw=EltConst.Nw; Nnode=EltConst.Nnode;
 if point(4)==3;of_mk('buildndn3d',EltConst,nodeE);Ndof=3;ke=zeros(6*nw,Nnode*3);
 elseif point(4)==2;of_mk('buildndn2d',EltConst,nodeE);Ndof=2;ke=zeros(3*nw,Nnode*2);
 end
 for jConst=1:size(rule,1) % loop on elements of rule
  coef=constit(rule(jConst,3)+1);
  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

 
 ke=of_mk('StressObserve',EltConst,int32(rule),constit,nodeE,point);


	*/

        mxArray       *field; 
        double        *nodeE, *NDN, *Nr, *Ns, *Nt, *N, *w, *jdet, *ke,
                      *constit, coef;
        int           Nw, Nnode, Nshape, *point, *rule, nNDN,
                      jConst, jN, jW, Mrule, i1, i2, Ndof, Nstress;

        /* inputs */
        if (nrhs!=6) { mexErrMsgTxt("6 input arguments needed."); }
        field = mxGetField(prhs[1], 0,"Nr");
        Nw = mxGetM(field); Nshape= mxGetN(field);/* number of weight points, nodes*/
        Nr =  mxGetPr (field);
        Ns =  mxGetPr (mxGetField(prhs[1], 0,"Ns"));   
        field=mxGetField(prhs[1], 0,"Nt");
        if (field !=NULL) Nt =  mxGetPr(field); else Nt=NULL;
        N =  mxGetPr(mxGetField(prhs[1], 0,"N"));
        w =  mxGetPr (mxGetField(prhs[1], 0,"w"));

        jdet        = mxGetPr (mxGetField(prhs[1], 0,"jdet"));
        NDN         = mxGetPr (mxGetField(prhs[1], 0,"NDN"));
        nNDN        = mxGetM (mxGetField(prhs[1], 0,"NDN"));

        nodeE = mxGetPr (prhs[4]); Nnode = mxGetM(prhs[4]);

        if (!mxIsInt32(prhs[5])) { mexErrMsgTxt("point must be int32"); }
        point =  mxGetData(prhs[5]); 
        rule =  mxGetData(prhs[2]); 
        Mrule = mxGetM(prhs[2]); /* size(rule,1) */
        Nstress=rule[Mrule-1]+1;
        if (point[3]==3) {
           BuildNDN3d(w,N,Nr,Ns,Nt,nodeE,jdet,NDN,Nw,Nnode,Nshape);
           plhs[0] = mxCreateDoubleMatrix(Nstress*Nw,Nnode*3,mxREAL); 
           ke = mxGetPr(plhs[0]); Ndof=3;
           for (jConst=0;jConst<Nstress*Nw*Nnode*3;jConst++) ke[jConst]=0.;
	} else if (point[3]==2) {
           BuildNDN2d(w,N,Nr,Ns,nodeE,jdet,NDN,Nw,Nnode,Nshape); 
           plhs[0] = mxCreateDoubleMatrix(Nstress*Nw,Nnode*2,mxREAL); 
           ke = mxGetPr(plhs[0]); Ndof=2;
           for (jConst=0;jConst<Nstress*Nw*Nnode;jConst++) ke[jConst]=0.;
	} else mexErrMsgTxt("Unknown element type : 2D or 3D?");

        constit = mxGetPr(prhs[3]);

        for (jConst=0;jConst<Mrule ;jConst++) {
          coef=constit[rule[jConst+2*Mrule]+point[6]];
          for (jN=0;jN<Nnode ;jN++) {
            for (jW=0;jW<Nw ;jW++) {
             i1=jW*rule[jConst+4*Mrule]+rule[jConst];
             i2=rule[jConst+3*Mrule]+jN*Ndof;
	     ke[i1+Nw*Nstress*i2] += NDN[jN+ nNDN * (rule[jConst+Mrule]+jW) ]*coef; 
	    } /* jW */
	  } /* jN */
	} /* jConst */

/* ------------------------------------------------------------------------*/
	/*
addpath c:/balmes/sdt.cur/6.5
ofutil of_mk
k=zeros(10); x=[1 2 3]; of_mk('k<-k+a*x*y',k,x,x,1,int32([0 0 0 size(k,1)  3 3 1 1]));k
k=zeros(10); x=[1 2 3;4 5 6]'; 
of_mk('k<-k+a*x*y',k,x,x,1,int32([0 0 0 size(k,1)  3 3 1 1]));k

	 */
} else if (!strcmp(CAM,"k<-k+a*x*y"))  {

   double        *a, *x, *y, *k;
   int          *opt;

   opt=mxGetData(prhs[5]);
   k=mxGetPr(prhs[1])+opt[0];  x=mxGetPr(prhs[2])+opt[1];  
   y=mxGetPr(prhs[3])+opt[2];  a=mxGetPr(prhs[4]);
   /* mexPrintf("a=%g %i %i",a[0],opt[6],opt[7]); */
   /* opt=[offk[0] offx[1] offy[2] size(k,1)[3] m[4] n[5]  incx[6] incy[7]] */
   if (mxGetM(prhs[4])>1) {
     int ji,jj;
     for (ji=0;ji<mxGetM(prhs[4]);ji++) {for (jj=0;jj<mxGetN(prhs[4]);jj++) {
       if (a[0]!=0) {
        of_dger(&opt[4],&opt[5],a,
          x+opt[4]*ji,&opt[6],
          y+opt[5]*jj,&opt[7],k,&opt[3]);
       }
       a++;
     }}
   } else { of_dger(&opt[4],&opt[5],a,x,&opt[6],y,&opt[7],k,&opt[3]);}


/* ------------------------------------------------------------------------*/
} else if (!strcmp(CAM,"Mecha3DInteg"))  {

  int* opt;
  opt=mxGetData(prhs[8]); 

 /* of_mk('Mecha3DInteg',k,F,d2wde2,Sigma,w,jdet,NDN,int32([Nnode Ndof Nw jW])); */
 Mecha3DInteg(mxGetPr(prhs[1]),mxGetPr(prhs[2]),mxGetPr(prhs[3]),
              mxGetPr(prhs[4]),mxGetPr(prhs[5]),mxGetPr(prhs[6]),mxGetPr(prhs[7]),
              opt[0],opt[1],opt[2],opt[3]);

/*---------------------------------------------------------  */
/* insert a matrix in a sparse one                           */
/*---------------------------------------------------------  */
}  else if (!strcmp(CAM,"asmsparse"))  {

      /*
k1=k+0.;
      sp_util('asmsparse',k,int32(i1),ke,elmap);
      of_mk('asmsparse',k1,int32(i1),ke,elmap);
size(find(k-k1))

       */      
    double      *pr,  *val, *mv; 
    int         nVal, ie, lr, lc, gr, gc, ii, *keind, *elmap, jj, *vir, *vjc,
                *ir, *jc, NDDL, IsSymK, IsSymVal, *opt; 
    mxClassID   typeK1;

    if (nrhs < 4 ) { mexErrMsgTxt( "At least 4 args required!" );  }

    ir = mxGetIr( prhs[1] ); jc = mxGetJc( prhs[1] ); pr = mxGetPr( prhs[1] ); 
    opt = mxGetData( prhs[2] ); keind = opt + 3;
    NDDL = opt[0]; IsSymVal = opt[1]; IsSymK = opt[2];  /* get options */
    nVal = mxGetNumberOfElements( prhs[3] );          /* length of values */

    if ((IsSymVal && (nVal != (NDDL*(NDDL+1)/2)))) {
      mexPrintf("%i %i %i",nVal,NDDL,IsSymVal);
      mexErrMsgTxt( "Number of values does not match triangle!" );
    }
    if (((!IsSymVal) && (nVal != (NDDL*NDDL)))) {
      mexPrintf("%i %i %i",nVal,NDDL,IsSymVal);
      mexErrMsgTxt( "Number of values does not match full matrix!" );
    }
    typeK1=mxGetClassID(prhs[3]);
    switch (typeK1) {
    case mxDOUBLE_CLASS:
      vir=NULL; vjc=NULL; break;
    case mxSPARSE_CLASS:
      vir = mxGetIr( prhs[3] ); vjc = mxGetJc( prhs[3] ); break;
    default:
      mexErrMsgTxt( "4th arg must be a sparse or double matrix!" );
    } 
    
 
    if (nrhs<5) {elmap=NULL;} else {elmap = mxGetData( prhs[4] );}

    AssembleSparse( mxGetIr(prhs[1]), mxGetJc(prhs[1]), mxGetPr(prhs[1]), 
                    opt[0],opt[1],opt[2],opt+3,elmap,
		    vir,vjc,mxGetPr(prhs[3]) );

    /* other input checking should be here !!! */

    /* Build output args. */
#ifdef FORMATLAB
    if (AllowMissing == 1) {
      if (nMiss > 0) {
	plhs[0] = mxCreateDoubleMatrix( 3, nMiss, mxREAL ); val = mxGetPr( plhs[0] );
	for (ii = 0; ii < nMiss; ii++) {
	  val[3*ii+0] = mij[2*ii+0];
	  val[3*ii+1] = mij[2*ii+1];
	  val[3*ii+2] = mv[ii];
	}
      } else { plhs[0] = mxCreateDoubleMatrix( 0, 0, mxREAL );  }
      nMiss = 0;
    }
#else
    plhs[0] = prhs[1];
    if (AllowMissing == 1) {
      if (nMiss > 0) {
	plhs[1] = mxCreateDoubleMatrix( 3, nMiss, mxREAL ); val = mxGetPr( plhs[0] );
	for (ii = 0; ii < nMiss; ii++) {
	  val[3*ii+0] = mij[2*ii+0];
	  val[3*ii+1] = mij[2*ii+1];
	  val[3*ii+2] = mv[ii];
	}
      } else { plhs[1] = mxCreateDoubleMatrix( 0, 0, mxREAL ); }
      nMiss = 0;
    }
#endif

/* initmissing -----------------------------------------------------------
     sp_mk( 'initmissing', type, size );
 ------------------------------------------------------------------ */
}  else if (!strcmp(CAM,"initmissing"))  {

    AllowMissing = (int)mxGetScalar( prhs[1] );
    nMissAlloc = (int)mxGetScalar( prhs[2] );

    mexPrintf( "missing: mode: %d, alloc: %d\n", AllowMissing, nMissAlloc );

    nMiss = 0;
    mij = mxCalloc( 2 * nMissAlloc, sizeof( int ) );
    mv = mxCalloc( nMissAlloc, sizeof( double ) );
    mexLock();
    mexMakeMemoryPersistent( mij );
    mexMakeMemoryPersistent( mv );
  /* getmissing -----------------------------------------------------------
     mtx = sp_mk( 'getmissing' );
     ------------------------------------------------------------------ */
}  else if (!strcmp(CAM,"getmissing"))  {
    int ii;
    double *val;

    mexPrintf( "miss: %d\n", nMiss );
    if (AllowMissing == 2) {
      if (nMiss > 0) {
	plhs[0] = mxCreateDoubleMatrix( 3, nMiss, mxREAL );
	val = mxGetPr( plhs[0] ); 
	for (ii = 0; ii < nMiss; ii++) {
	  val[3*ii+0] = mij[2*ii+0];
	  val[3*ii+1] = mij[2*ii+1];
	  val[3*ii+2] = mv[ii];
	}
      } else {
	plhs[0] = mxCreateDoubleMatrix( 0, 0, mxREAL );
      }
    } else {
      plhs[0] = mxCreateDoubleMatrix( 0, 0, mxREAL );
    }
    nMiss = nMissAlloc = 0;
    mxFree( mij );
    mxFree( mv );
    mexUnlock();

  /* fillval -----------------------------------------------------------

     set all values of a (sparse) matrix or a vector to a constant

     k=sparse(ones(10,10));of_mk('fillval',k,0);

     ------------------------------------------------------------------ */
}  else if (!strcmp(CAM,"fillvalue"))  {

    int nVal;
    double *pr, *pr0, val;

    pr0 = pr = mxGetPr( prhs[1] ); 
    if (mxIsSparse( prhs[1] )) { nVal = mxGetNzmax( prhs[1] );
    } else if (mxIsDouble( prhs[1] )) { nVal = mxGetNumberOfElements( prhs[1] );
    } else {
      mexErrMsgTxt( "Unsupported input type!" );
    }
    val = mxGetScalar( prhs[2] ); /*      mexPrintf( "%f %d\n", val, nVal ); */

    while (pr < (pr0 + nVal)) {
      *pr = val; pr++;
    }
#ifndef FORMATLAB
    plhs[0] = prhs[1];
#endif  
 
  /* asmFull ----------------------------------------------------------
     Dangerous - no bound checks!
     ------------------------------------------------------------------ */
}  else if (!strcmp(CAM,"asm_full"))  {

    int nVal, ii, ir, *pm;
    double *pr, *pv;

    pr = mxGetPr( prhs[1] ); 
    pm = mxGetData( prhs[2] ); 
    pv = mxGetPr( prhs[3] ); 
    
    nVal = mxGetNumberOfElements( prhs[2] );
    for (ii = 0; ii < nVal; ii++) {
      ir = pm[ii];
      if (ir == -1) continue;
      pr[ir] += pv[ii];
    }

  /* insert ----------------------------------------------------------       
     Offsets in C sense!
     ------------------------------------------------------------------ */
}  else if (!strcmp(CAM,"insertval"))  {

    int offR = (int) (*mxGetPr( prhs[2] ));
    int offA = (int) (*mxGetPr( prhs[4] ));
    int nVal = (int) (*mxGetPr( prhs[5] ));
    double *pr, *pa;

/*      mexPrintf( "offR: %d, offA: %d, nVal: %d\n", offR, offA, nVal ); */

    pr = mxGetPr( prhs[1] ) + offR;
    pa = mxGetPr( prhs[3] ) + offA;
    
    if ((mxGetNumberOfElements( prhs[1] ) - offR) < nVal) {
      mexErrMsgTxt( "insert(): result vector too small!" );
    }
    if ((mxGetNumberOfElements( prhs[3] ) - offA) < nVal) {
      mexErrMsgTxt( "insert(): source vector too small!" );
    }
    memcpy( pr, pa, nVal * sizeof( double ) );

  /* mesh_meshGraph ------------------------------------------------------ 
     ------------------------------------------------------------------ */
}  else if (!strcmp(CAM,"meshGraph"))  {

    mxArray *out;
    int32 nNod, *nEl, nGr, nEPMax, *nEP, *conn, *ptrGr;
    int32 nnz, ii;
    int32 *prow, *icol, *aux;
    float64 *val;

    nrhs--; prhs++; /* !!! */
    nNod = (int32) mxGetScalar( prhs[0] );
    nGr = (int32) mxGetScalar( prhs[1] );
    nEPMax = (int32) mxGetScalar( prhs[2] );
    ptrGr = (int32 *) mxGetPr( prhs[3] );
    nEl = (int32 *) mxGetPr( prhs[4] );
    nEP = (int32 *) mxGetPr( prhs[5] );
    conn = (int32 *) mxGetPr( prhs[6] );

    mesh_meshGraph( &nnz, &prow, &icol, nNod, nGr, nEPMax,
		    ptrGr, nEl, nEP, conn );

    out = mxCreateSparse( nNod, nNod, nnz, mxREAL );

    val = mxGetPr( out );
    for (ii = 0; ii < nnz; ii++) {
      val[ii] = 1.0;
    }
    aux = mxGetIr( out );
    if (aux) mxFree( aux );
    mxSetIr( out, icol );

    aux = mxGetJc( out );
    if (aux) mxFree( aux );
    mxSetJc( out, prow );

    plhs[0] = out;

}  else if (!strcmp(CAM,"dofGraph"))  {

    mxArray *out;
    int32 nNod, nEq, *dpn, *dofOffset, *eq;
    int32 dnnz, ii;
    int32 *dprow, *dicol, *nprow, *nicol, *aux;
    float64 *val;

    nrhs--; prhs++; /* !!! */
    nNod = mxGetM( prhs[0] );
    nprow = mxGetJc( prhs[0] );
    nicol = mxGetIr( prhs[0] );

    nEq = (int32) mxGetScalar( prhs[1] );
    dpn = (int32 *) mxGetPr( prhs[2] );
    dofOffset = (int32 *) mxGetPr( prhs[3] );
    eq = (int32 *) mxGetPr( prhs[4] );

    mesh_dofGraph( &dnnz, &dprow, &dicol,
		   nNod, nprow, nicol, nEq, dpn, dofOffset, eq );

    out = mxCreateSparse( nEq, nEq, dnnz, mxREAL );

    val = mxGetPr( out );
    for (ii = 0; ii < dnnz; ii++) {
      val[ii] = 1.0;
    }
    aux = mxGetIr( out );
    if (aux) mxFree( aux );
    mxSetIr( out, dicol );

    aux = mxGetJc( out );
    if (aux) mxFree( aux );
    mxSetJc( out, dprow );

    plhs[0] = out;

/*end of commands */

/* ------------------------------------------------------------------------*/
}  else {mexPrintf("%s\n",CAM);mexErrMsgTxt("Not a valid of_mk call");}
return;

  }
/* -----------------------------------------------------------------------*/
/* calls with commands with <8 characters MODULEF family of elements      */
/* -----------------------------------------------------------------------*/

  if (!mxIsInt32(prhs[2])) {mexErrMsgTxt("integ must be int32");}
  if (!mxIsInt32(prhs[1])) {mexErrMsgTxt("point must be int32");}
  
  point = mxGetData(prhs[1]); 
  integ = mxGetData(prhs[2]);
  constit = mxGetPr(prhs[3]);
  node = mxGetPr(prhs[4]);

  if (nrhs>5)  estate = mxGetPr(prhs[5]);
  else         estate = node;
  if (nrhs>6)  defe = mxGetPr(prhs[6]);
  else         defe = node;
  if (nrhs>7)  eltconst = mxGetPr(prhs[7]);
  else         eltconst = node;
  if (nrhs>8)  InfoAtNode = mxGetPr(prhs[8]);
  else         InfoAtNode = node;


  integ[point[5]+4]=point[4]; /* output type in integ */
  typ=point[4];

  if (point[0]==0 ) { /*     DEAL WITH DEFAULT OUTPUT SIZES */
     if (typ == 0) {
	    point[0] = integ[point[5] +2]*(integ[point[5]+2]+1)/2;
	    point[1] = point[0];
    } else if (typ == 1 || typ == 2 || typ == 3) { /* symmetric one matrix only*/
	    point[0] = integ[point[5] +2]*(integ[point[5]+2]+1)/2;
    } else if (typ == 100) {/*  RHS is given*/
	    point[0] = integ[point[5] +2];
    }
  }

  plhs[0] = mxCreateDoubleMatrix(point[0],1,mxREAL); out = mxGetPr(plhs[0]);
  
  if (nlhs>1)  {
      plhs[1] = mxCreateDoubleMatrix(point[1],1,mxREAL);out1=mxGetPr(plhs[1]);
  } else out1 = out;

  out2 = out;
  
  ierr[0]=' '; ierr[1]=0;
  io_switch(CAM,ierr,point,integ,constit,node,estate,defe,eltconst,out,
      out1,out2,InfoAtNode);

  if (ierr[0]!=' ') mexPrintf("%s\n",ierr);

  return;

}


