Line data Source code
1 : //* This file is part of the MOOSE framework
2 : //* https://www.mooseframework.org
3 : //*
4 : //* All rights reserved, see COPYRIGHT for full restrictions
5 : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 : //*
7 : //* Licensed under LGPL 2.1, please see LICENSE for details
8 : //* https://www.gnu.org/licenses/lgpl-2.1.html
9 :
10 : #include "ComputeBlockOrientationByRotation.h"
11 : #include "MooseMesh.h"
12 :
13 : #include "libmesh/mesh_tools.h"
14 :
15 : registerMooseObject("SolidMechanicsApp", ComputeBlockOrientationByRotation);
16 :
17 : InputParameters
18 8 : ComputeBlockOrientationByRotation::validParams()
19 : {
20 8 : InputParameters params = ComputeBlockOrientationBase::validParams();
21 16 : params.addParam<unsigned int>("bins", 20, "Number of bins to segregate quaternions");
22 16 : params.addParam<Real>("L_norm", 1, "Specifies the type of average the user intends to perform");
23 8 : params.addClassDescription(
24 : "This object determines the orientation of each grain (block) by identifying the most common "
25 : "orientation direction among the material points within the grain.");
26 8 : return params;
27 0 : }
28 :
29 4 : ComputeBlockOrientationByRotation::ComputeBlockOrientationByRotation(
30 4 : const InputParameters & parameters)
31 : : ComputeBlockOrientationBase(parameters),
32 4 : _updated_rotation(getMaterialProperty<RankTwoTensor>("updated_rotation")),
33 8 : _bins(getParam<unsigned int>("bins")),
34 12 : _L_norm(getParam<Real>("L_norm"))
35 : {
36 4 : }
37 :
38 : void
39 17 : ComputeBlockOrientationByRotation::initialize()
40 : {
41 17 : ComputeBlockOrientationBase::initialize();
42 :
43 : _block_ea_values.clear();
44 : _quat.clear();
45 17 : }
46 :
47 : void
48 416 : ComputeBlockOrientationByRotation::execute()
49 : {
50 416 : RankTwoTensor rot;
51 416 : MathUtils::mooseSetToZero(rot);
52 :
53 3744 : for (unsigned int _qp = 0; _qp < _qrule->n_points(); _qp++)
54 3328 : rot += _JxW[_qp] * _coord[_qp] * _updated_rotation[_qp];
55 416 : rot /= _current_elem_volume;
56 :
57 : // transform RankTwoTensor to Eigen::Matrix
58 : Eigen::Matrix<Real, 3, 3> rot_mat;
59 :
60 1664 : for (unsigned int i = 0; i < 3; ++i)
61 4992 : for (unsigned int j = 0; j < 3; ++j)
62 3744 : rot_mat(i, j) = rot(i, j);
63 :
64 : // SVD-based projection to SO(3)
65 : Eigen::JacobiSVD<Eigen::Matrix<Real, 3, 3>> svd(rot_mat,
66 416 : Eigen::ComputeFullU | Eigen::ComputeFullV);
67 :
68 : Eigen::Matrix<Real, 3, 3> U = svd.matrixU();
69 : Eigen::Matrix<Real, 3, 3> V = svd.matrixV();
70 :
71 : // Enforce proper rotation (det = +1)
72 416 : Eigen::Matrix<Real, 3, 3> R = U * V.transpose();
73 416 : if (R.determinant() < 0.0)
74 : {
75 0 : Eigen::Matrix<Real, 3, 3> D = Eigen::Matrix<Real, 3, 3>::Identity();
76 0 : D(2, 2) = -1.0;
77 0 : R = U * D * V.transpose();
78 : }
79 :
80 : // compute Quaternion from rotation matrix
81 : Eigen::Quaternion<Real> q(R);
82 : // store quaternion values for the current element
83 832 : _quat[_current_elem->subdomain_id()].push_back(std::make_tuple(q.w(), q.x(), q.y(), q.z()));
84 416 : }
85 :
86 : void
87 17 : ComputeBlockOrientationByRotation::finalize()
88 : {
89 153 : for (const auto & block : _fe_problem.mesh().meshSubdomains())
90 : {
91 : // Use allgather to sync data from all processors
92 136 : _communicator.allgather(_quat[block]);
93 136 : _block_ea_values[block] = computeSubdomainEulerAngles(block);
94 : }
95 17 : }
96 :
97 : EulerAngles
98 136 : ComputeBlockOrientationByRotation::computeSubdomainEulerAngles(const SubdomainID & sid)
99 : {
100 : // creating a map to store the quaternion count for each bin index
101 : std::map<std::tuple<int, int, int, int>, unsigned int> feature_weights;
102 :
103 : // precompute the scaling from the bins once and reuse
104 136 : const Real bin_scale = 0.5 * _bins;
105 :
106 680 : for (const auto & [w, x, y, z] : _quat[sid])
107 : {
108 544 : const auto bin = std::make_tuple<int, int, int, int>(std::floor(w * bin_scale),
109 544 : std::floor(x * bin_scale),
110 544 : std::floor(y * bin_scale),
111 544 : std::floor(z * bin_scale));
112 544 : feature_weights[bin]++;
113 : }
114 :
115 : // quaternion average matrix Q*w*Q^T
116 : typedef Eigen::Matrix<Real, 4, 4> Matrix4x4;
117 : Matrix4x4 quat_mat = Matrix4x4::Zero();
118 : typedef Eigen::Matrix<Real, 4, 1> Vector4;
119 :
120 : Real total_weight = 0.0;
121 :
122 680 : for (const auto & [w, x, y, z] : _quat[sid])
123 : {
124 : Vector4 v(w, x, y, z);
125 :
126 544 : const auto bin = std::make_tuple<int, int, int, int>(std::floor(w * bin_scale),
127 544 : std::floor(x * bin_scale),
128 544 : std::floor(y * bin_scale),
129 544 : std::floor(z * bin_scale));
130 544 : const auto bin_size = feature_weights[bin];
131 544 : const auto weight = std::pow(bin_size, _L_norm);
132 544 : total_weight += weight;
133 :
134 544 : quat_mat += v * weight * v.transpose();
135 : }
136 :
137 136 : quat_mat *= 1.0 / total_weight;
138 :
139 : // compute eigenvalues and eigenvectors
140 136 : Eigen::EigenSolver<Matrix4x4> EigenSolver(quat_mat);
141 136 : Vector4 eigen_values = EigenSolver.eigenvalues().real();
142 136 : Matrix4x4 eigen_vectors = EigenSolver.eigenvectors().real();
143 :
144 : // Selecting eigenvector corresponding to max eigenvalue to compute average Euler angle
145 136 : Vector4::Index max_index = 0;
146 : eigen_values.maxCoeff(&max_index);
147 136 : const auto max_vec = eigen_vectors.col(max_index) / eigen_vectors.col(max_index).norm();
148 : const Eigen::Quaternion<Real> q(max_vec(0), max_vec(1), max_vec(2), max_vec(3));
149 :
150 272 : return EulerAngles(q);
151 : }
|
