/****************************************************************************
*
*                                 M U E S L I   v 1.8
*
*
*     Copyright 2020 IMDEA Materials Institute, Getafe, Madrid, Spain
*     Contact: muesli.materials@imdea.org
*     Author: Ignacio Romero (ignacio.romero@imdea.org)
*
*     This file is part of MUESLI.
*
*     MUESLI 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, either version 3 of the License, or
*     (at your option) any later version.
*
*     MUESLI 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 MUESLI.  If not, see <http://www.gnu.org/licenses/>.
*
****************************************************************************/

#include "interface_abaqus.h"
#include "muesli/muesli.h"

using namespace muesli;



// matlabel(in)   material label for muesli library. This is a label that
//                is set in the file below, and decided in muesli
//                when using abaqus, the label must be indicated in the .inp file
// STRESS: (in/out) in is Cauchy stress (11,22,33,12,13,23) at tn
//                  output is updated Cauchy stress at tn+1.
// STATEV: (in/out) in is state variables in the last converged solution tn
//                  output is updated state variables
// CMNAME  (in)  Material name
// NOEL    (in)  Element number
// NPT     (in)  Gauss point number
// COORDS  (in)  Gauss point coordinates
// PROPS   (in)  Material properties
// NPROPS  (in)  Dimensions of PROPS
// CELENT  (in)  Characteristic elmt length
// TIME    (in)  TIME[0]: time in current step, TIME[1]: total time
// DTIME   (in)  Time step increment (in)
// PNEWDT  (out) 1: if step converged, <1: dt reduction requested (output)
// STRAN:  (in)  Strain at tn (11,22,33,12,13,23) for small strain problems
//                  Corrotational strain (11,22,33,12,13,23) for large strain problems
// DSTRAN: (in)  Mechanical strain increment (11,22,33,12,13,23) for small strain problems
//                  Mechanical corrotational strain increment (11,22,33,12,13,23) for large strain problems
// DFGRD0: (in)  Deformation gradient at time tn
// DFGRD1: (in)  Deformation gradient at current time tn+1
// TEMP:   (in)  Temperature at time tn
// DTEMP:  (in)  Temperature increment from tn to tn+1
// NTENS   (in)  size of stress array (e.g., 6, 3)
// DDSDDE  (in/out) dimensions NTENS x NTENS
//               Cabaqus_ijkl = Cspatial_ijkl + delta_ij sigma_jl + delta_jl sigma_ik
// NSTATEV:(in)  Number of state variables
// SSE:    (in/out) in input, elastic energy in the last converged solution tn
//                  in output, current elastic energy
// SPD:    (in/out) in input, dissipated energy in the last converged solution tn
//                  in output, current dissipated energy
// KSTEP:  (in)     step number
// KINC:   (in)     increment number, within the current step
// RPL:    (out)    heat generated by unit time
// DDSDDT: (out)    stress rate. Dimensions NTENS
// DRPLDE  (out)    Derivative of heat generation w.r.t. strain. Dimensions NTENS
// DRPLDT: (out)    Derivative of heat generation w.r.t. temperature
// PREDEF: (in)     Prededefined field at the beginning of the step (array)
// DPRED:  (in)     Prededefined field increments (array)
// NDI:    (in)     Number of direct stress components (normal stress / longitudinal strains)
// NSHR    (out)    Number of shear stress components (angular strains)
// NTENS   (in)     NDI+NSHR
// DROT    (in)     Rotation incrementation
// LAYER   (in)     Number of layers for composited shells or layered solids
// KSPT    (in)     Number current section point number within layer

//


static void abaqusFiniteStrainTangent(const istensor& Cauchy, const itensor4& cspatial, itensor4& cabaqus);



