////////////////////////////////////////////////////////////
// set up the initial u for main program
// inputs: u_interior :  values of deformation on each node 
// Written by Weigang Zhong  11/15/2003
////////////////////////////////////////////////////////////

#include"martensite.h"
//#include"time.h"

#define eps 1.0e-2


// the first choice of initial guess of deformation
void set_init_1 (double *u_interior)
{
  int i, j;
  double *x_local, **F;
  double r;
  
  x_local = vector (1, 2);
  F = dmatrix(1, 2, 1, 2);

  F[1][1] = 1.0;
  F[1][2] = (1.0 - 2.0 * lambda) * epsilon;
  F[2][1] = 0.0;
  F[2][2] = 1.0;
	 
  for (i = 1; i < N; i++)
    {
      for (j = 1; j < N; j++)
	{
	  x_local[1] = X0 + H * i;
	  x_local[2] = Y0 + H * j;
	  r = randf();
//	  r = .8;
          //srand( (unsigned)time( NULL ) ); 
	  //r = (double)(rand()/(double)RAND_MAX);
	 
	  *(u_interior + 2 * (i - 1) * (N - 1) + 2 * (j - 1) + 1) = 
		  F[1][1] * x_local[1] + F[1][2] * x_local[2] + eps * (r - .5);
	  *(u_interior + 2 * (i - 1) * (N - 1) + 2 * (j - 1) + 2) = 
	   	  F[2][1] * x_local[1] + F[2][2] * x_local[2] + eps * (r - .5);
	}
    }

  return;
}

// alternate choice of inital guess of deformation
void set_init_2 (double *u_interior)
{
  int i, j;

  for (i = 1; i < N; i++)
    {
      for (j = 1; j < N; j++)
	{
	  u_interior[2 * (i - 1) * (N - 1) + 2 * (j - 1) + 1] = (H * j + .5*i);
	  u_interior[2 * (i - 1) * (N - 1) + 2 * (j - 1) + 2] = (.5 * H * i + 1);
	}
    }

  return;
}

