Line data Source code
1 : //
2 : // This class implements elasticity for a linear material with hexagonal symmetry.
3 : // The elastic modulus tensor :math:`\mathbb{C}` has major and minor symmetry but not cubic symmetry,
4 : // so the symmetry class for this material is :ref:`Set::Sym::MajorMinor`.
5 : // The five elastic constants in the unrotated system are :math:`C_{11},C_{12},C_{13},C_{33},C_{44}`,
6 : // and the resulting elastic modulus tensor is
7 : //
8 : // .. math::
9 : //
10 : // \mathbb{C} &=
11 : // \begin{bmatrix}
12 : // C_{0000} & C_{0011} & C_{0022} & C_{0012} & C_{0020} & C_{0001} \\
13 : // C_{1100} & C_{1111} & C_{1122} & C_{1112} & C_{1120} & C_{1101} \\
14 : // C_{2200} & C_{2211} & C_{2222} & C_{2212} & C_{2220} & C_{2201} \\
15 : // C_{1200} & C_{1211} & C_{1222} & C_{1212} & C_{1220} & C_{1201} \\
16 : // C_{2000} & C_{2011} & C_{2022} & C_{2012} & C_{2020} & C_{2001} \\
17 : // C_{0100} & C_{0111} & C_{0122} & C_{0112} & C_{0120} & C_{0101}
18 : // \end{bmatrix} \\
19 : // &=
20 : // \begin{bmatrix}
21 : // C_{11} & C_{12} & C_{12} & & & \\
22 : // C_{12} & C_{11} & C_{13} & & & \\
23 : // C_{12} & C_{13} & C_{33} & & & \\
24 : // & & & C_{44} & & \\
25 : // & & & & 0.5(C_{11}-C_{12}) & \\
26 : // & & & & & 0.5(C_{11}-C_{12})
27 : // \end{bmatrix}
28 : //
29 : // Rotations can be specified using Bunge euler angles to generate a rotation matrix :math:`\mathbf{R}`.
30 : // The rotated modulus tensor :math:`\mathbb{C}^{rot}` is then
31 : //
32 : // .. math::
33 : //
34 : // \mathbb{C}^{rot}_{psqt} = \mathbb{C}_{ijkl}R_{pi}R_{qj}R_{sk}R_{tl}
35 : //
36 :
37 : #ifndef MODEL_SOLID_LINEAR_HEXAGONAL_H_
38 : #define MODEL_SOLID_LINEAR_HEXAGONAL_H_
39 :
40 : #include "Model/Solid/Solid.H"
41 : #include "IO/ParmParse.H"
42 :
43 : namespace Model
44 : {
45 : namespace Solid
46 : {
47 : namespace Linear
48 : {
49 : class Hexagonal : public Solid<Set::Sym::MajorMinor>
50 : {
51 : public:
52 :
53 1081 : Hexagonal() {};
54 : Hexagonal(Solid<Set::Sym::MajorMinor> base) : Solid<Set::Sym::MajorMinor>(base) {};
55 1432 : virtual ~Hexagonal() {};
56 :
57 : void
58 219 : Define(Set::Scalar C11, Set::Scalar C12, Set::Scalar C13, Set::Scalar C33, Set::Scalar C44, Set::Scalar phi1, Set::Scalar Phi, Set::Scalar phi2)
59 : {
60 219 : Eigen::Matrix3d m;
61 219 : m = Eigen::AngleAxisd(phi2, Eigen::Vector3d::UnitX()) *
62 438 : Eigen::AngleAxisd(Phi, Eigen::Vector3d::UnitZ()) *
63 438 : Eigen::AngleAxisd(phi1, Eigen::Vector3d::UnitX());
64 219 : Define(C11,C12,C13,C33,C44,m);
65 219 : }
66 : void
67 219 : Define(Set::Scalar C11, Set::Scalar C12, Set::Scalar C13, Set::Scalar C33, Set::Scalar C44, Eigen::Matrix3d R)
68 : {
69 :
70 : amrex::Real Ctmp[3][3][3][3];
71 219 : ddw = Set::Matrix4<AMREX_SPACEDIM,Set::Sym::MajorMinor>::Zero();
72 :
73 876 : for(int i = 0; i < 3; i++)
74 2628 : for(int j = 0; j < 3; j++)
75 7884 : for(int k = 0; k < 3; k++)
76 23652 : for(int l = 0; l < 3; l++)
77 : {
78 17739 : if(i == j && j == k && k == l)
79 : {
80 657 : if (i==0) Ctmp[i][j][k][l] = C11;
81 438 : else if (i==1) Ctmp[i][j][k][l] = C11;
82 219 : else if (i==2) Ctmp[i][j][k][l] = C33;
83 : }
84 17082 : else if (i==k && j==l)
85 : {
86 1314 : if ((i==0 && j==1) || (i==1 && j==0))
87 438 : Ctmp[i][j][k][l] = 0.5*(C11-C12);
88 : else
89 876 : Ctmp[i][j][k][l] = C44;
90 : }
91 15768 : else if (i==j && k==l)
92 : {
93 1314 : if ((i==0 && k==1) || (i==1 && k==0)) Ctmp[i][j][k][l] = C12;
94 876 : else if ((i==0 && k==2) || (i==2 && k==0)) Ctmp[i][j][k][l] = C13;
95 438 : else if ((i==1 && k==2) || (i==2 && k==1)) Ctmp[i][j][k][l] = C13;
96 : }
97 14454 : else Ctmp[i][j][k][l] = 0.0;
98 : }
99 830 : for(int p = 0; p < AMREX_SPACEDIM; p++)
100 2352 : for(int q = 0; q < AMREX_SPACEDIM; q++)
101 6780 : for(int s = 0; s < AMREX_SPACEDIM; s++)
102 19788 : for(int t = 0; t < AMREX_SPACEDIM; t++)
103 : {
104 14749 : ddw(p,q,s,t) = 0.0;
105 58996 : for(int i = 0; i < 3; i++)
106 176988 : for(int j = 0; j < 3; j++)
107 530964 : for(int k = 0; k < 3; k++)
108 1592888 : for(int l = 0; l < 3; l++)
109 2389342 : ddw(p,q,s,t) += R(p,i)*R(s,k)*Ctmp[i][j][k][l]*R(q,j)*R(t,l);
110 : }
111 219 : }
112 520 : Set::Scalar W(const Set::Matrix & gradu) const override
113 : {
114 1040 : return ( 0.5 * gradu.transpose() * (ddw*gradu) ).trace();
115 : }
116 28920 : Set::Matrix DW(const Set::Matrix & gradu) const override
117 : {
118 57840 : return ddw*gradu;
119 : }
120 12540 : Set::Matrix4<AMREX_SPACEDIM,Set::Sym::MajorMinor> DDW(const Set::Matrix & /*gradu*/) const override
121 : {
122 12540 : return ddw;
123 : }
124 0 : virtual void Print(std::ostream &out) const override
125 : {
126 0 : out << ddw;
127 0 : }
128 :
129 : public:
130 : Set::Matrix4<AMREX_SPACEDIM,Set::Sym::MajorMinor> ddw;
131 : static const KinematicVariable kinvar = KinematicVariable::gradu;
132 :
133 : AMREX_FORCE_INLINE
134 : static Hexagonal Combine(const std::vector<Hexagonal> &models, const std::vector<Set::Scalar> &eta)
135 : {
136 : Hexagonal ret;
137 : ret.ddw = Set::Matrix4<AMREX_SPACEDIM,Set::Sym::MajorMinor>::Zero();
138 : Set::Scalar etasum = 0.;
139 : for (unsigned int n = 0 ; n < models.size(); n++) etasum += eta[n];
140 : for (unsigned int n = 0 ; n < models.size(); n++)
141 : {
142 : ret.ddw += models[n].ddw * (eta[n] / etasum);
143 : }
144 : return ret;
145 : }
146 :
147 130 : static Hexagonal Zero()
148 : {
149 130 : Hexagonal ret;
150 130 : ret.Define(0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0);
151 130 : return ret;
152 :
153 : }
154 44 : static Hexagonal Random()
155 : {
156 44 : return Random(Util::Random(), Util::Random(), Util::Random(),Util::Random(),Util::Random());
157 : }
158 88 : static Hexagonal Random(Set::Scalar C11, Set::Scalar C12, Set::Scalar C13, Set::Scalar C33, Set::Scalar C44)
159 : {
160 88 : Hexagonal ret;
161 88 : Set::Scalar phi1 = 2.0*Set::Constant::Pi * Util::Random();
162 88 : Set::Scalar Phi = 2.0*Set::Constant::Pi * Util::Random();
163 88 : Set::Scalar phi2 = 2.0*Set::Constant::Pi * Util::Random();
164 88 : ret.Define(C11,C12,C13,C33,C44,phi1,Phi,phi2);
165 88 : return ret;
166 : }
167 :
168 1 : static void Parse(Hexagonal & value, IO::ParmParse & pp)
169 : {
170 1 : Set::Scalar C11 = NAN, C12 = NAN, C13 = NAN, C33 = NAN, C44 = NAN;
171 1 : pp_query_required("C11",C11); // Elastic constant
172 1 : pp_query_required("C12",C12); // Elastic constant
173 1 : pp_query_required("C13",C13); // Elastic constant
174 1 : pp_query_required("C33",C33); // Elastic constant
175 1 : pp_query_required("C44",C44); // Elastic constant
176 :
177 1 : if (pp.contains("random"))
178 : {
179 0 : value = Hexagonal::Random(C11,C12,C13,C33,C44);
180 0 : return;
181 : }
182 :
183 1 : Set::Scalar phi1 = NAN, Phi = NAN, phi2 = NAN;
184 1 : pp_query_default("phi1",phi1,0.0); // Bunge Euler angle :math:`\phi_1`
185 1 : pp_query_default("Phi", Phi, 0.0); // Bunge Euler angle :math:`\Phi`
186 1 : pp_query_default("phi2",phi2,0.0); // Bunge Euler angle :math:`\phi_2`
187 1 : value.Define(C11,C12,C13,C33,C44,phi1,Phi,phi2);
188 : }
189 :
190 : #define OP_CLASS Hexagonal
191 : #define OP_VARS X(ddw)
192 : #include "Model/Solid/InClassOperators.H"
193 : };
194 : #include "Model/Solid/ExtClassOperators.H"
195 :
196 : }
197 : }
198 : }
199 :
200 : #endif
|
