Line data Source code
1 : #ifndef MODEL_INTERFACE_GB_H
2 : #define MODEL_INTERFACE_GB_H
3 :
4 : #include <AMReX.H>
5 : #include <AMReX_AmrCore.H>
6 :
7 : #include <eigen3/Eigen/Eigenvalues>
8 :
9 : #include <iostream>
10 : #include <fstream>
11 :
12 : namespace Model
13 : {
14 : namespace Interface
15 : {
16 : namespace GB
17 : {
18 : class GB
19 : {
20 : public:
21 12224 : GB() {};
22 12224 : virtual ~GB() {};
23 : virtual Set::Scalar W(const Set::Scalar theta) const = 0;
24 : virtual Set::Scalar DW(const Set::Scalar theta) const = 0;
25 : virtual Set::Scalar DDW(const Set::Scalar theta) const = 0;
26 1362712 : virtual Set::Scalar W(const Set::Vector& a_n) const { return W(atan2(a_n(1), a_n(0))); };
27 0 : virtual Set::Scalar DW(const Set::Vector&, const Set::Vector&) const { return NAN; };
28 0 : virtual Set::Scalar DDW(const Set::Vector&, const Set::Vector&) const { return NAN; };
29 :
30 : void ExportToFile(std::string filename, amrex::Real dTheta)
31 : {
32 : std::ofstream outFile;
33 : outFile.open(filename);
34 :
35 : for (amrex::Real theta = 0; theta < 2 * pi; theta = theta + dTheta)
36 : {
37 : outFile << theta << " " << W(theta) << std::endl;
38 : }
39 : outFile.close();
40 : }
41 :
42 : enum Regularization { Wilhelm, K23 };
43 :
44 : // Return anisotropic driving force, regularization
45 : std::tuple<Set::Scalar, Set::Scalar>
46 1362712 : DrivingForce(const Set::Vector& Deta, const Set::Matrix& DDeta, const Set::Matrix4<AMREX_SPACEDIM, Set::Sym::Full>& DDDDeta)
47 : {
48 : #if AMREX_SPACEDIM == 2
49 1362712 : Set::Scalar Theta = atan2(Deta(1), Deta(0));
50 1362712 : Set::Scalar sigma = /*pf.l_gb*0.75**/W(Theta);
51 1362712 : Set::Scalar Dsigma = /*pf.l_gb*0.75**/DW(Theta);
52 1362712 : Set::Scalar DDsigma = /*pf.l_gb*0.75**/DDW(Theta);
53 1362712 : Set::Scalar sinTheta = sin(Theta);
54 1362712 : Set::Scalar cosTheta = cos(Theta);
55 1362712 : Set::Scalar sin2Theta = sinTheta * sinTheta;
56 1362712 : Set::Scalar cos2Theta = cosTheta * cosTheta;
57 1362712 : Set::Scalar cosThetasinTheta = cosTheta * sinTheta;
58 :
59 : Set::Scalar boundary_term =
60 1362712 : -(sigma * DDeta.trace() +
61 1362712 : Dsigma * (cos(2.0 * Theta) * DDeta(0, 1) + 0.5 * sin(2.0 * Theta) * (DDeta(1, 1) - DDeta(0, 0))) +
62 1362712 : 0.5 * DDsigma * (sin2Theta * DDeta(0, 0) - 2. * cosThetasinTheta * DDeta(0, 1) + cos2Theta * DDeta(1, 1)));
63 :
64 : Set::Scalar curvature_term =
65 1362712 : DDDDeta(0, 0, 0, 0) * (sin2Theta * sin2Theta) +
66 1362712 : DDDDeta(0, 0, 0, 1) * (4.0 * sin2Theta * cosThetasinTheta) +
67 1362712 : DDDDeta(0, 0, 1, 1) * (6.0 * sin2Theta * cos2Theta) +
68 1362712 : DDDDeta(0, 1, 1, 1) * (4.0 * cosThetasinTheta * cos2Theta) +
69 1362712 : DDDDeta(1, 1, 1, 1) * (cos2Theta * cos2Theta);
70 2725424 : return std::make_tuple(boundary_term, curvature_term);
71 :
72 : #elif AMREX_SPACEDIM == 3
73 :
74 0 : Set::Vector normal = Deta / Deta.lpNorm<2>();
75 :
76 : // GRAHM-SCHMIDT PROCESS to get orthonormal basis
77 0 : const Set::Vector e1(1, 0, 0), e2(0, 1, 0), e3(0, 0, 1);
78 0 : Set::Vector _t2, _t3;
79 0 : if (fabs(normal(0)) > fabs(normal(1)) && fabs(normal(0)) > fabs(normal(2)))
80 : {
81 0 : _t2 = e2 - normal.dot(e2) * normal; _t2 /= _t2.lpNorm<2>();
82 0 : _t3 = e3 - normal.dot(e3) * normal - _t2.dot(e3) * _t2; _t3 /= _t3.lpNorm<2>();
83 : }
84 0 : else if (fabs(normal(1)) > fabs(normal(0)) && fabs(normal(1)) > fabs(normal(2)))
85 : {
86 0 : _t2 = e1 - normal.dot(e1) * normal; _t2 /= _t2.lpNorm<2>();
87 0 : _t3 = e3 - normal.dot(e3) * normal - _t2.dot(e3) * _t2; _t3 /= _t3.lpNorm<2>();
88 : }
89 : else
90 : {
91 0 : _t2 = e1 - normal.dot(e1) * normal; _t2 /= _t2.lpNorm<2>();
92 0 : _t3 = e2 - normal.dot(e2) * normal - _t2.dot(e2) * _t2; _t3 /= _t3.lpNorm<2>();
93 : }
94 :
95 : // Compute Hessian projected into tangent space (spanned by _t1,_t2)
96 0 : Eigen::Matrix2d DDeta2D;
97 0 : DDeta2D <<
98 0 : _t2.dot(DDeta * _t2), _t2.dot(DDeta * _t3),
99 0 : _t3.dot(DDeta * _t2), _t3.dot(DDeta * _t3);
100 0 : Eigen::SelfAdjointEigenSolver<Eigen::Matrix2d> eigensolver(2);
101 0 : eigensolver.computeDirect(DDeta2D);
102 0 : Eigen::Matrix2d eigenvecs = eigensolver.eigenvectors();
103 :
104 : // Compute tangent vectors embedded in R^3
105 0 : Set::Vector t2 = _t2 * eigenvecs(0, 0) + _t3 * eigenvecs(0, 1),
106 0 : t3 = _t2 * eigenvecs(1, 0) + _t3 * eigenvecs(1, 1);
107 :
108 : // Compute components of second Hessian in t2,t3 directions
109 0 : Set::Scalar DH2 = 0.0, DH3 = 0.0;
110 0 : Set::Scalar DH23 = 0.0;
111 0 : for (int p = 0; p < 3; p++)
112 0 : for (int q = 0; q < 3; q++)
113 0 : for (int r = 0; r < 3; r++)
114 0 : for (int s = 0; s < 3; s++)
115 : {
116 0 : DH2 += DDDDeta(p, q, r, s) * t2(p) * t2(q) * t2(r) * t2(s);
117 0 : DH3 += DDDDeta(p, q, r, s) * t3(p) * t3(q) * t3(r) * t3(s);
118 0 : DH23 += DDDDeta(p, q, r, s) * t2(p) * t2(q) * t3(r) * t3(s);
119 : }
120 :
121 0 : Set::Scalar sigma = W(normal);
122 0 : Set::Scalar DDK2 = DDW(normal, _t2);
123 0 : Set::Scalar DDK3 = DDW(normal, _t3);
124 :
125 : // GB energy anisotropy term
126 0 : Set::Scalar gbenergy_df = -sigma * DDeta.trace() - DDK2 * DDeta2D(0, 0) - DDK3 * DDeta2D(1, 1);
127 :
128 : // Second order curvature term
129 0 : Set::Scalar reg_df = NAN;
130 0 : if (regularization == Regularization::Wilhelm) reg_df = DH2 + DH3 + 2.0 * DH23;
131 0 : else if (regularization == Regularization::K23) reg_df = DH2 + DH3;
132 :
133 0 : return std::make_tuple(gbenergy_df, reg_df);
134 :
135 : #endif
136 : }
137 :
138 :
139 : protected:
140 : static constexpr amrex::Real pi = 3.14159265359;
141 : Regularization regularization = Regularization::Wilhelm;
142 : };
143 : }
144 : }
145 : }
146 :
147 : #endif
|