//  laminate guess
void set_init_3 (double *u_interior)
{
  int i, j, k, m, n;
  double *x_local, *x_boundary, *x_internal, *u;
  double **F, **Ua, **Ub;
  double r;
  
  x_local    = vector (1, 2);
  x_boundary = vector (1, 2);
  x_internal = vector (1, 2);
  F          = dmatrix(1, 2, 1, 2);
  Ua         = dmatrix(1, 2, 1, 2);
  Ub         = dmatrix(1, 2, 1, 2);
  u          = vector (1, 2 * (N + 1) * (N + 1));

  F[1][1] = 1.0;
  F[1][2] = (1.0 - 2.0 * lambda) * epsilon;
  F[2][1] = 0.0;
  F[2][2] = 1.0;

  Ua[1][1] = 1.0;
  Ua[1][2] = epsilon;
  Ua[2][1] = 0.0;
  Ua[2][2] = 1.0;
  Ub[1][1] = 1.0;
  Ub[1][2] = -1 * epsilon;
  Ub[2][1] = 0.0;
  Ub[2][2] = 1.0;

  k = delta1 + delta2;
  n = (N - 2 * delta - 1) / k;
  cout << "\nnumber is " << n << "\n";
  
  for (i = delta + 1; i < N - delta; i++)
    {
      for (m = 0; m < n; m++)
        {
          for (j = delta + 1 + m * k; j <= delta + delta1 + m * k; j++)
            {
	      x_local[1] = X0 + H * i;
	      x_local[2] = Y0 + H * j;
	      
	      *(u_interior + 2 * (i - 1) * (N - 1) + 2 * (j - 1) + 1) =
  		  Ua[1][1] * x_local[1] + Ua[1][2] * x_local[2] ;
	      *(u_interior + 2 * (i - 1) * (N - 1) + 2 * (j - 1) + 2) =
  		  Ua[2][1] * x_local[1] + Ua[2][2] * x_local[2] ;
            }
          for (j = delta + delta1 + m * k + 1; j <= delta + (m + 1) * k; j++)
            {
	      x_local[1] = X0 + H * i;
	      x_local[2] = Y0 + H * j;
	      
	      *(u_interior + 2 * (i - 1) * (N - 1) + 2 * (j - 1) + 1) =
  		  Ub[1][1] * x_local[1] + Ub[1][2] * x_local[2] ;
	      *(u_interior + 2 * (i - 1) * (N - 1) + 2 * (j - 1) + 2) =
  		  Ub[2][1] * x_local[1] + Ub[2][2] * x_local[2] ;
            } 
	}
      
      for (j = delta + k * n + 1; j < N - delta; j++) 
        {
	  x_local[1] = X0 + H * i;
	  x_local[2] = Y0 + H * j;

          *(u_interior + 2 * (i - 1) * (N - 1) + 2 * (j - 1) + 1) =
  	      Ua[1][1] * x_local[1] + Ua[1][2] * x_local[2] ;
	  *(u_interior + 2 * (i - 1) * (N - 1) + 2 * (j - 1) + 2) =
  	      Ua[2][1] * x_local[1] + Ua[2][2] * x_local[2] ;
	}	
  } 

  u = vector (1, 2 * (N + 1) * (N + 1)); 
  set_boundary (u_interior, u);
  
  for (i = delta + 1; i < N - delta; i++)
    {	  
      x_boundary[1] = X0 + H * i;
      x_boundary[2] = Y0;
      x_internal[1] = X0 + H * i;
      x_internal[2] = Y0 + H * (delta + 1);
      
      for (j = 1; j<=delta; j++)
        { 
	  x_local[1] = X0 + H * i;
          x_local[2] = Y0 + H * j;  
	  
          *(u_interior + 2 * (i - 1) * (N - 1) + 2 * (j - 1) + 1) =
		  (u[2 * i * (N + 1) + 1] - u[2 * i * (N + 1) + 2 * delta + 3]) /
		    (sqrt((DSQR(x_local[2] - x_internal[2]) + DSQR(x_local[1] - x_internal[1])) / 
		      (DSQR(x_local[2] - x_boundary[2]) + DSQR(x_local[1] - x_boundary[1]))) - 1) +
		  u[2 * i * (N + 1) + 1];
          
	  *(u_interior + 2 * (i - 1) * (N - 1) + 2 * (j - 1) + 2) =
		  (u[2 * i * (N + 1) + 2] - u[2 * i * (N + 1) + 2 * delta + 4]) * 
		    ((*(u_interior + 2 * (i - 1) * (N - 1) + 2 * (j - 1) + 1)) - u[2 * i * (N + 1) + 2 * delta + 3]) / 
		      (u[2 * i * (N + 1) + 1] - u[2 * i * (N + 1) + 2 * delta + 3]) + 
		  u[2 * i * (N + 1) + 2 * delta + 4];

	  *(u + 2 * i * (N + 1) + 2 * j + 1) = 
	          *(u_interior + 2 * (i - 1) * (N - 1) + 2 * (j - 1) + 1);
	  
          *(u + 2 * i * (N + 1) + 2 * j + 2) = 
	          *(u_interior + 2 * (i - 1) * (N - 1) + 2 * (j - 1) + 2);
        }
      
      x_boundary[1] = X0 + H * i;
      x_boundary[2] = Y0 + H * N;
      x_internal[1] = X0 + H * i;
      x_internal[2] = Y0 + H * (N - delta - 1);
  	
      for (j = N - delta; j <= N - 1; j++)
        {
          x_local[1] = X0 + H * i;
          x_local[2] = Y0 + H * j; 
	  
	  *(u_interior + 2 * (i - 1) * (N - 1) + 2 * (j - 1) + 1) =
	           (u[2 * i * (N + 1) + 2 * N + 1] - u[2 * i * (N + 1) + 2 * N - 2 * delta - 1]) /
		     (sqrt((DSQR(x_local[2] - x_internal[2]) + DSQR(x_local[1] - x_internal[1])) /
			(DSQR(x_local[2] - x_boundary[2]) + DSQR(x_local[1] - x_boundary[1]))) - 1) +
		   u[2 * i * (N + 1) + 2 * N + 1];
		            
          *(u_interior + 2 * (i - 1) * (N - 1) + 2 * (j - 1) + 2) =
	           (u[2 * i * (N + 1) + 2 * N + 2] - u[2 * i * (N + 1) + 2 * N - 2 * delta]) * 
		     ((*(u_interior + 2 * (i - 1) * (N - 1) + 2 * (j - 1) + 1)) - u[2 * i * (N + 1) + 2 * N - 2 * delta - 1]) / 
		       (u[2 * i * (N + 1) + 2 * N + 1] - u[2 * i * (N + 1) + 2 * N - 2 * delta - 1]) + 
		   u[2 * i * (N + 1) + 2 * N - 2 * delta];

	  *(u + 2 * i * (N + 1) + 2 * j + 1) =                     
                  *(u_interior + 2 * (i - 1) * (N - 1) + 2 * (j - 1) + 1);

	  *(u + 2 * i * (N + 1) + 2 * j + 2) =
	          *(u_interior + 2 * (i - 1) * (N - 1) + 2 * (j - 1) + 2);			  
	}
    }

  
  for (j = delta + 1; j < N - delta; j++)
    {	  
      x_boundary[1] = X0;
      x_boundary[2] = Y0 + H * j;
      x_internal[1] = X0 + H * (delta + 1);
      x_internal[2] = Y0 + H * j;
      
      for (i = 1; i<=delta; i++)
        { 
	  x_local[1] = X0 + H * i;
          x_local[2] = Y0 + H * j;  
	  
          *(u_interior + 2 * (i - 1) * (N - 1) + 2 * (j - 1) + 1) =
		  (u[2 * j + 1] - u[2 * (delta + 1) * (N + 1) + 2 * j + 1]) /
		    (sqrt((DSQR(x_local[2] - x_internal[2]) + DSQR(x_local[1] - x_internal[1])) / 
		      (DSQR(x_local[2] - x_boundary[2]) + DSQR(x_local[1] - x_boundary[1]))) - 1) +
		  u[2 * j + 1];
          
	  *(u_interior + 2 * (i - 1) * (N - 1) + 2 * (j - 1) + 2) =
		  (u[2 * j + 2] - u[2 * (delta + 1) * (N + 1) + 2 * j + 2]) * 
		    ((*(u_interior + 2 * (i - 1) * (N - 1) + 2 * (j - 1) + 1)) - u[2 * (delta + 1) * (N + 1) + 2 * j + 1]) / 
		      (u[2 * j + 1] - u[2 * (delta + 1) * (N + 1) + 2 * j + 1]) + 
		  u[2 * (delta + 1) * (N + 1) + 2 * j + 2];

	  *(u + 2 * i * (N + 1) + 2 * j + 1) = 
	          *(u_interior + 2 * (i - 1) * (N - 1) + 2 * (j - 1) + 1);
	  
          *(u + 2 * i * (N + 1) + 2 * j + 2) = 
	          *(u_interior + 2 * (i - 1) * (N - 1) + 2 * (j - 1) + 2);
        }
      
      x_boundary[1] = X0 + H * N;
      x_boundary[2] = Y0 + H * j;
      x_internal[1] = X0 + H * (N - delta - 1);
      x_internal[2] = Y0 + H * j;
  	
      for (i = N - delta; i <= N - 1; i++)
        {
          x_local[1] = X0 + H * i;
          x_local[2] = Y0 + H * j; 
	  
	  *(u_interior + 2 * (i - 1) * (N - 1) + 2 * (j - 1) + 1) =
		  (u[2 * N * (N + 1) + 2 * j + 1] - u[2 * (N - delta - 1) * (N + 1) + 2 * j + 1]) /
		    (sqrt((DSQR(x_local[2] - x_internal[2]) + DSQR(x_local[1] - x_internal[1])) / 
		      (DSQR(x_local[2] - x_boundary[2]) + DSQR(x_local[1] - x_boundary[1]))) - 1) +
		  u[2 * N * (N + 1) + 2 * j + 1];
          
	  *(u_interior + 2 * (i - 1) * (N - 1) + 2 * (j - 1) + 2) =
		  (u[2 * N * (N + 1) + 2 * j + 2] - u[2 * (N - delta - 1) * (N + 1) + 2 * j + 2]) * 
		    ((*(u_interior + 2 * (i - 1) * (N - 1) + 2 * (j - 1) + 1)) - u[2 * (N - delta - 1) * (N + 1) + 2 * j + 1]) / 
		      (u[2 * N * (N + 1) + 2 * j + 1] - u[2 * (N - delta - 1) * (N + 1) + 2 * j + 1]) + 
		  u[2 * (N - delta - 1) * (N + 1) + 2 * j + 2];

	  *(u + 2 * i * (N + 1) + 2 * j + 1) = 
	          *(u_interior + 2 * (i - 1) * (N - 1) + 2 * (j - 1) + 1);
	  
          *(u + 2 * i * (N + 1) + 2 * j + 2) = 
	          *(u_interior + 2 * (i - 1) * (N - 1) + 2 * (j - 1) + 2);
	}
    }



  for (i = N - delta; i <= N - 1; i++)
    for (j = N - delta; j <= N - 1; j++)
      {
        x_local[1] = X0 + H * i;
	x_local[2] = Y0 + H * j;
	r = randf();
	cout << "\nr is " << r << "\n";
        *(u_interior + 2 * (i - 1) * (N - 1) + 2 * (j - 1) + 1) =
                F[1][1] * x_local[1] + F[1][2] * x_local[2] + eps * (r - .5);
        *(u_interior + 2 * (i - 1) * (N - 1) + 2 * (j - 1) + 2) =
                F[2][1] * x_local[1] + F[2][2] * x_local[2] + eps * (r - .5);
      }

  
  return;
}