Line data Source code
1 : #ifndef MODEL_SOLID_FINITE_PSEUDOLINEAR_CUBIC_H_
2 : #define MODEL_SOLID_FINITE_PSEUDOLINEAR_CUBIC_H_
3 :
4 : #include "IO/ParmParse.H"
5 : #include "Model/Solid/Solid.H"
6 :
7 : namespace Model
8 : {
9 : namespace Solid
10 : {
11 : namespace Finite
12 : {
13 : namespace PseudoLinear
14 : {
15 : class Cubic : public Solid<Set::Sym::Major>
16 : {
17 : public:
18 140 : Cubic() {};
19 : //Cubic(Solid<Set::Sym::Major> base) : Solid<Set::Sym::Major>(base) {};
20 258 : virtual ~Cubic() {};
21 :
22 960 : Set::Scalar W(const Set::Matrix & F) const override
23 : {
24 960 : q.normalize();
25 : Set::Matrix4<AMREX_SPACEDIM,Set::Sym::MajorMinor> CC
26 960 : = Set::Matrix4<AMREX_SPACEDIM,Set::Sym::MajorMinor>::Cubic(C11,C12,C44,q);
27 960 : Set::Matrix E = 0.5 * (F.transpose() * F - Set::Matrix::Identity());
28 1920 : return 0.5 * (E * (CC*E).transpose()).trace();
29 : }
30 3920 : Set::Matrix DW(const Set::Matrix & F) const override
31 : {
32 3920 : q.normalize();
33 : Set::Matrix4<AMREX_SPACEDIM,Set::Sym::MajorMinor> CC
34 3920 : = Set::Matrix4<AMREX_SPACEDIM,Set::Sym::MajorMinor>::Cubic(C11,C12,C44,q);
35 3920 : Set::Matrix C = F.transpose() * F;
36 3920 : Set::Matrix dw = Set::Matrix::Zero();
37 15020 : for (int i = 0; i < AMREX_SPACEDIM; i++)
38 43080 : for (int J = 0; J < AMREX_SPACEDIM; J++)
39 125280 : for (int A = 0; A < AMREX_SPACEDIM; A++)
40 367920 : for (int R = 0; R < AMREX_SPACEDIM; R++)
41 1087920 : for (int S = 0; S < AMREX_SPACEDIM; S++)
42 : {
43 1626600 : dw(i,J) += 0.5 * F(i,A)*CC(J,A,R,S) * (C(R,S) - ((R==S) ? 1.0 : 0.0));
44 : }
45 7840 : return dw;
46 : }
47 40 : Set::Matrix4<AMREX_SPACEDIM,Set::Sym::Major> DDW(const Set::Matrix & F) const override
48 : {
49 40 : q.normalize();
50 : Set::Matrix4<AMREX_SPACEDIM,Set::Sym::MajorMinor> CC
51 40 : = Set::Matrix4<AMREX_SPACEDIM,Set::Sym::MajorMinor>::Cubic(C11,C12,C44,q);
52 40 : Set::Matrix4<AMREX_SPACEDIM,Set::Sym::Major> ddw;
53 40 : Set::Matrix C = F.transpose() * F;
54 40 : Set::Matrix I = Set::Matrix::Identity();
55 :
56 140 : for (int i = 0; i < AMREX_SPACEDIM; i++)
57 360 : for (int J = 0; J < AMREX_SPACEDIM; J++)
58 960 : for (int k = 0; k < AMREX_SPACEDIM; k++)
59 2640 : for (int L = 0; L < AMREX_SPACEDIM; L++)
60 : {
61 1940 : ddw(i,J,k,L) = 0.0;
62 7440 : for (int A = 0; A < AMREX_SPACEDIM; A++)
63 21360 : for (int B = 0; B < AMREX_SPACEDIM; B++)
64 : {
65 31720 : ddw(i,J,k,L) += 0.5*I(i,k)*CC(J,L,A,B) * (C(A,B) - I(A,B));
66 47580 : ddw(i,J,k,L) += F(i,A)*F(k,B)*CC(A,J,B,L);
67 15860 : if (std::isnan(ddw(i,J,k,L))) Util::Abort(INFO);
68 : }
69 : }
70 :
71 80 : return ddw;
72 : }
73 0 : virtual void Print(std::ostream &out) const override
74 : {
75 : //out << "CC = " << CC;
76 0 : out << "C11=" << C11 << " C12=" << C12 << " C44=" << C44 <<
77 0 : " q=" << q.w() << " + " << q.x() << "i + " << q.y() << "j + " << q.z() << "k";
78 0 : }
79 :
80 : public:
81 : //Set::Matrix4<AMREX_SPACEDIM,Set::Sym::MajorMinor> CC;
82 : static constexpr KinematicVariable kinvar = KinematicVariable::F;
83 : Set::Scalar C11=1.68, C12=1.21, C44=0.75;
84 : mutable Set::Quaternion q = Set::Quaternion(1.0, 0.0, 0.0, 0.0);
85 :
86 : public:
87 4 : static Cubic Zero()
88 : {
89 4 : Cubic ret;
90 4 : ret.C11 = 0.0;
91 4 : ret.C12 = 0.0;
92 4 : ret.C44 = 0.0;
93 4 : ret.q = Set::Quaternion(0.0,0.0,0.0,0.0);
94 : //ret.CC = Set::Matrix4<AMREX_SPACEDIM,Set::Sym::MajorMinor>::Zero();
95 4 : return ret;
96 0 : }
97 128 : static Cubic Random()
98 : {
99 128 : return Random(Util::Random(), Util::Random(), Util::Random());
100 : }
101 128 : static Cubic Random(Set::Scalar C11, Set::Scalar C12, Set::Scalar C44)
102 : {
103 128 : Cubic ret;
104 128 : ret.C11 = C11;
105 128 : ret.C12 = C12;
106 128 : ret.C44 = C44;
107 128 : ret.q = Set::Quaternion::UnitRandom();
108 : //Set::Scalar phi1 = 2.0*Set::Constant::Pi * Util::Random();
109 : //Set::Scalar Phi = Set::Constant::Pi * Util::Random();
110 : //Set::Scalar phi2 = 2.0*Set::Constant::Pi * Util::Random();
111 : //ret.CC = Set::Matrix4<AMREX_SPACEDIM,Set::Sym::MajorMinor>::Cubic(C11,C12,C44,phi1,Phi,phi2);
112 : //if (ret.CC.contains_nan()) Util::Abort(INFO);
113 128 : return ret;
114 0 : }
115 0 : static void Parse(Cubic & value, IO::ParmParse & pp)
116 : {
117 0 : pp_query_default("C11",value.C11, 1.68); // Elastic constant
118 0 : pp_query_default("C12",value.C12, 1.21); // Elastic constant
119 0 : pp_query_default("C44",value.C44, 0.75); // Elastic constant
120 :
121 0 : if (pp.contains("random"))
122 : {
123 0 : value = Cubic::Random(value.C11,value.C12,value.C44);
124 0 : return;
125 : }
126 :
127 0 : Set::Scalar phi1 = 0.0, Phi = 0.0, phi2 = 0.0;
128 0 : pp_query_default("phi1",phi1,0.0); // Bunge Euler angle :math:`\phi_1`
129 0 : pp_query_default("Phi",Phi,0.0); // Bunge Euler angle :math:`\Phi`
130 0 : pp_query_default("phi2",phi2,0.0); // Bunge Euler angle :math:`\phi_2`
131 0 : Eigen::Matrix3d R;
132 0 : R = Eigen::AngleAxisd(phi2, Eigen::Vector3d::UnitX()) *
133 0 : Eigen::AngleAxisd(Phi, Eigen::Vector3d::UnitZ()) *
134 0 : Eigen::AngleAxisd(phi1, Eigen::Vector3d::UnitX());
135 0 : value.q = Set::Quaternion(R);
136 : }
137 :
138 : #define OP_CLASS Cubic
139 : #define OP_VARS X(C11) X(C12) X(C44) X(q)
140 : #include "Model/Solid/InClassOperators.H"
141 : };
142 : #include "Model/Solid/ExtClassOperators.H"
143 :
144 : }
145 : }
146 : }
147 : }
148 :
149 : #endif
|
