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