Cijkl.hh

#ifndef Cal_Cijkl_hh
#define Cal_Cijkl_hh
 
#include "../MathFunctions/transpose.hh"
#include "Material.hh"
 
namespace Dune{
    namespace Composites{
 
        Dune::FieldMatrix<double,6,6> Cijkl(Dune::Composites::Material& material){
            Dune::FieldMatrix<double,6,6> C(0.0); // Initialise Tensor as Zeros
            int type = material.Type;
 
            if (type == 0){ // Isotropic
 
                double E = material.matprop[0];
                double nu = material.matprop[1];
 
                double tmp = E / ( (1.0 + nu) * ( 1.0 - 2.0 * nu ) );
                double G = E / (2.0 * (1 + nu));
 
                C[0][0] = ( 1.0 - nu ) * tmp;    C[0][1] = nu * tmp; C[0][2] = nu * tmp;
                C[1][0] = nu * tmp; C[1][1] = ( 1.0 - nu ) * tmp; C[1][2] = nu * tmp;
                C[2][0] = nu * tmp; C[2][1] = nu * tmp; C[2][2] = ( 1.0 - nu ) * tmp;
                C[3][3] = G;
                C[4][4] = G;
                C[5][5] = G;
            }
 
            if (type == 1){ // Composite
 
                double E1 = material.matprop[0];
                double E2 = material.matprop[1];
                double E3 = material.matprop[2];
                double nu12 = material.matprop[3];
                double nu21 = nu12 * (E2 / E1);
                double nu13 = material.matprop[4];
                double nu31 = nu13 * (E3 / E1);
                double nu23 = material.matprop[5];
                double nu32 = nu23 * (E3 / E2);
                double G12 = material.matprop[6];
                double G13 = material.matprop[7];
                double G23 = material.matprop[8];
 
                Dune::FieldMatrix<double,3,3> Sinv = {{1./E1,-nu21/E2,-nu31/E3},{-nu12/E1,1./E2,-nu32/E3},{-nu13/E1,-nu23/E2,1./E3}}; 
                Sinv.invert();
 
                C[0][0] = Sinv[0][0];
                C[1][1] = Sinv[1][1];
                C[2][2] = Sinv[2][2];
                C[0][1] = Sinv[0][1]; C[1][0] = Sinv[1][0];
                C[0][2] = Sinv[0][2]; C[2][0] = Sinv[2][0];
                C[1][2] = Sinv[1][2]; C[2][1] = Sinv[2][1];
                C[3][3] = G23;
                C[4][4] = G13;
                C[5][5] = G12;
            }
            return C;
        }
    }
}
 
 
#endif