/*     CalculiX - A 3-dimensional finite element program                   */
/*              Copyright (C) 1998-2011 Guido Dhondt                          */
/*     This program is free software; you can redistribute it and/or     */
/*     modify it under the terms of the GNU General Public License as    */
/*     published by the Free Software Foundation(version 2);    */
/*                    */

/*     This program is distributed in the hope that it will be useful,   */
/*     but WITHOUT ANY WARRANTY; without even the implied warranty of    */ 
/*     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the      */
/*     GNU General Public License for more details.                      */

/*     You should have received a copy of the GNU General Public License */
/*     along with this program; if not, write to the Free Software       */
/*     Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.         */

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <string.h>
#include "CalculiX.h"

#ifdef SPOOLES
   #include "spooles.h"
#endif
#ifdef SGI
   #include "sgi.h"
#endif
#ifdef TAUCS
   #include "tau.h"
#endif
#ifdef PARDISO
   #include "pardiso.h"
#endif

void storecontactdof(int *nope,int *nactdof, int *mt, int *konl, int **ikactcontp, 
          int *nactcont,int *nactcont_, double *bcont, double *fnl, 
          int *ikmpc, int *nmpc, int *ilmpc,int *ipompc, int *nodempc, 
          double *coefmpc){

  int j,j1,jdof,id,k,l,ist,index,node,ndir,*ikactcont=*ikactcontp;

  for(j=0;j<*nope;j++){
      for(j1=0;j1<3;j1++){
	  jdof=nactdof[*mt*(konl[j]-1)+j1+1];
	  if(jdof!=0){
	      
	      jdof--;
	      FORTRAN(nident,(ikactcont,&jdof,nactcont,&id));
	      do{
		  if(id>0){
		      if(ikactcont[id-1]==jdof){
			  break;
		      }
		  }
		  (*nactcont)++;
		  if(*nactcont>*nactcont_){
		      *nactcont_=(int)(1.1**nactcont_);
		      RENEW(ikactcont,int,*nactcont_);
		  }
		  k=*nactcont-1;
		  l=k-1;
		  while(k>id){
		      ikactcont[k--]=ikactcont[l--];
		  }
		  ikactcont[id]=jdof;
		  break;
	      }while(1);
	      
	      bcont[jdof]-=fnl[3*j+j1];
	  }else{
	      jdof=8*(konl[j]-1)+j1+1;
	      FORTRAN(nident,(ikmpc,&jdof,nmpc,&id));
	      if(id>0){
		  if(ikmpc[id-1]==jdof){
		      id=ilmpc[id-1];
		      ist=ipompc[id-1];
		      index=nodempc[3*ist-1];
		      if(index==0) continue;
		      do{
			  node=nodempc[3*index-3];
			  ndir=nodempc[3*index-2];
			  jdof=nactdof[*mt*(node-1)+ndir];
			  if(jdof!=0){
			      
			      jdof--;
			      FORTRAN(nident,(ikactcont,&jdof,nactcont,&id));
			      do{
				  if(id>0){
				      if(ikactcont[id-1]==jdof){
					  break;
				      }
				  }
				  (*nactcont)++;
				  if(*nactcont>*nactcont_){
				      *nactcont_=(int)(1.1**nactcont_);
				      RENEW(ikactcont,int,*nactcont_);
				  }
				  k=*nactcont-1;
				  l=k-1;
				  while(k>id){
				      ikactcont[k--]=ikactcont[l--];
				  }
				  ikactcont[id]=jdof;
				  break;
			      }while(1);
			      
/*			      bcont[jdof]+=coefmpc[index-1]*
			      fnl[3*j+j1]/coefmpc[ist-1];*/
			      bcont[jdof]-=coefmpc[index-1]*
				  fnl[3*j+j1]/coefmpc[ist-1];
			  }
			  index=nodempc[3*index-1];
			  if(index==0) break;
		      }while(1);
		  }
	      }
	  }
      }
  }

  *ikactcontp=ikactcont;
  
  return;
}

