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          57 : ComputeExternalGrainForceAndTorque::validParams()
      19             : {
      20          57 :   InputParameters params = ShapeElementUserObject::validParams();
      21          57 :   params.addClassDescription("Userobject for calculating force and torque acting on a grain");
      22         114 :   params.addParam<MaterialPropertyName>("force_density", "force_density", "Force density material");
      23         114 :   params.addParam<UserObjectName>("grain_data", "center of mass of grains");
      24         114 :   params.addCoupledVar("c", "Concentration field");
      25         114 :   params.addCoupledVar("etas", "Array of coupled order parameters");
      26          57 :   return params;
      27           0 : }
      28             : 
      29          30 : ComputeExternalGrainForceAndTorque::ComputeExternalGrainForceAndTorque(
      30          30 :     const InputParameters & parameters)
      31             :   : DerivativeMaterialInterface<ShapeElementUserObject>(parameters),
      32             :     GrainForceAndTorqueInterface(),
      33          30 :     _c_name(coupledName("c", 0)),
      34          30 :     _c_var(coupled("c")),
      35          30 :     _dF_name(getParam<MaterialPropertyName>("force_density")),
      36          30 :     _dF(getMaterialPropertyByName<std::vector<RealGradient>>(_dF_name)),
      37          30 :     _dFdc(getMaterialPropertyByName<std::vector<RealGradient>>(
      38          30 :         derivativePropertyNameFirst(_dF_name, _c_name))),
      39          30 :     _op_num(coupledComponents("etas")),
      40          30 :     _grain_tracker(getUserObject<GrainTrackerInterface>("grain_data")),
      41          30 :     _vals_var(_op_num),
      42          30 :     _vals_name(_op_num),
      43          60 :     _dFdeta(_op_num)
      44             : {
      45         108 :   for (unsigned int i = 0; i < _op_num; ++i)
      46             :   {
      47          78 :     _vals_var[i] = coupled("etas", i);
      48         156 :     _vals_name[i] = coupledName("etas", i);
      49          78 :     _dFdeta[i] = &getMaterialPropertyByName<std::vector<RealGradient>>(
      50         156 :         derivativePropertyNameFirst(_dF_name, _vals_name[i]));
      51             :   }
      52          30 : }
      53             : 
      54             : void
      55         170 : ComputeExternalGrainForceAndTorque::initialize()
      56             : {
      57         170 :   _grain_num = _grain_tracker.getTotalFeatureCount();
      58         170 :   _ncomp = 6 * _grain_num;
      59             : 
      60         170 :   _force_values.resize(_grain_num);
      61         170 :   _torque_values.resize(_grain_num);
      62         170 :   _force_torque_store.assign(_ncomp, 0.0);
      63             : 
      64         170 :   if (_fe_problem.currentlyComputingJacobian())
      65             :   {
      66          40 :     _total_dofs = _subproblem.es().n_dofs();
      67          40 :     _force_torque_c_jacobian_store.assign(_ncomp * _total_dofs, 0.0);
      68          40 :     _force_torque_eta_jacobian_store.resize(_op_num);
      69             : 
      70         120 :     for (unsigned int i = 0; i < _op_num; ++i)
      71          80 :       _force_torque_eta_jacobian_store[i].assign(_ncomp * _total_dofs, 0.0);
      72             :   }
      73         170 : }
      74             : 
      75             : void
      76       29750 : ComputeExternalGrainForceAndTorque::execute()
      77             : {
      78       29750 :   const auto & op_to_grains = _grain_tracker.getVarToFeatureVector(_current_elem->id());
      79             : 
      80       89250 :   for (unsigned int i = 0; i < _grain_num; ++i)
      81      178500 :     for (unsigned int j = 0; j < _op_num; ++j)
      82      119000 :       if (i == op_to_grains[j])
      83             :       {
      84       20468 :         const auto centroid = _grain_tracker.getGrainCentroid(i);
      85      102340 :         for (_qp = 0; _qp < _qrule->n_points(); ++_qp)
      86       81872 :           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       81872 :                 _JxW[_qp] * _coord[_qp] *
      90       81872 :                 (_current_elem->vertex_average() - centroid).cross(_dF[_qp][j]);
      91       81872 :             _force_torque_store[6 * i + 0] += _JxW[_qp] * _coord[_qp] * _dF[_qp][j](0);
      92       81872 :             _force_torque_store[6 * i + 1] += _JxW[_qp] * _coord[_qp] * _dF[_qp][j](1);
      93       81872 :             _force_torque_store[6 * i + 2] += _JxW[_qp] * _coord[_qp] * _dF[_qp][j](2);
      94       81872 :             _force_torque_store[6 * i + 3] += compute_torque(0);
      95       81872 :             _force_torque_store[6 * i + 4] += compute_torque(1);
      96       81872 :             _force_torque_store[6 * i + 5] += compute_torque(2);
      97             :           }
      98             :       }
      99       29750 : }
     100             : 
     101             : void
     102       28000 : ComputeExternalGrainForceAndTorque::executeJacobian(unsigned int jvar)
     103             : {
     104       28000 :   const auto & op_to_grains = _grain_tracker.getVarToFeatureVector(_current_elem->id());
     105             : 
     106       28000 :   if (jvar == _c_var)
     107       84000 :     for (unsigned int i = 0; i < _grain_num; ++i)
     108      168000 :       for (unsigned int j = 0; j < _op_num; ++j)
     109      112000 :         if (i == op_to_grains[j])
     110             :         {
     111       19264 :           const auto centroid = _grain_tracker.getGrainCentroid(i);
     112       96320 :           for (_qp = 0; _qp < _qrule->n_points(); ++_qp)
     113       77056 :             if (_dFdc[_qp][j](0) != 0.0 || _dFdc[_qp][j](1) != 0.0 || _dFdc[_qp][j](2) != 0.0)
     114             :             {
     115       77056 :               const Real factor = _JxW[_qp] * _coord[_qp] * _phi[_j][_qp];
     116             :               const RealGradient compute_torque_jacobian_c =
     117       77056 :                   factor * (_current_elem->vertex_average() - centroid).cross(_dFdc[_qp][j]);
     118       77056 :               _force_torque_c_jacobian_store[(6 * i + 0) * _total_dofs + _j_global] +=
     119       77056 :                   factor * _dFdc[_qp][j](0);
     120       77056 :               _force_torque_c_jacobian_store[(6 * i + 1) * _total_dofs + _j_global] +=
     121       77056 :                   factor * _dFdc[_qp][j](1);
     122       77056 :               _force_torque_c_jacobian_store[(6 * i + 2) * _total_dofs + _j_global] +=
     123       77056 :                   factor * _dFdc[_qp][j](2);
     124       77056 :               _force_torque_c_jacobian_store[(6 * i + 3) * _total_dofs + _j_global] +=
     125             :                   compute_torque_jacobian_c(0);
     126       77056 :               _force_torque_c_jacobian_store[(6 * i + 4) * _total_dofs + _j_global] +=
     127             :                   compute_torque_jacobian_c(1);
     128       77056 :               _force_torque_c_jacobian_store[(6 * i + 5) * _total_dofs + _j_global] +=
     129             :                   compute_torque_jacobian_c(2);
     130             :             }
     131             :         }
     132             : 
     133       84000 :   for (unsigned int i = 0; i < _op_num; ++i)
     134       56000 :     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       28000 : }
     163             : 
     164             : void
     165         153 : ComputeExternalGrainForceAndTorque::finalize()
     166             : {
     167         153 :   gatherSum(_force_torque_store);
     168         459 :   for (unsigned int i = 0; i < _grain_num; ++i)
     169             :   {
     170         306 :     _force_values[i](0) = _force_torque_store[6 * i + 0];
     171         306 :     _force_values[i](1) = _force_torque_store[6 * i + 1];
     172         306 :     _force_values[i](2) = _force_torque_store[6 * i + 2];
     173         306 :     _torque_values[i](0) = _force_torque_store[6 * i + 3];
     174         306 :     _torque_values[i](1) = _force_torque_store[6 * i + 4];
     175         306 :     _torque_values[i](2) = _force_torque_store[6 * i + 5];
     176             :   }
     177             : 
     178         153 :   if (_fe_problem.currentlyComputingJacobian())
     179             :   {
     180          36 :     gatherSum(_force_torque_c_jacobian_store);
     181         108 :     for (unsigned int i = 0; i < _op_num; ++i)
     182          72 :       gatherSum(_force_torque_eta_jacobian_store[i]);
     183             :   }
     184         153 : }
     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          93 : ComputeExternalGrainForceAndTorque::getForceValues() const
     204             : {
     205          93 :   return _force_values;
     206             : }
     207             : 
     208             : const std::vector<RealGradient> &
     209          93 : ComputeExternalGrainForceAndTorque::getTorqueValues() const
     210             : {
     211          93 :   return _torque_values;
     212             : }
     213             : 
     214             : const std::vector<Real> &
     215          66 : ComputeExternalGrainForceAndTorque::getForceCJacobians() const
     216             : {
     217          66 :   return _force_torque_c_jacobian_store;
     218             : }
     219             : const std::vector<std::vector<Real>> &
     220          66 : ComputeExternalGrainForceAndTorque::getForceEtaJacobians() const
     221             : {
     222          66 :   return _force_torque_eta_jacobian_store;
     223             : }