/***************************************************************
File:   3ds.c - 3-Dimensional Stress Tensor
Desc:   Given a stress tensor, with normal stresses Sxx, Syy,
        and Szz, shear stresses, Sxy, Sxz and Syz, the stress
        transformation is calculated to give principal stresses 
        and the corresponding transformation matrix.
        Other information is provided such as von Mises stress,
        octahedral stress, max shear and the deviatoric stress.

        David Edwards
        Rev 1.1 - Sat Dec 30 15:07:53 EST 2000
***************************************************************/

#include <stdio.h>
#include <math.h>
#define NRANSI
#include "nr.h"
#include "nrutil.h"

void input (char *str, double *addr)
{
    printf(str);
    scanf("%lf", addr);
}

double sq (double x)
{
    return (x*x);
}

double r2d (double rad)
{
    return (rad*180.0/M_PI);
}

int main(void)
{
    double sxx,syy,szz,sxy,sxz,syz;
    double **s,*sp,**vcol,**v;
    double d,Sv,sum3,Tauoct,sum;
    int nrot,i,j;

    s=dmatrix(1,3,1,3);     // Input stress tensor
    sp=dvector(1,3);        // Principal stresses
    vcol=dmatrix(1,3,1,3);  // Transformation matrix (column filled)
    v=dmatrix(1,3,1,3);     // Transformation matrix (row filled)

    input("\n  Tensile Stress, Sxx [ksi]: ",&sxx);
    input("  Tensile Stress, Syy [ksi]: ",&syy);
    input("  Tensile Stress, Szz [ksi]: ",&szz);
    input("\n  Shear Stress,   Sxy [ksi]: ",&sxy);
    input("  Shear Stress,   Sxz [ksi]: ",&sxz);
    input("  Shear Stress,   Syz [ksi]: ",&syz);

    s[1][1]=sxx; s[2][2]=syy; s[3][3]=szz;
    s[1][2]=s[2][1]=sxy; s[1][3]=s[3][1]=sxz; s[2][3]=s[3][2]=syz;
    
    // Print the input stress tensor
    printf("\n  Input stress tensor [ksi]:");
    for (i=1;i<=3;i++) {
        printf("\n    %9.2lf %9.2lf %9.2lf",s[i][1],s[i][2],s[i][3]);
    }
    
    // Find the principal stresses and transformation matrix
    jacobi(s,3,sp,vcol,&nrot);
    eigsrt(sp,vcol,3);

    // v = transpose(vcol)
    for (i=1;i<=3;i++) {
        for (j=1;j<=3;j++) 
            v[i][j]=vcol[j][i];
    }
    
    // Make sure v is proper orhtogonal (i.e. obeys right hand rule)
    d=det3(v);
    if ((1-fabs(d)) > 1.0e-8)
            nrerror("V is not proper orthogonal");

    if (d < 0.0) for (j=1;j<=3;j++) v[3][j] = -v[3][j];


    // Reset s for deviatoric tensor (jacobi messed it up)
    s[1][1]=sxx; s[2][2]=syy; s[3][3]=szz;
    s[1][2]=s[2][1]=sxy; s[1][3]=s[3][1]=sxz; s[2][3]=s[3][2]=syz;

    //
    // Print all the results
    //
    printf("\n\n  Principal Stresses [ksi]:");
    printf("\n    S1 = %.2lf  S2 = %.2lf  S3 = %.2lf",sp[1],sp[2],sp[3]);

    printf("\n\n  Maximum Shear Stress [ksi]:     %7.2lf",.5*(sp[1]-sp[3]));
        
    Sv=sqrt(.5*(sq(sp[1]-sp[2])+sq(sp[1]-sp[3])+sq(sp[2]-sp[3])));
    printf("\n\n  von Mises Stress [ksi]:         %7.2lf",Sv);


    sum3=(sp[1]+sp[2]+sp[3])/3.0;
    Tauoct=sqrt(sq(sp[1]-sp[2])+sq(sp[1]-sp[3])+sq(sp[2]-sp[3]))/3.0;
    printf("\n\n  Octahedral Normal Stress [ksi]: %7.2lf",sum3);
    printf("\n  Octahedral Shear  Stress [ksi]: %7.2lf",Tauoct);

    // Create deviatoric stress tensor
    for (i=1;i<=3;i++) {
        s[i][i] -= sum3;
    }
    
    sum=0.0;
    for (i=1;i<=3;i++) {
        for (j=1;j<=3;j++) 
            sum += s[i][j]*s[i][j];
    }
    printf("\n\n  Deviatoric Stress [ksi]:        %7.2lf",sqrt(sum));

    // Output the transformation matrix in two forms
    printf("\n\n Principal Transformation Matrix (direction cosines):");
    for (i=1;i<=3;i++) {
        printf("\n    %9.4lf %9.4lf %9.4lf",v[i][1],v[i][2],v[i][3]);
    }

    printf("\n\n  Transformation angles [deg]:\n"
           "  (NOTE: Vij = angle between ei' and ej)");
    for (i=1;i<=3;i++) {
        printf("\n    ");
        for (j=1;j<=3;j++) 
            printf("%9.2lf ",r2d(acos(v[i][j])));
    }
    
    printf("\n\n");

    free_dmatrix(s,1,3,1,3);
    free_dvector(sp,1,3);
    free_dmatrix(vcol,1,3,1,3);
    free_dmatrix(v,1,3,1,3);
}