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         103 :     Cubic() {};
      23             :     Cubic(Solid<Set::Sym::MajorMinor> base) : Solid<Set::Sym::MajorMinor>(base) {};
      24         218 :     virtual ~Cubic() {};
      25             : 
      26             :     void
      27          95 :     Define(Set::Scalar C11, Set::Scalar C12, Set::Scalar C44, Set::Scalar phi1, Set::Scalar Phi, Set::Scalar phi2)
      28             :     {
      29          95 :         ddw = Set::Matrix4<AMREX_SPACEDIM,Set::Sym::MajorMinor>::Cubic(C11,C12,C44,phi1,Phi,phi2);
      30          95 :     }
      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           7 :     static Cubic Zero()
      72             :     {
      73           7 :         Cubic ret;
      74           7 :         ret.Define(0.0,0.0,0.0,0.0,0.0,0.0);
      75           7 :         return ret;
      76             :         
      77             :     }
      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             :     }
      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             :     }
     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