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 "ComputeBlockOrientationByMisorientation.h"
      11             : #include "MooseMesh.h"
      12             : 
      13             : #include "libmesh/mesh_tools.h"
      14             : 
      15             : registerMooseObject("SolidMechanicsApp", ComputeBlockOrientationByMisorientation);
      16             : 
      17             : InputParameters
      18           8 : ComputeBlockOrientationByMisorientation::validParams()
      19             : {
      20           8 :   InputParameters params = ComputeBlockOrientationBase::validParams();
      21           8 :   params.addClassDescription("This object computes the orientation of each grain (block) by "
      22             :                              "calculating the maximum misorientation within the grain.");
      23           8 :   return params;
      24           0 : }
      25             : 
      26           4 : ComputeBlockOrientationByMisorientation::ComputeBlockOrientationByMisorientation(
      27           4 :     const InputParameters & parameters)
      28             :   : ComputeBlockOrientationBase(parameters),
      29           4 :     _updated_rotation(getMaterialProperty<RankTwoTensor>("updated_rotation")),
      30          12 :     _misorient(getMaterialProperty<Real>("misorientation"))
      31             : {
      32           4 : }
      33             : 
      34             : void
      35          17 : ComputeBlockOrientationByMisorientation::initialize()
      36             : {
      37          17 :   ComputeBlockOrientationBase::initialize();
      38             : 
      39             :   _block_ea_values.clear();
      40             :   _grain_misorientation.clear();
      41          17 : }
      42             : 
      43             : void
      44         416 : ComputeBlockOrientationByMisorientation::execute()
      45             : {
      46         416 :   RankTwoTensor rot;
      47         416 :   Real misorient = 0.0;
      48             : 
      49        3744 :   for (unsigned int _qp = 0; _qp < _qrule->n_points(); _qp++)
      50             :   {
      51        3328 :     rot += _JxW[_qp] * _coord[_qp] * _updated_rotation[_qp];
      52        3328 :     misorient += _JxW[_qp] * _coord[_qp] * _misorient[_qp];
      53             :   }
      54         416 :   rot /= _current_elem_volume;
      55         416 :   misorient /= _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             : 
      83             :   // compute EulerAngle from Quaternion
      84         416 :   EulerAngles ea = EulerAngles(q);
      85             : 
      86             :   // store value for the current misorientation in the subdomain
      87             :   // save EulerAngle in tuple so that we can gather the data from all processors
      88         416 :   _grain_misorientation[_current_elem->subdomain_id()].emplace_back(
      89             :       misorient, ea.phi1, ea.Phi, ea.phi2);
      90         416 : }
      91             : 
      92             : void
      93          17 : ComputeBlockOrientationByMisorientation::finalize()
      94             : {
      95         153 :   for (const auto & block : _fe_problem.mesh().meshSubdomains())
      96             :   {
      97         136 :     _communicator.allgather(_grain_misorientation[block]);
      98         136 :     _block_ea_values[block] = computeSubdomainEulerAngles(block);
      99             :   }
     100          17 : }
     101             : 
     102             : EulerAngles
     103         136 : ComputeBlockOrientationByMisorientation::computeSubdomainEulerAngles(const SubdomainID & sid)
     104             : {
     105             :   Real max_misorientation = 0.0; // misorientation values should within [0, pi]
     106         136 :   EulerAngles ea;
     107             :   bool has_update = false;
     108         680 :   for (const auto & [misorientation, phi1, Phi, phi2] : _grain_misorientation[sid])
     109             :   {
     110         544 :     if (misorientation > max_misorientation)
     111             :     {
     112             :       max_misorientation = misorientation;
     113         337 :       ea = EulerAngles(phi1, Phi, phi2);
     114             :       has_update = true;
     115             :     }
     116             :   }
     117             : 
     118             :   // check if we actually update EulerAngle based on misorientation
     119             :   // if not, grab it from the previous step
     120         136 :   if (!has_update)
     121           0 :     return _block_ea_values[sid];
     122             : 
     123         136 :   return ea;
     124             : }