Line data    Source code

       1             : //* This file is part of the MOOSE framework
       2             : //* https://mooseframework.inl.gov
       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 "ComputeGrainForceAndTorque.h"
      11             : #include "GrainTrackerInterface.h"
      12             : 
      13             : #include "libmesh/quadrature.h"
      14             : 
      15             : registerMooseObject("PhaseFieldApp", ComputeGrainForceAndTorque);
      16             : 
      17             : InputParameters
      18          60 : ComputeGrainForceAndTorque::validParams()
      19             : {
      20          60 :   InputParameters params = ShapeElementUserObject::validParams();
      21          60 :   params.addClassDescription("Userobject for calculating force and torque acting on a grain");
      22         120 :   params.addParam<MaterialPropertyName>("force_density", "force_density", "Force density material");
      23         120 :   params.addParam<UserObjectName>("grain_data", "center of mass of grains");
      24         120 :   params.addCoupledVar("c", "Concentration field");
      25         120 :   params.addCoupledVar("etas", "Array of coupled order parameters");
      26          60 :   return params;
      27           0 : }
      28             : 
      29          32 : ComputeGrainForceAndTorque::ComputeGrainForceAndTorque(const InputParameters & parameters)
      30             :   : DerivativeMaterialInterface<ShapeElementUserObject>(parameters),
      31             :     GrainForceAndTorqueInterface(),
      32          32 :     _c_name(coupledName("c", 0)),
      33          32 :     _c_var(coupled("c")),
      34          32 :     _dF_name(getParam<MaterialPropertyName>("force_density")),
      35          32 :     _dF(getMaterialPropertyByName<std::vector<RealGradient>>(_dF_name)),
      36          32 :     _dFdc(getMaterialPropertyByName<std::vector<RealGradient>>(
      37          32 :         derivativePropertyNameFirst(_dF_name, _c_name))),
      38          32 :     _op_num(coupledComponents("etas")),
      39          32 :     _grain_tracker(getUserObject<GrainTrackerInterface>("grain_data")),
      40          32 :     _vals_var(_op_num),
      41          32 :     _vals_name(_op_num),
      42          64 :     _dFdgradeta(_op_num)
      43             : {
      44          96 :   for (unsigned int i = 0; i < _op_num; ++i)
      45             :   {
      46          64 :     _vals_var[i] = coupled("etas", i);
      47         128 :     _vals_name[i] = coupledName("etas", i);
      48          64 :     _dFdgradeta[i] = &getMaterialPropertyByName<std::vector<Real>>(
      49         128 :         derivativePropertyNameFirst(_dF_name, _vals_name[i]));
      50             :   }
      51          32 : }
      52             : 
      53             : void
      54         416 : ComputeGrainForceAndTorque::initialize()
      55             : {
      56         416 :   _grain_num = _grain_tracker.getTotalFeatureCount();
      57         416 :   _ncomp = 6 * _grain_num;
      58             : 
      59         416 :   _force_values.resize(_grain_num);
      60         416 :   _torque_values.resize(_grain_num);
      61         416 :   _force_torque_store.assign(_ncomp, 0.0);
      62             : 
      63         416 :   if (_fe_problem.currentlyComputingJacobian())
      64             :   {
      65         181 :     _total_dofs = _subproblem.es().n_dofs();
      66         181 :     _force_torque_c_jacobian_store.assign(_ncomp * _total_dofs, 0.0);
      67         181 :     _force_torque_eta_jacobian_store.resize(_op_num);
      68             : 
      69         543 :     for (unsigned int i = 0; i < _op_num; ++i)
      70         362 :       _force_torque_eta_jacobian_store[i].assign(_ncomp * _total_dofs, 0.0);
      71             :   }
      72         416 : }
      73             : 
      74             : void
      75       47755 : ComputeGrainForceAndTorque::execute()
      76             : {
      77       47755 :   const auto & op_to_grains = _grain_tracker.getVarToFeatureVector(_current_elem->id());
      78             : 
      79      141265 :   for (unsigned int i = 0; i < _grain_num; ++i)
      80      280530 :     for (unsigned int j = 0; j < _op_num; ++j)
      81      187020 :       if (i == op_to_grains[j])
      82             :       {
      83       28016 :         const auto centroid = _grain_tracker.getGrainCentroid(i);
      84      140080 :         for (_qp = 0; _qp < _qrule->n_points(); ++_qp)
      85      112064 :           if (_dF[_qp][j](0) != 0.0 || _dF[_qp][j](1) != 0.0 || _dF[_qp][j](2) != 0.0)
      86             :           {
      87             :             const RealGradient compute_torque =
      88        2010 :                 _JxW[_qp] * _coord[_qp] *
      89        2010 :                 (_current_elem->vertex_average() - centroid).cross(_dF[_qp][j]);
      90        2010 :             _force_torque_store[6 * i + 0] += _JxW[_qp] * _coord[_qp] * _dF[_qp][j](0);
      91        2010 :             _force_torque_store[6 * i + 1] += _JxW[_qp] * _coord[_qp] * _dF[_qp][j](1);
      92        2010 :             _force_torque_store[6 * i + 2] += _JxW[_qp] * _coord[_qp] * _dF[_qp][j](2);
      93        2010 :             _force_torque_store[6 * i + 3] += compute_torque(0);
      94        2010 :             _force_torque_store[6 * i + 4] += compute_torque(1);
      95        2010 :             _force_torque_store[6 * i + 5] += compute_torque(2);
      96             :           }
      97             :       }
      98       47755 : }
      99             : 
     100             : void
     101      246780 : ComputeGrainForceAndTorque::executeJacobian(unsigned int jvar)
     102             : {
     103      246780 :   const auto & op_to_grains = _grain_tracker.getVarToFeatureVector(_current_elem->id());
     104             : 
     105      246780 :   if (jvar == _c_var)
     106      250740 :     for (unsigned int i = 0; i < _grain_num; ++i)
     107      501480 :       for (unsigned int j = 0; j < _op_num; ++j)
     108      334320 :         if (i == op_to_grains[j])
     109             :         {
     110       50064 :           const auto centroid = _grain_tracker.getGrainCentroid(i);
     111      250320 :           for (_qp = 0; _qp < _qrule->n_points(); ++_qp)
     112      200256 :             if (_dFdc[_qp][j](0) != 0.0 || _dFdc[_qp][j](1) != 0.0 || _dFdc[_qp][j](2) != 0.0)
     113             :             {
     114        3624 :               const Real factor = _JxW[_qp] * _coord[_qp] * _phi[_j][_qp];
     115             :               const RealGradient compute_torque_jacobian_c =
     116        3624 :                   factor * (_current_elem->vertex_average() - centroid).cross(_dFdc[_qp][j]);
     117        3624 :               _force_torque_c_jacobian_store[(6 * i + 0) * _total_dofs + _j_global] +=
     118        3624 :                   factor * _dFdc[_qp][j](0);
     119        3624 :               _force_torque_c_jacobian_store[(6 * i + 1) * _total_dofs + _j_global] +=
     120        3624 :                   factor * _dFdc[_qp][j](1);
     121        3624 :               _force_torque_c_jacobian_store[(6 * i + 2) * _total_dofs + _j_global] +=
     122        3624 :                   factor * _dFdc[_qp][j](2);
     123        3624 :               _force_torque_c_jacobian_store[(6 * i + 3) * _total_dofs + _j_global] +=
     124             :                   compute_torque_jacobian_c(0);
     125        3624 :               _force_torque_c_jacobian_store[(6 * i + 4) * _total_dofs + _j_global] +=
     126             :                   compute_torque_jacobian_c(1);
     127        3624 :               _force_torque_c_jacobian_store[(6 * i + 5) * _total_dofs + _j_global] +=
     128             :                   compute_torque_jacobian_c(2);
     129             :             }
     130             :         }
     131             : 
     132      740340 :   for (unsigned int i = 0; i < _op_num; ++i)
     133      493560 :     if (jvar == _vals_var[i])
     134      489600 :       for (unsigned int j = 0; j < _grain_num; ++j)
     135      979200 :         for (unsigned int k = 0; k < _op_num; ++k)
     136      652800 :           if (j == op_to_grains[k])
     137             :           {
     138       95904 :             const auto centroid = _grain_tracker.getGrainCentroid(j);
     139      479520 :             for (_qp = 0; _qp < _qrule->n_points(); ++_qp)
     140      383616 :               if ((*_dFdgradeta[i])[_qp][j] != 0.0)
     141             :               {
     142        7248 :                 const Real factor = _JxW[_qp] * _coord[_qp] * (*_dFdgradeta[i])[_qp][k];
     143             :                 const RealGradient compute_torque_jacobian_eta =
     144        7248 :                     factor * (_current_elem->vertex_average() - centroid).cross(_grad_phi[_j][_qp]);
     145        7248 :                 _force_torque_eta_jacobian_store[i][(6 * j + 0) * _total_dofs + _j_global] +=
     146        7248 :                     factor * _grad_phi[_j][_qp](0);
     147        7248 :                 _force_torque_eta_jacobian_store[i][(6 * j + 1) * _total_dofs + _j_global] +=
     148        7248 :                     factor * _grad_phi[_j][_qp](1);
     149        7248 :                 _force_torque_eta_jacobian_store[i][(6 * j + 2) * _total_dofs + _j_global] +=
     150        7248 :                     factor * _grad_phi[_j][_qp](2);
     151        7248 :                 _force_torque_eta_jacobian_store[i][(6 * j + 3) * _total_dofs + _j_global] +=
     152             :                     compute_torque_jacobian_eta(0);
     153        7248 :                 _force_torque_eta_jacobian_store[i][(6 * j + 4) * _total_dofs + _j_global] +=
     154             :                     compute_torque_jacobian_eta(1);
     155        7248 :                 _force_torque_eta_jacobian_store[i][(6 * j + 5) * _total_dofs + _j_global] +=
     156             :                     compute_torque_jacobian_eta(2);
     157             :               }
     158             :           }
     159      246780 : }
     160             : 
     161             : void
     162         363 : ComputeGrainForceAndTorque::finalize()
     163             : {
     164         363 :   gatherSum(_force_torque_store);
     165        1077 :   for (unsigned int i = 0; i < _grain_num; ++i)
     166             :   {
     167         714 :     _force_values[i](0) = _force_torque_store[6 * i + 0];
     168         714 :     _force_values[i](1) = _force_torque_store[6 * i + 1];
     169         714 :     _force_values[i](2) = _force_torque_store[6 * i + 2];
     170         714 :     _torque_values[i](0) = _force_torque_store[6 * i + 3];
     171         714 :     _torque_values[i](1) = _force_torque_store[6 * i + 4];
     172         714 :     _torque_values[i](2) = _force_torque_store[6 * i + 5];
     173             :   }
     174             : 
     175         363 :   if (_fe_problem.currentlyComputingJacobian())
     176             :   {
     177         158 :     gatherSum(_force_torque_c_jacobian_store);
     178         474 :     for (unsigned int i = 0; i < _op_num; ++i)
     179         316 :       gatherSum(_force_torque_eta_jacobian_store[i]);
     180             :   }
     181         363 : }
     182             : 
     183             : void
     184          53 : ComputeGrainForceAndTorque::threadJoin(const UserObject & y)
     185             : {
     186             :   const auto & pps = static_cast<const ComputeGrainForceAndTorque &>(y);
     187         677 :   for (unsigned int i = 0; i < _ncomp; ++i)
     188         624 :     _force_torque_store[i] += pps._force_torque_store[i];
     189          53 :   if (_fe_problem.currentlyComputingJacobian())
     190             :   {
     191      558203 :     for (unsigned int i = 0; i < _ncomp * _total_dofs; ++i)
     192      558180 :       _force_torque_c_jacobian_store[i] += pps._force_torque_c_jacobian_store[i];
     193          69 :     for (unsigned int i = 0; i < _op_num; ++i)
     194     1116406 :       for (unsigned int j = 0; j < _ncomp * _total_dofs; ++j)
     195     1116360 :         _force_torque_eta_jacobian_store[i][j] += pps._force_torque_eta_jacobian_store[i][j];
     196             :   }
     197          53 : }
     198             : 
     199             : const std::vector<RealGradient> &
     200         298 : ComputeGrainForceAndTorque::getForceValues() const
     201             : {
     202         298 :   return _force_values;
     203             : }
     204             : 
     205             : const std::vector<RealGradient> &
     206         298 : ComputeGrainForceAndTorque::getTorqueValues() const
     207             : {
     208         298 :   return _torque_values;
     209             : }
     210             : 
     211             : const std::vector<Real> &
     212       84312 : ComputeGrainForceAndTorque::getForceCJacobians() const
     213             : {
     214       84312 :   return _force_torque_c_jacobian_store;
     215             : }
     216             : const std::vector<std::vector<Real>> &
     217      168552 : ComputeGrainForceAndTorque::getForceEtaJacobians() const
     218             : {
     219      168552 :   return _force_torque_eta_jacobian_store;
     220             : }