int interface_abaqus_(int*    matlabel,
                      double* STRESS,
                      double* STATEV,
                      double* DDSDDE,
                      double* SSE,
                      double* SPD,
                      double* SCD,
                      double* RPL,
                      double* DDSDDT,
                      double* DRPLDE,
                      double* DRPLDT,
                      double* STRAN,
                      double* DSTRAN,
                      double* TIME,
                      double* DTIME,
                      double* TEMP,
                      double* DTEMP,
                      double* PREDEF,
                      double* DPRED,
                      char    CMNAME[80],
                      int*    NDI,
                      int*    NSHR,
                      int*    NTENS,
                      int*    NSTATEV,
                      double* PROPS,
                      int*    NPROPS,
                      double* COORDS,
                      double* DROT,
                      double* PNEWDT,
                      double* CELENT,
                      double* DFGRD0,
                      double* DFGRD1,
                      int*    NOEL,
                      int*    NPT,
                      int*    LAYER,
                      int*    KSPT,
                      int*    KSTEP,
                      int*    KINC)
{
    double tangent[6][6];
    
    switch (*matlabel)
    {
        case 1:
        {
            // linear elastic material
            //
            // matproperties[0] : E
            // matproperties[1] : nu
            const double E     = PROPS[0];
            const double nu    = PROPS[1];
            const double rho   = 1.0;

            elasticIsotropicMaterial *theMaterial = new elasticIsotropicMaterial("an elastic material", E, nu, rho);
            elasticIsotropicMP       *thePoint    = dynamic_cast<elasticIsotropicMP*>(theMaterial->createMaterialPoint());

            // set converged and current states
            istensor epsc(    STRAN[0]+DSTRAN[0],      STRAN[1]+DSTRAN[1],      STRAN[2]+DSTRAN[2],
                          .5*(STRAN[5]+DSTRAN[5]), .5*(STRAN[4]+DSTRAN[4]), .5*(STRAN[3]+DSTRAN[3]));

            double timec(TIME[0]);
            thePoint->updateCurrentState(timec, epsc);
            thePoint->commitCurrentState();

            // compute all the required fields
            double stressMuesli[6];
            istensor stressT;
            thePoint->stress(stressT);
            ContraContraSymTensorToVector(stressT, stressMuesli);

            STRESS[0] = stressMuesli[0];
            STRESS[1] = stressMuesli[1];
            STRESS[2] = stressMuesli[2];
            STRESS[3] = stressMuesli[5];
            STRESS[4] = stressMuesli[4];
            STRESS[5] = stressMuesli[3];

            // tangent
            itensor4 tg;
            thePoint->tangentTensor(tg);

            for (unsigned i=0; i<6; i++)
            {
                for (unsigned j=0; j<6; j++)
                {
                    tangent[i][j] = tg(voigt(0,i), voigt(1,i), voigt(0,j), voigt(1,j));
                }
            }

            delete thePoint;
            delete theMaterial;

        }
            break;

        case 2:
        {
            // small strain elastoplastic material
            //
            // material constants
            //
            const double E     = PROPS[0];   // Young modulus
            const double nu    = PROPS[1];   // Poisson ratio
            const double Y0    = PROPS[2];   // yield stress
            const double Hi    = PROPS[3];   // isotropic hardening modulus
            const double Hk    = PROPS[4];   // kinematic hardening modulus
            const double rho   = 1.0;        // density

            // create material and material point
            splasticMaterial *theMaterial = new splasticMaterial("a splastic material", E, nu, rho, Hi, Hk, Y0, 0.0, "mises");
            splasticMP       *thePoint    = dynamic_cast<splasticMP*>(theMaterial->createMaterialPoint());

            // set converged states
            istensor epsn( STRAN[0], STRAN[1], STRAN[2], .5*STRAN[5], .5*STRAN[4], .5*STRAN[3]);
            istensor ep_n;
            double   dg_n;
            double   xi_n;
            istensor Xi_n;

            double tf(TIME[1]);
            if ( tf == 0.0 )
            {
                ep_n.setZero();
                dg_n = 0.0;
                xi_n = 0.0;
                Xi_n.setZero();
            }
            else
            {
                istensor ep_t( STATEV[0], STATEV[1], STATEV[2], STATEV[3], STATEV[4], STATEV[5]);
                ep_n = ep_t;
                dg_n = STATEV[6];
                xi_n = STATEV[7];
                istensor Xi_t( STATEV[8], STATEV[9], STATEV[10], STATEV[11], STATEV[12], STATEV[13]);
                Xi_n = Xi_t;
            }

            // set current states
            istensor epsc(    STRAN[0]+DSTRAN[0],      STRAN[1]+DSTRAN[1],      STRAN[2]+DSTRAN[2],
                          .5*(STRAN[5]+DSTRAN[5]), .5*(STRAN[4]+DSTRAN[4]), .5*(STRAN[3]+DSTRAN[3]));

            double timen(TIME[0]), timec(TIME[0]+*DTIME);

            thePoint->setConvergedState(timen, epsn, dg_n, ep_n, xi_n, Xi_n);
            thePoint->updateCurrentState(timec, epsc);

            // compute all the required fields
            double stressMuesli[6];
            istensor stressT;
            thePoint->stress(stressT);
            ContraContraSymTensorToVector(stressT, stressMuesli);

            STRESS[0] = stressMuesli[0];
            STRESS[1] = stressMuesli[1];
            STRESS[2] = stressMuesli[2];
            STRESS[3] = stressMuesli[5];
            STRESS[4] = stressMuesli[4];
            STRESS[5] = stressMuesli[3];

            // tangent
            itensor4 tg;
            thePoint->tangentTensor(tg);

            for (unsigned i=0; i<6; i++)
            {
                for (unsigned j=0; j<6; j++)
                {
                    tangent[i][j] = tg( voigt(0,i), voigt(1,i), voigt(0,j), voigt(1,j));
                }
            }

            // recover and store state variables
            double xi_c;
            istensor ep_c, Xi_c;

            materialState state_c = thePoint->getCurrentState();

            xi_c = state_c.theDouble[1];
            ep_c = state_c.theStensor[1];
            Xi_c = state_c.theStensor[2];

            // store the updated state variables
            STATEV[0] = ep_c(0,0);
            STATEV[1] = ep_c(1,1);
            STATEV[2] = ep_c(2,2);
            STATEV[3] = ep_c(1,2);
            STATEV[4] = ep_c(2,0);
            STATEV[5] = ep_c(0,1);
            //STATEV[6] = dg_c;
            STATEV[7] = xi_c;
            STATEV[8] = Xi_c(0,0);
            STATEV[9] = Xi_c(1,1);
            STATEV[10] = Xi_c(2,2);
            STATEV[11] = Xi_c(1,2);
            STATEV[12] = Xi_c(2,0);
            STATEV[13] = Xi_c(0,1);
            
            delete thePoint;
            delete theMaterial;
        }
            

        case 3:
        {
            // finite strain neohookean
            //
            // material constants
            //
            const double E     = PROPS[0];   // Young modulus
            const double nu    = PROPS[1];   // Poisson ratio
            const double rho   = 1.0;        // density

            // create material and material point
            neohookeanMaterial *theMaterial = new neohookeanMaterial("a neohookean material", E, nu, rho);
            neohookeanMP       *thePoint    = dynamic_cast<neohookeanMP*>(theMaterial->createMaterialPoint());

            // set converged states
            itensor Fn;
            double tf(TIME[1]);
            if (tf == 0.0)
                Fn.setZero();
            else
                muesli::vectorToTensor(DFGRD0, Fn);

            // set current state
            itensor Fc;
            muesli::vectorToTensor(DFGRD1, Fc);
            double timen(TIME[0]), timec(TIME[0]+*DTIME);

            thePoint->setConvergedState(timen, Fn);
            thePoint->updateCurrentState(timec, Fc);

            // compute Cauchy stress and pass it to Abaqus voigt
            istensor sigma;
            thePoint->CauchyStress(sigma);

            double stressMuesli[6];
            ContraContraSymTensorToVector(sigma, stressMuesli);
            STRESS[0] = stressMuesli[0];
            STRESS[1] = stressMuesli[1];
            STRESS[2] = stressMuesli[2];
            STRESS[3] = stressMuesli[5];
            STRESS[4] = stressMuesli[4];
            STRESS[5] = stressMuesli[3];

            // tangent
            itensor4 cspatial;
            thePoint->spatialTangent(cspatial);
            itensor4 cabaqus;
            abaqusFiniteStrainTangent(sigma, cspatial, cabaqus);


            for (unsigned i=0; i<6; i++)
            {
                for (unsigned j=0; j<6; j++)
                {
                    tangent[i][j] = cabaqus(voigt(0,i), voigt(1,i), voigt(0,j), voigt(1,j));
                }
            }

            // recover and store state variables
            // nothing to be done
            
            delete thePoint;
            delete theMaterial;
        }


        default:
            for (unsigned i=0; i<6; i++)
            {
                for (unsigned j=0; j<6; j++)
                {
                    tangent[i][j] = 0.0;
                }
            }
            break;
    }

    int k=0;
    for (int a=0; a< *NTENS; a++)
        for (int b=0; b< *NTENS; b++)
        {
            DDSDDE[k] = tangent[a][b];
            k++;
        }

    return 0;
}




void abaqusFiniteStrainTangent(const istensor& sigma, const itensor4& csp, itensor4& cabaqus)
{
    istensor id = istensor::identity();
    for (unsigned i=0; i<3; i++)
        for (unsigned j=0; j<3; j++)
            for (unsigned k=0; k<3; k++)
                for (unsigned l=0; l<3; l++)
                    cabaqus(i,j,k,l) = csp(i,j,k,l) + id(i,k)*sigma(j,l) + id(j,l)*sigma(i,k);
}



