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