
/* ------------------------------------------------------------------------*/
/* Step0 declarations at the beginning of 'MatrixIntegration' */
#if  (MatrixIntegrationStep==0)
 double *fe, *X;
#endif

/* ------------------------------------------------------------------------*/
/* STEP1 prepare for strategy */ 
#if  (MatrixIntegrationStep==1)
/* Follower pressure (fluid/structure coupling) */
else if (!strcmp("fs_matrix",CAM)) { 
       StrategyType=4;
     if (point[4]==5) {
       fe=mxCalloc(3*Nnode,sizeof(double));
       X=mxCalloc(3*Nnode,sizeof(double));
     } else {fe=NULL; X=NULL;}
}
#endif

/* ------------------------------------------------------------------------*/
/* STEP2 switch for strategy implementation  */ 
#if  (MatrixIntegrationStep==2)
case 4: {

if (point[4]==5) {
       /* follower pressure, see implementation in elem0.m*/
       int ji,jj,jk,jl,jq,jqq;
       double Axr[9], Axs[9], *RHS, p, axr[3], axs[3], res[3];

       for (j1=0;j1<9;j1++)  {Axr[j1] = 0.; Axs[j1] = 0.; } 
       for (j1=0;j1<3*Nnode;j1++) fe[j1] = 0.;
       for (j1=0;j1<pow(4*Nnode,2);j1++) ke[j1] = 0.;
       
       for (ji=0;ji<3;ji++) { for(jq=0;jq<Nnode;jq++) {
	   X[ji+3*jq] = def[CurDofPos[ji+4*jq]]+nodeE[jq+Nnode*ji];
       }}

       for (jw=0;jw<Nw;jw++) {/* LOOP ON INTEGRATION POINTS */
	 p=0.;
	 for (jq=0;jq<Nnode;jq++) p+=N[jw+Nw*jq]*def[CurDofPos[3+4*jq]];
	 for (j1=0;j1<3;j1++) {axr[j1] = 0; axs[j1] = 0;}
	 for (ji=0;ji<3;ji++) { for(jq=0;jq<Nnode;jq++) {
	     axr[ji]+= Nr[jw+Nw*jq]*X[ji+3*jq];
	     axs[ji]+= Ns[jw+Nw*jq]*X[ji+3*jq];
         }}
	 Axr[5] = axr[0]; Axr[6] = axr[1]; Axr[1] = axr[2];
	 Axr[7] = -axr[0]; Axr[2] = -axr[1]; Axr[3] = -axr[2];
	 Axs[5] = axs[0]; Axs[6] = axs[1]; Axs[1] = axs[2];
	 Axs[7] = -axs[0]; Axs[2] = -axs[1]; Axs[3] = -axs[2]; 

	 cross(axr,axs,res);

	 for (jq=0;jq<Nnode;jq++) { for (jk=0;jk<3;jk++) {
           fe[3*jq+jk] += p*w[jw]*res[jk]*N[jw+Nw*jq];
	   for (jqq=0;jqq<Nnode;jqq++) { for (jl=0;jl<3;jl++) {
	      ke[(jq+Nnode*jk)+4*Nnode*(jqq+Nnode*jl)] +=
	                     p*w[jw]*N[jw+Nw*jq]*(Axr[3*jl+jk]*Ns[jw+Nw*jqq]
                                                  - Axs[3*jl+jk]*Nr[jw+Nw*jqq]);
           }}
	 }}
       } /* end of the loop on integration points */
       /* assemble RHS*/
       RHS=def+Mdef;
       for (jq=0;jq<Nnode;jq++) { for (jj=0;jj<3;jj++) {
	   RHS[CurDofPos[jj+4*jq]] += fe[jj+3*jq];
	 }}

	 } else {

	   /* this needs to be compiled
 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
 mexErrMsgTxt(" Linear fs_coupling still in elem0.m");*/

 double  coef;
 int     ji, jj, jk, jw;

 for (j1=0;j1<pow(4*Nnode,2);j1++) ke[j1] = 0.;

 for (jw=0;jw<Nw;jw++) {
   for (ji=0;ji<3;ji++) {
    coef=bas[6+ji+jw*9]*jdet[jw]*w[jw];
    if (coef!=0.) {
      for (jj=0; jj<Nnode;jj++) { for (jk=0; jk<Nnode;jk++){ 
        if (point[4]==1) {
         ke[Nnode*ji+jj+4*Nnode*(jk+Nnode*3)] -= coef*N[jw+Nw*jj]*N[jw+Nw*jk];
        } else if (point[4]==2) {
         ke[jk+Nnode*3+4*Nnode*(Nnode*ji+jj)] += coef*N[jw+Nw*jj]*N[jw+Nw*jk];
        } /* if */
      }} /* jk, jj */
    } /* if coef */

 }} /* ji, jw */



	 }
}
break; 

#endif

/* ------------------------------------------------------------------------*/
/* STEP3 possible cleanup if needed */

#if  (MatrixIntegrationStep==3)
 case 4: {
   if (fe!=NULL) mxFree(fe); 
   if (X!=NULL) mxFree(X);
   break;
 }
#endif
