#_______________________________________________________________________________
#
# COMPARACIÓN DE LA DESCOMPOSICIÓN POLAR Y LA DESCOMPOSICIÓN POLAR INFINITESIMAL
#
#   Este programa realiza la descomposición polar y su versión infinitesimal
#   para un tensor gradiente de movimientos (F=I+J) en 2D/3d arbitrario.
#
#   Se introduce como dato el tensor gradiente de desplazamientos (J).   
#
#   NOTA IMPORTANTE:
# 
#      Los resultados de la descomposición infinitesimal sólo son fiables cuando
#      la norma del tensor J es mucho menor que la norma de la matriz identidad.
#
#_______________________________________________________________________________



# Lectura del tensor J

J = input("Introduce el tensor J = [,,,] = ");

[nf, nc] = size(J);

if (nf != nc)
   printf("Error: Solo se admiten matrices cuadradas")
   return
end

n = nf;
if ((n < 2) || (n > 4))
  printf("Error: Solo se admiten matrices 2x2 o 3x3\n")
  return
end

printf("\n")

I = eye(n);

F = I+J;


# Descomposicion Polar

printf("DESCOMPOSICIÓN POLAR:\n")

f = det(F)                             % coeficiente de expansión volúmetrica
if (f <= 0)
  printf("Error: el determinante de F = I + J tiene que ser positivo\n")
  return
end

C = F' * F;                            % C = Tensor derecho de Cauchy-Green
[U, Lambda2] = eig(C);                 % Autovectores y autovalores de C
Lambda = sqrt(Lambda2);                % Raíz cuadrada de los autovalores de C
E = U * (Lambda - I) * U'              % Tensor de deformación de Biot
R = F * U * inv(Lambda) * U'           % Tensor de rotación
cdl = f^(1/n)                          % Coeficiente de dilatación lineal
Ed = 1/cdl * (E + (1 - cdl) * I)       % Tensor de distorsión
In =   cdl * I                         % Tensor de inflación


# Descomposicion Polar Infinitesimal

printf("DESCOMPOSICIÓN POLAR INFINITESIMAL:\n")

trj   = trace(J)
f_inf = 1 + trj                         % coeficiente de expansión volumétrica infinitesimal

E_inf =     (J+J')/2                    % Tensor de deformación infinitesimal
R_inf = I + (J-J')/2                    % Tensor de rotación infinitesimal
cdl_inf = 1+trj/n                       % Coeficiente de dilatación lineal infinitesimal
Ed_inf = E_inf - trj/n * I              % Tensor de distorsión infinitesimal
In_inf = (1+trj/n) * I                  % Tensor de inflación infinitesimal



