Line data    Source code

       1              : //
       2              : // This class implements basic cubic elasticity.
       3              : // For a discussion on cubic elasticity, `please see this link <http://solidmechanics.org/text/Chapter3_2/Chapter3_2.htm#Sect3_2_16>`_.
       4              : //
       5              : 
       6              : #ifndef MODEL_SOLID_LINEAR_CUBIC_H_
       7              : #define MODEL_SOLID_LINEAR_CUBIC_H_
       8              : 
       9              : #include "Model/Solid/Solid.H"
      10              : #include "IO/ParmParse.H"
      11              : 
      12              : namespace Model
      13              : {
      14              : namespace Solid
      15              : {
      16              : namespace Linear
      17              : {
      18              : class Cubic : public Solid<Set::Sym::MajorMinor>
      19              : {
      20              : public:
      21              : 
      22          100 :     Cubic() {};
      23              :     Cubic(Solid<Set::Sym::MajorMinor> base) : Solid<Set::Sym::MajorMinor>(base) {};
      24          198 :     virtual ~Cubic() {};
      25              : 
      26              :     void
      27           92 :     Define(Set::Scalar C11, Set::Scalar C12, Set::Scalar C44, Set::Scalar phi1, Set::Scalar Phi, Set::Scalar phi2)
      28              :     {
      29           92 :         ddw = Set::Matrix4<AMREX_SPACEDIM,Set::Sym::MajorMinor>::Cubic(C11,C12,C44,phi1,Phi,phi2);
      30           92 :     }
      31              :     void
      32              :     Define(Set::Scalar C11, Set::Scalar C12, Set::Scalar C44, Eigen::Matrix3d R)
      33              :     {
      34              :         ddw = Set::Matrix4<AMREX_SPACEDIM,Set::Sym::MajorMinor>::Cubic(C11,C12,C44,R);
      35              :     }
      36          520 :     Set::Scalar W(const Set::Matrix & gradu) const override
      37              :     {
      38         1040 :         return ( 0.5 * gradu.transpose() * (ddw*gradu) ).trace();
      39              :     }
      40         3920 :     Set::Matrix DW(const Set::Matrix & gradu) const override
      41              :     {
      42         7840 :         return ddw*gradu;
      43              :     }
      44           40 :     Set::Matrix4<AMREX_SPACEDIM,Set::Sym::MajorMinor> DDW(const Set::Matrix & /*gradu*/) const override
      45              :     {
      46           40 :         return ddw;
      47              :     }
      48            0 :     virtual void Print(std::ostream &out) const override 
      49              :     {
      50            0 :         out << ddw;
      51            0 :     }
      52              : 
      53              : public:
      54              :     Set::Matrix4<AMREX_SPACEDIM,Set::Sym::MajorMinor> ddw;
      55              :     static const KinematicVariable kinvar = KinematicVariable::gradu;
      56              : 
      57              :     AMREX_FORCE_INLINE
      58              :     static Cubic Combine(const std::vector<Cubic> &models, const std::vector<Set::Scalar> &eta)
      59              :     {
      60              :         Cubic ret;
      61              :         ret.ddw = Set::Matrix4<AMREX_SPACEDIM,Set::Sym::MajorMinor>::Zero();
      62              :         Set::Scalar etasum = 0.;
      63              :         for (unsigned int n = 0 ; n < models.size(); n++) etasum += eta[n];
      64              :         for (unsigned int n = 0 ; n < models.size(); n++)
      65              :         {
      66              :             ret.ddw += models[n].ddw * (eta[n] / etasum);
      67              :         }
      68              :         return ret;
      69              :     }
      70              : 
      71            4 :     static Cubic Zero()
      72              :     {
      73            4 :         Cubic ret;
      74            4 :         ret.Define(0.0,0.0,0.0,0.0,0.0,0.0);
      75            4 :         return ret;
      76              :         
      77            0 :     }
      78           44 :     static Cubic Random()
      79              :     {
      80           44 :         return Random(Util::Random(), Util::Random(), Util::Random());
      81              :     }
      82           88 :     static Cubic Random(Set::Scalar C11, Set::Scalar C12, Set::Scalar C44)
      83              :     {
      84           88 :         Cubic ret;
      85           88 :         Set::Scalar phi1 = 2.0*Set::Constant::Pi * Util::Random();
      86           88 :         Set::Scalar Phi  = 2.0*Set::Constant::Pi * Util::Random();
      87           88 :         Set::Scalar phi2 = 2.0*Set::Constant::Pi * Util::Random();
      88           88 :         ret.Define(C11,C12,C44,phi1,Phi,phi2);
      89           88 :         return ret;
      90            0 :     }
      91              : 
      92            0 :     static void Parse(Cubic & value, IO::ParmParse & pp)
      93              :     {
      94            0 :         Set::Scalar C11 = NAN, C12 = NAN, C44 = NAN;
      95            0 :         pp_query_default("C11",C11,1.68); // Elastic constant
      96            0 :         pp_query_default("C12",C12,1.21); // Elastic constant
      97            0 :         pp_query_default("C44",C44,0.75); // Elastic constant
      98              : 
      99            0 :         if (pp.contains("random"))
     100              :         {
     101            0 :             value = Cubic::Random(C11,C12,C44);
     102            0 :             return;
     103              :         }
     104              : 
     105            0 :         std::string anglefmt;
     106              :         // specify whether using radians or degrees
     107            0 :         pp_query_validate("anglefmt",anglefmt,{"radians","degrees"});
     108              : 
     109            0 :         Set::Scalar phi1 = NAN, Phi = NAN, phi2 = NAN;
     110            0 :         pp_query_default("phi1",phi1,0.0);  // Bunge Euler angle :math:`\phi_1` about x axis
     111            0 :         pp_query_default("Phi", Phi, 0.0);  // Bunge Euler angle :math:`\Phi`   about z axis
     112            0 :         pp_query_default("phi2",phi2,0.0);  // Bunge Euler angle :math:`\phi_2` about x axis
     113              : 
     114            0 :         if (anglefmt == "degrees")
     115              :         {
     116            0 :             phi1 *= Set::Constant::Pi / 180.0;
     117            0 :             Phi  *= Set::Constant::Pi / 180.0;
     118            0 :             phi2 *= Set::Constant::Pi / 180.0;
     119              :         }
     120              : 
     121            0 :         value.Define(C11,C12,C44,phi1,Phi,phi2);
     122            0 :     }
     123              : 
     124              :     #define OP_CLASS Cubic
     125              :     #define OP_VARS X(ddw)
     126              :     #include "Model/Solid/InClassOperators.H"
     127              : };
     128              : #include "Model/Solid/ExtClassOperators.H"
     129              : 
     130              : }
     131              : }
     132              : }
     133              : 
     134              : #endif