/****************************************************************************
*
*                                 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 "tensor.h"



namespace muesli
{
    unsigned voigt(unsigned i, unsigned j)
    {
        const unsigned voigtmap[2][6] = { {0,1,2,1,2,0}, {0,1,2,2,0,1}};
        return voigtmap[i][j];
    }


    void matrixToTensor(const double Cm[6][6], itensor4& Ct)
    {
        for (unsigned i=0; i<6; i++)
        {
            Ct( voigt(0,i), voigt(1,i), voigt(0,i), voigt(1,i)) = Cm[i][i];
            Ct( voigt(1,i), voigt(0,i), voigt(0,i), voigt(1,i)) = Cm[i][i];
            Ct( voigt(0,i), voigt(1,i), voigt(1,i), voigt(0,i)) = Cm[i][i];
            Ct( voigt(1,i), voigt(0,i), voigt(1,i), voigt(0,i)) = Cm[i][i];

            for (unsigned j=i+1; j<6; j++)
            {
                Ct( voigt(0,i), voigt(1,i), voigt(0,j), voigt(1,j)) = Cm[i][j];
                Ct( voigt(1,i), voigt(0,i), voigt(0,j), voigt(1,j)) = Cm[i][j];
                Ct( voigt(0,i), voigt(1,i), voigt(1,j), voigt(0,j)) = Cm[i][j];
                Ct( voigt(1,i), voigt(0,i), voigt(1,j), voigt(0,j)) = Cm[i][j];

                Ct( voigt(0,j), voigt(1,j), voigt(0,i), voigt(1,i)) = Cm[i][j];
                Ct( voigt(1,j), voigt(0,j), voigt(0,i), voigt(1,i)) = Cm[i][j];
                Ct( voigt(0,j), voigt(1,j), voigt(1,i), voigt(0,i)) = Cm[i][j];
                Ct( voigt(1,j), voigt(0,j), voigt(1,i), voigt(0,i)) = Cm[i][j];
            }
        }
    }

    
    void ContraContraSymTensorToVector(const istensor& St, double Sv[6])
    {
        for (unsigned i=0; i<6; i++)
            Sv[i] = St( voigt(0,i),voigt(1,i));
    }



    void CovaCovaSymTensorToVector(const istensor& St, double Sv[6])
    {
        for (unsigned i=0; i<3; i++)
            Sv[i] = St( voigt(0,i),voigt(1,i));

        for (unsigned i=3; i<6; i++)
            Sv[i] = 2.0*St( voigt(0,i),voigt(1,i));
    }



    void tensorToVector(const itensor& Ft,  double Fv[9])
    {
        for (unsigned i=0; i<3; i++)
            for (unsigned j=0; j<3; j++)
            {
                Fv[3*i+j] = Ft(i,j);
            }
    }



    void tensorToMatrix(const itensor4& Ct, double Cm[6][6])
    {
        for (unsigned i=0; i<6; i++)
        {
            Cm[i][i] = Ct( voigt(0,i), voigt(1,i), voigt(0,i), voigt(1,i));
            for (unsigned j=i+1; j<6; j++)
            {
                Cm[i][j] = Cm[j][i] = Ct( voigt(0,i), voigt(1,i), voigt(0,j), voigt(1,j));
            }
        }
    }


    
    void vectorToContraContraSymTensor(const double Sv[6], istensor& St)
    {
        for (unsigned i=0; i<6; i++)
        {
            St( voigt(0,i),voigt(1,i)) = St( voigt(1,i),voigt(0,i)) = Sv[i];
        }
    }
    
    
    
    void vectorToCovaCovaSymTensor(const double Sv[6], istensor& St)
    {
        for (unsigned i=0; i<3; i++)
            St( voigt(0,i),voigt(1,i)) = St( voigt(1,i),voigt(0,i)) = Sv[i];

        for (unsigned i=3; i<6; i++)
            St( voigt(0,i),voigt(1,i)) = St( voigt(1,i),voigt(0,i)) = 0.5*Sv[i];
    }



    void vectorToTensor(const double Fv[9], itensor& Ft)
    {
        for (unsigned i=0; i<3; i++)
            for (unsigned j=0; j<3; j++)
            {
                Ft(i,j) = Fv[3*i+j];
            }
    }
}

