/****************************************************************************
*
*                                 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_lsdyna.h"
#include "../muesli.h"


using namespace muesli;



// matlabel:      label referring to material type
// matproperties: material properties including thermal expansion, Delta(volumetric strain)=alpha*(tempc-temp0) 
// epsilon:       strain increment since last converged solution
// time:          current time      
// sigma:         stress tensor
// history:       history variables in the last converged solution
//*****************Symmetric tensor notation correspondance between LS_DYNA and tensor_MUeSLi************************/
/****************LS_DYNA************************************************************MUesLI***************************/
//         |1 4 6|                              |                             |1 6 5|                                /
//         |  2 5|                              |                             |  2 4|                                /
//         |    3|                              |                             |    3|                                /
//*******************************************************************************************************************/

int interface_lsdyna_(float* matproperties,
                     float* epsilon,
                     float* sigma,
                     float* history,
                     float* time,
                     float* epsilon_p,
                     float* temperature,
                     int * failel
                     )
{
    int mn = static_cast<int>(matproperties[4]);
    ////////////////mn->material number/////////////
    // mn = 1, elastic small strain                /
    // mn = 2, plastic Von Mises small strain      /
    // mn = 3, viscoplastic Von Mises small strain /
    // mn = 4, Neohookean finite strain            /
    ////////////////////////////////////////////////
    switch (mn)
    {
        // linear elastic material
            //
            // matproperties[0] : E
            // matproperties[1] : nu
            // matproperties[2] : K
            // matproperties[3] : G
            // matproperties[4] : Material type(1)
            // matproperties[5] : rho
            // matproperties[6] : alpha
            // matproperties[7] : temp0 
            // history[0]       : strain_n(1,1)
            // history[1]       : strain_n(2,2)
            // history[2]       : strain_n(3,3)
            // history[3]       : strain_n(1,2)
            // history[4]       : strain_n(2,3)
            // history[5]       : strain_n(3,1)
            // epsilon[0]       : Dstrain_c(1,1)
            // epsilon[1]       : Dstrain_c(2,2)
            // epsilon[2]       : Dstrain_c(3,3)
            // epsilon[3]       : Dstrain_c(1,2)
            // epsilon[4]       : Dstrain_c(2,3)
            // epsilon[5]       : Dstrain_c(1,3)
            // sigma[0]         : sigma_c(1,1)
            // sigma[1]         : sigma_c(2,2)
            // sigma[2]         : sigma_c(3,3)
            // sigma[3]         : sigma_c(1,2)
            // sigma[4]         : sigma_c(2,3)
            // sigma[5]         : sigma_c(1,3)

        case 1:
        {
            // recover material data from interface arrays
            const double E     = matproperties[0];
            const double nu    = matproperties[1];
            const double rho   = matproperties[5];

            elasticIsotropicMaterial *theMaterial = new elasticIsotropicMaterial("an elastic material", E, nu, rho);
            elasticIsotropicMP       *thePoint    = dynamic_cast<elasticIsotropicMP*>(theMaterial->createMaterialPoint());
                   
            // set converged and current states
            istensor epsc(history[0] + epsilon[0],history[1] + epsilon[1],history[2] + epsilon[2],
                          history[4] + epsilon[4],history[5] + epsilon[5],history[3] + epsilon[3]);
          
            
            double  timec(*time);
           
            thePoint->updateCurrentState(timec, epsc);
            
            // compute all the required fields
            double stress[6];
            istensor stressT;
            thePoint->stress(stressT);
            ContraContraSymTensorToVector(stressT, stress);
            //for (unsigned k=0; k<6; k++) sigma[k] = stress[k];!!Different Voigt convention
            sigma[0] = (float)stress[0];
            sigma[1] = (float)stress[1];
            sigma[2] = (float)stress[2];
            sigma[4] = (float)stress[3];
            sigma[5] = (float)stress[4];
            sigma[3] = (float)stress[5];
            
            
            // update state in lsdyna data
            history[0] += epsilon[0];
            history[1] += epsilon[1];
            history[2] += epsilon[2];
            history[3] += epsilon[3];
            history[4] += epsilon[4];
            history[5] += epsilon[5];
            
            delete thePoint;
            delete theMaterial;
        }
            break;
            
        // small strain elastoplastic material-Von mises
        //
        // matproperties[0]  : E
        // matproperties[1]  : nu
        // matproperties[2]  : K
        // matproperties[3]  : G
        // matproperties[4]  : material type(2)
        // matproperties[5]  : rho
        // matproperties[6]  : alpha
        // matproperties[7]  : temp0
        // matproperties[8]  : Hiso
        // matproperties[9]  : Hkine
        // matproperties[10] : Yield stress
        // history[0]        : strain_n(1,1) //Total strain in state n
        // history[1]        : strain_n(2,2) //..
        // history[2]        : strain_n(3,3) //.. 
        // history[3]        : strain_n(1,2) //..
        // history[4]        : strain_n(2,3) //.. 
        // history[5]        : strain_n(3,1) //.. 
        // history[6]        : ep_n(1,1) //Plastic strain in state n 
        // history[7]        : ep_n(2,2) //..
        // history[8]        : ep_n(3,3) //..
        // history[9]        : ep_n(1,2) //..
        // history[10]       : ep_n(2,3) //..
        // history[11]       : ep_n(3,1) //..
        // history[12]       : dg_n //plastic slip in time n
        // history[13]       : xi_n //strain-like isotropic hardening variable
        // history[14]       : time_n
        // history[15]       : Xi_n(1,1) //strain-like kinematic hardening
        // history[16]       : Xi_n(2,2) //...
        // history[17]       : Xi_n(3,3) //...
        // history[18]       : Xi_n(1,2) //...
        // history[19]       : Xi_n(2,3) //...
        // history[20]       : Xi_n(3,1) //...   
        // epsilon[0]        : Dstrain_c(1,1) //strain increment
        // epsilon[1]        : Dstrain_c(2,2) //...
        // epsilon[2]        : Dstrain_c(3,3) //...
        // epsilon[3]        : Dstrain_c(1,2) //...
        // epsilon[4]        : Dstrain_c(2,3) //...
        // epsilon[5]        : Dstrain_c(1,3) //...
        // sigma[0]          : sigma_c(1,1)
        // sigma[1]          : sigma_c(2,2)
        // sigma[2]          : sigma_c(3,3)
        // sigma[3]          : sigma_c(1,2)
        // sigma[4]          : sigma_c(2,3)
        // sigma[5]          : sigma_c(1,3)
        case 2:
        {
            const double E     = matproperties[0];
            const double nu    = matproperties[1];
            const double rho   = matproperties[5];
            const double Hiso  = matproperties[8];
            const double Hkine = matproperties[9];
            const double Y0    = matproperties[10];
            const std::string  plst  = "mises";    
            splasticMaterial *theMaterial = new splasticMaterial("an splastic material", E, nu, rho, Hiso, Hkine, Y0, 0.0, plst);
            splasticMP       *thePoint    = dynamic_cast<splasticMP*>(theMaterial->createMaterialPoint());

            // set converged and current states
            istensor epsc(history[0] + epsilon[0],history[1] + epsilon[1],history[2] + epsilon[2],
                          history[4] + epsilon[4],history[5] + epsilon[5],history[3] + epsilon[3]);
            istensor epsn(history[0],history[1],history[2],
                          history[4],history[5],history[3]);
            istensor ep_n(history[6], history[7], history[8], history[10], history[11], history[9]);

            istensor Xi_n(history[15], history[16], history[17], history[19], history[20], history[18]);

            double  timen = history[14];
            double  dg_n  = history[12];
            double  xi_n  = history[13];
            double  timec(*time);
            thePoint->setConvergedState(timen, epsn, dg_n, ep_n, xi_n, Xi_n);
            thePoint->updateCurrentState(timec, epsc);
            thePoint->commitCurrentState();//no need for commitCurrentState(){Sn->Sc}, since the point will be deleted

            materialState state_n = thePoint->getCurrentState();
            xi_n = state_n.theDouble[1];
            ep_n = state_n.theStensor[1];
            Xi_n = state_n.theStensor[2];


            // compute all the required fields
            double stress[6];
            istensor stressT;
            thePoint->stress(stressT);
            ContraContraSymTensorToVector(stressT, stress);
            *epsilon_p = (float) thePoint->plasticSlip();
        
            //Stress update
            sigma[0] = (float)stress[0];
            sigma[1] = (float)stress[1];
            sigma[2] = (float)stress[2];
            sigma[4] = (float)stress[3];
            sigma[5] = (float)stress[4];
            sigma[3] = (float)stress[5];
            
            // update state in lsdyna data
            history[0] += epsilon[0];
            history[1] += epsilon[1];
            history[2] += epsilon[2];
            history[3] += epsilon[3];
            history[4] += epsilon[4];
            history[5] += epsilon[5];

            history[6]  = (float)ep_n(0,0);
            history[7]  = (float)ep_n(1,1);
            history[8]  = (float)ep_n(2,2);
            history[9]  = (float)ep_n(0,1);
            history[10] = (float)ep_n(1,2);
            history[11] = (float)ep_n(2,0);
            history[12] = (float)dg_n;
            history[13] = (float)xi_n;
            history[14] = (float)*time;
            history[15] = (float)Xi_n(0,0);
            history[16] = (float)Xi_n(1,1);
            history[17] = (float)Xi_n(2,2);
            history[18] = (float)Xi_n(0,1);
            history[19] = (float)Xi_n(1,2);
            history[20] = (float)Xi_n(2,0);
                
            delete thePoint;
            delete theMaterial;     
        }
        break;



        // small strain visco-elastoplastic material-Von mises
        //
        // matproperties[0]  : E
        // matproperties[1]  : nu
        // matproperties[2]  : K
        // matproperties[3]  : G
        // matproperties[4]  : material type(3)
        // matproperties[5]  : rho
        // matproperties[6]  : alpha
        // matproperties[7]  : temp0
        // matproperties[8]  : Hiso
        // matproperties[9]  : Hkine
        // matproperties[10] : Yield stress
        // matproperties[11] : Fuidity/Viscosity
        // history[0]        : strain_n(1,1) //Total strain in state n
        // history[1]        : strain_n(2,2) //..
        // history[2]        : strain_n(3,3) //.. 
        // history[3]        : strain_n(1,2) //..
        // history[4]        : strain_n(2,3) //.. 
        // history[5]        : strain_n(3,1) //.. 
        // history[6]        : ep_n(1,1) //Plastic strain in state n 
        // history[7]        : ep_n(2,2) //..
        // history[8]        : ep_n(3,3) //..
        // history[9]        : ep_n(1,2) //..
        // history[10]       : ep_n(2,3) //..
        // history[11]       : ep_n(3,1) //..
        // history[12]       : dg_n //plastic slip in time n
        // history[13]       : xi_n //strain-like isotropic hardening variable
        // history[14]       : time_n
        // history[15]       : Xi_n(1,1) //strain-like kinematic hardening
        // history[16]       : Xi_n(2,2) //...
        // history[17]       : Xi_n(3,3) //...
        // history[18]       : Xi_n(1,2) //...
        // history[19]       : Xi_n(2,3) //...
        // history[20]       : Xi_n(3,1) //...   
        // epsilon[0]        : Dstrain_c(1,1) //strain increment
        // epsilon[1]        : Dstrain_c(2,2) //...
        // epsilon[2]        : Dstrain_c(3,3) //...
        // epsilon[3]        : Dstrain_c(1,2) //...
        // epsilon[4]        : Dstrain_c(2,3) //...
        // epsilon[5]        : Dstrain_c(1,3) //...
        // sigma[0]          : sigma_c(1,1)
        // sigma[1]          : sigma_c(2,2)
        // sigma[2]          : sigma_c(3,3)
        // sigma[3]          : sigma_c(1,2)
        // sigma[4]          : sigma_c(2,3)
        // sigma[5]          : sigma_c(1,3)
         case 3:
        {
      
            const double E     = matproperties[0];
            const double nu    = matproperties[1];
            const double rho   = matproperties[5];
            const double Hiso  = matproperties[8];
            const double Hkine = matproperties[9];
            const double Y0    = matproperties[10];
            const double eta   = matproperties[11];
            const std::string  plst  = "mises";    

            viscoplasticMaterial *theMaterial = new viscoplasticMaterial("a viscoplastic material", E, nu, rho, Hiso, Hkine, Y0, plst, eta);
            viscoplasticMP       *thePoint    = dynamic_cast<viscoplasticMP*>(theMaterial->createMaterialPoint());

            // set converged and current states
            istensor epsc(history[0] + epsilon[0],history[1] + epsilon[1],history[2] + epsilon[2],
                          history[4] + epsilon[4],history[5] + epsilon[5],history[3] + epsilon[3]);
            istensor epsn(history[0],history[1],history[2],
                          history[4],history[5],history[3]);
            istensor ep_n(history[6], history[7], history[8], history[10], history[11], history[9]);

            istensor Xi_n(history[15], history[16], history[17], history[19], history[20], history[18]);

            double  timen = history[14];
            double  dg_n  = history[12];
            double  xi_n  = history[13];
            double  timec(*time);
            thePoint->setConvergedState(timen, dg_n, ep_n, xi_n, Xi_n, epsn);
            thePoint->updateCurrentState(timec, epsc);
            thePoint->commitCurrentState();//no need for commitCurrentState(){Sn->Sc}, since the point will be deleted

            materialState state_n = thePoint->getCurrentState();
            dg_n = state_n.theDouble[0];
            xi_n = state_n.theDouble[1];
            ep_n = state_n.theStensor[1];
            Xi_n = state_n.theStensor[2];


            // compute all the required fields
            double stress[6];
            istensor stressT;
            thePoint->stress(stressT);
            ContraContraSymTensorToVector(stressT, stress);
            *epsilon_p = (float) thePoint->plasticSlip();
        
            //Stress update
            sigma[0] = (float) stress[0];
            sigma[1] = (float) stress[1];
            sigma[2] = (float) stress[2];
            sigma[4] = (float) stress[3];
            sigma[5] = (float) stress[4];
            sigma[3] = (float) stress[5];
            
            // update state in lsdyna data
            history[0] += (float) epsilon[0];
            history[1] += (float) epsilon[1];
            history[2] += (float) epsilon[2];
            history[3] += (float) epsilon[3];
            history[4] += (float) epsilon[4];
            history[5] += (float) epsilon[5];

            history[6]  = (float) ep_n(0,0);
            history[7]  = (float) ep_n(1,1);
            history[8]  = (float) ep_n(2,2); 
            history[9]  = (float) ep_n(0,1); 
            history[10] = (float) ep_n(1,2); 
            history[11] = (float) ep_n(2,0); 
            history[12] = (float) dg_n;
            history[13] = (float) xi_n;
            history[14] = (float) *time;
            history[15] = (float) Xi_n(0,0);
            history[16] = (float) Xi_n(1,1);
            history[17] = (float) Xi_n(2,2);
            history[18] = (float) Xi_n(0,1); 
            history[19] = (float) Xi_n(1,2);
            history[20] = (float) Xi_n(2,0);
                
            delete thePoint;
            delete theMaterial;     
        }
        break;

        // Nehookean elastic material
        //
        // matproperties[0] : E
        // matproperties[1] : nu
        // matproperties[2] : K
        // matproperties[3] : G
        // matproperties[4] : Material type(4)
        // matproperties[5] : alpha
            // matproperties[6] : temp0 
            // history[0]       : F(1,1)
            // history[1]       : F(2,1)
            // history[2]       : F(3,1)
            // history[3]       : F(1,2)
            // history[4]       : F(2,2)
            // history[5]       : F(3,2)
            // history[6]       : F(1,3)
            // history[7]       : F(2,3)
            // history[8]       : F(3,3)
            // sigma[0]         : sigma_c(1,1)
            // sigma[1]         : sigma_c(2,2)
            // sigma[2]         : sigma_c(3,3)
            // sigma[3]         : sigma_c(1,2)
            // sigma[4]         : sigma_c(2,3)
            // sigma[5]         : sigma_c(1,3)
        case 4:
        {
            // recover material data from interface arrays
            const double E     = matproperties[0];
            const double nu    = matproperties[1];
            const double rho   = 0.0; //missing

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

            double   timec(*time);
            itensor  tF;
            tF(0,0) = history[0];
            tF(0,1) = history[3];
            tF(0,2) = history[6];
            tF(1,0) = history[1];
            tF(1,1) = history[4];
            tF(1,2) = history[7];
            tF(2,0) = history[2];
            tF(2,1) = history[5];
            tF(2,2) = history[8];

        
            thePoint->updateCurrentState(timec, tF);
            thePoint->commitCurrentState();
            istensor stress;
            thePoint->CauchyStress(stress);

            //!!Different Voigt convention
            sigma[0] = (float) stress(0,0);
            sigma[1] = (float) stress(1,1);
            sigma[2] = (float) stress(2,2);
            sigma[4] = (float) stress(1,2);
            sigma[5] = (float) stress(0,2);
            sigma[3] = (float) stress(0,1);
            
            delete thePoint;
            delete theMaterial;
        }
            break;


       default:
            std::cout<<"Material with number, "<<mn<<", not yet defined!!"<<std::endl;
            break;
    }

    return 0;
}
