#include <stdio.h>
#include <math.h>
//Resuelve la ec de la elasticidad:
// Lapl(u) + A Grad(Div u)+ B =0
//u, B = vectores real
// A = 1/(1-2*nu) =2 para nu = 0.25
//Ponemos H = 1 m  
// Bx = 0  
//Hacemos Dx=Dy=H/Ny=> L = 2*(Nx-1) H/Ny
//By =-Dx^2*ro*g*2*(1+nu)/Y: gravedad
// Ny = nº de intervalos en y
// Nº intervalos x = 2*(Nx-1) (simetria) 
main()
{
  double ux[101][101],uy[101][101];
  double nu=0.25,Y=1.e5,ro=1000.,A,B,C,g=9.8,H=1.,L;
  int i,j,Nx=7,Ny=12;// L = 2*(Nx-1)*Dx, H=Ny*Dx
  int t=0; // "tiempo"
  char carac;
  // Inicializa valores
  A=1./(1.-2.*nu);
  B = -ro*g*2.*(1.+nu)/Y*H*H/Ny/Ny; 
  L=2*(Nx-1)*H/Ny;
  for(i=0;i<=Nx+1;i++) for(j=0;j<=Ny+1;j++)
  {ux[i][j]=0;uy[i][j]=0;
  }
  // Ciclos para solucion del sistema, metodo Gauss-Seidel
   while(carac!='0')
  {
   //Comenzamos en x = L/2 e y =0 para estabilidad del método de Gauss-Seidel
   for(i=Nx;i>=1;i--) for(j=1;j<=Ny;j++)
   {  // Se fija ux(L/2,y)=0 por simetría.
    if(i<Nx) ux[i][j]=(1+A)/(4+2*A)*(ux[i+1][j]+ux[i-1][j])
             +1./(4+2*A)*(ux[i][j+1]+ux[i][j-1])
             + A/4/(4+2*A)*(uy[i+1][j+1]+uy[i-1][j-1]-uy[i+1][j-1]-uy[i-1][j+1]);
    uy[i][j]=(1+A)/(4+2*A)*(uy[i][j+1]+uy[i][j-1])
             +1./(4+2*A)*(uy[i+1][j]+uy[i-1][j])
             +A/4/(4+2*A)*(ux[i+1][j+1]+ux[i-1][j-1]-ux[i+1][j-1]-ux[i-1][j+1])
             +B/(4+2*A);            
   }
    //Contorno
    // Suelo, y =0
    for(i=0;i<=Nx+1;i++) uy[i][0]=0;
    for(i=Nx-1;i>=0;i--) ux[i][0]=ux[i+1][0]+uy[i][1]-uy[i][0];
    ux[Nx+1][0]=ux[Nx-1][0];
    // x = 0, lado libre
  
   for(j=1;j<=Ny+1;j++) ux[0][j]=ux[1][j]; //sigma xx =0 
   for(j=1;j<=Ny;j++) uy[0][j]=uy[0][j-1]-ux[1][j]+ux[0][j]; //sigma xy = 0
   
   //Parte de arriba: Ny+1
 
   for(i=Nx;i>=0;i--) uy[i][Ny+1]=uy[i][Ny]+ux[0][j]-ux[0][j-1];//sigma yy =0 
    for(i=Nx;i>=0;i--) ux[i][Ny+1]=ux[i][Ny]-uy[i+1][Ny+1]+uy[i][Ny+1];//sigma xy =0
      //centro: simetria, ux(L/2,y)=0 fijado al inicio y además:
   for(j=0;j<=Ny;j++)
   {uy[Nx+1][j]=uy[Nx-1][j]; 
    ux[Nx+1][j]=-ux[Nx-1][j]; //innecesario pero clarifica
   } 

      // escribe resultados
    if(t%(Nx*Nx)==0)
   { 
   printf("\n\n No de pasos: %i   datos de ux \n",t);
   for(j=Ny;j>=0;j--)
   { printf("\n");
     for(i=0;i<=Nx+1;i++)printf("%5.0lf\t",1.e5*ux[i][j]);
     for(i=Nx-2;i>=0;i--)printf("%5.0lf\t",-1.e5*ux[i][j]);
   }
     printf("\n\n teclea 0 y <ret> para acabar: ");
     carac=getchar();
   }
  
  t++;
  }
    printf("\n\n Resultados finales: ux* 1.e5");
    for(j=Ny;j>=0;j--)
   { printf("\n");
     for(i=1;i<=Nx;i++)printf("%5.0lf\t",1.e5*ux[i][j]);
     //for(i=Nx-1;i>=1;i--)printf("%5.0lf\t",-1.e5*ux[i][j]);
   }
    printf("\n\n Resultados finales: uy* 1.e5");
    for(j=Ny;j>=0;j--)
   { printf("\n");
     for(i=1;i<=Nx;i++)printf("%5.0lf\t",1.e5*uy[i][j]);
     //for(i=Nx-1;i>=1;i--)printf("%5.0lf\t",1.e5*uy[i][j]);
   }
   getchar();getchar();
   
   printf("\n\n Solucion analitica (Landau, Elasticidad  pág 28");
   printf("\n\n ux:"); 
   // nu*ro*g/Y*Dx*Dy = B*(A+1)/2/(A+3) 
   // nu = (A+1)/2
   C=nu*ro*g/Y;
      for(j=Ny;j>=0;j--)
   { printf("\n");
     for(i=1;i<=Nx;i++)printf("%5.0lf\t",1.e5*C*H*(Ny-j)/Ny*L*(i-Nx)/2./Nx);
   }
      printf("\n\n uy:"); 
    for(j=Ny;j>=0;j--)
   { printf("\n");
     for(i=1;i<=Nx;i++)
     printf("%5.0lf\t",-1.e5*C/2/nu*(H*H*(1.-(Ny-j)*(Ny-j)/Ny/Ny)-nu*L*(i-Nx)/2./Nx*L*(i-Nx)/2./Nx));
   }
getchar();getchar();
 
}