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          92 : ComputeGrainForceAndTorque::validParams()
      19             : {
      20          92 :   InputParameters params = ShapeElementUserObject::validParams();
      21          92 :   params.addClassDescription("Userobject for calculating force and torque acting on a grain");
      22         184 :   params.addParam<MaterialPropertyName>("force_density", "force_density", "Force density material");
      23         184 :   params.addParam<UserObjectName>("grain_data", "center of mass of grains");
      24         184 :   params.addCoupledVar("c", "Concentration field");
      25         184 :   params.addCoupledVar("etas", "Array of coupled order parameters");
      26          92 :   return params;
      27           0 : }
      28             : 
      29          48 : ComputeGrainForceAndTorque::ComputeGrainForceAndTorque(const InputParameters & parameters)
      30             :   : DerivativeMaterialInterface<ShapeElementUserObject>(parameters),
      31             :     GrainForceAndTorqueInterface(),
      32          48 :     _c_name(coupledName("c", 0)),
      33          48 :     _c_var(coupled("c")),
      34          48 :     _dF_name(getParam<MaterialPropertyName>("force_density")),
      35          48 :     _dF(getMaterialPropertyByName<std::vector<RealGradient>>(_dF_name)),
      36          48 :     _dFdc(getMaterialPropertyByName<std::vector<RealGradient>>(
      37          48 :         derivativePropertyNameFirst(_dF_name, _c_name))),
      38          48 :     _op_num(coupledComponents("etas")),
      39          48 :     _grain_tracker(getUserObject<GrainTrackerInterface>("grain_data")),
      40          48 :     _vals_var(_op_num),
      41          48 :     _vals_name(_op_num),
      42          96 :     _dFdgradeta(_op_num)
      43             : {
      44         144 :   for (unsigned int i = 0; i < _op_num; ++i)
      45             :   {
      46          96 :     _vals_var[i] = coupled("etas", i);
      47         192 :     _vals_name[i] = coupledName("etas", i);
      48          96 :     _dFdgradeta[i] = &getMaterialPropertyByName<std::vector<Real>>(
      49         192 :         derivativePropertyNameFirst(_dF_name, _vals_name[i]));
      50             :   }
      51          48 : }
      52             : 
      53             : void
      54         522 : ComputeGrainForceAndTorque::initialize()
      55             : {
      56         522 :   _grain_num = _grain_tracker.getTotalFeatureCount();
      57         522 :   _ncomp = 6 * _grain_num;
      58             : 
      59         522 :   _force_values.resize(_grain_num);
      60         522 :   _torque_values.resize(_grain_num);
      61         522 :   _force_torque_store.assign(_ncomp, 0.0);
      62             : 
      63         522 :   if (_fe_problem.currentlyComputingJacobian())
      64             :   {
      65         227 :     _total_dofs = _subproblem.es().n_dofs();
      66         227 :     _force_torque_c_jacobian_store.assign(_ncomp * _total_dofs, 0.0);
      67         227 :     _force_torque_eta_jacobian_store.resize(_op_num);
      68             : 
      69         681 :     for (unsigned int i = 0; i < _op_num; ++i)
      70         454 :       _force_torque_eta_jacobian_store[i].assign(_ncomp * _total_dofs, 0.0);
      71             :   }
      72         522 : }
      73             : 
      74             : void
      75       55765 : ComputeGrainForceAndTorque::execute()
      76             : {
      77       55765 :   const auto & op_to_grains = _grain_tracker.getVarToFeatureVector(_current_elem->id());
      78             : 
      79      164895 :   for (unsigned int i = 0; i < _grain_num; ++i)
      80      327390 :     for (unsigned int j = 0; j < _op_num; ++j)
      81      218260 :       if (i == op_to_grains[j])
      82             :       {
      83       32704 :         const auto centroid = _grain_tracker.getGrainCentroid(i);
      84      163520 :         for (_qp = 0; _qp < _qrule->n_points(); ++_qp)
      85      130816 :           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        2345 :                 _JxW[_qp] * _coord[_qp] *
      89        2345 :                 (_current_elem->vertex_average() - centroid).cross(_dF[_qp][j]);
      90        2345 :             _force_torque_store[6 * i + 0] += _JxW[_qp] * _coord[_qp] * _dF[_qp][j](0);
      91        2345 :             _force_torque_store[6 * i + 1] += _JxW[_qp] * _coord[_qp] * _dF[_qp][j](1);
      92        2345 :             _force_torque_store[6 * i + 2] += _JxW[_qp] * _coord[_qp] * _dF[_qp][j](2);
      93        2345 :             _force_torque_store[6 * i + 3] += compute_torque(0);
      94        2345 :             _force_torque_store[6 * i + 4] += compute_torque(1);
      95        2345 :             _force_torque_store[6 * i + 5] += compute_torque(2);
      96             :           }
      97             :       }
      98       55765 : }
      99             : 
     100             : void
     101      287940 : ComputeGrainForceAndTorque::executeJacobian(unsigned int jvar)
     102             : {
     103      287940 :   const auto & op_to_grains = _grain_tracker.getVarToFeatureVector(_current_elem->id());
     104             : 
     105      287940 :   if (jvar == _c_var)
     106      292620 :     for (unsigned int i = 0; i < _grain_num; ++i)
     107      585240 :       for (unsigned int j = 0; j < _op_num; ++j)
     108      390160 :         if (i == op_to_grains[j])
     109             :         {
     110       58440 :           const auto centroid = _grain_tracker.getGrainCentroid(i);
     111      292200 :           for (_qp = 0; _qp < _qrule->n_points(); ++_qp)
     112      233760 :             if (_dFdc[_qp][j](0) != 0.0 || _dFdc[_qp][j](1) != 0.0 || _dFdc[_qp][j](2) != 0.0)
     113             :             {
     114        4228 :               const Real factor = _JxW[_qp] * _coord[_qp] * _phi[_j][_qp];
     115             :               const RealGradient compute_torque_jacobian_c =
     116        4228 :                   factor * (_current_elem->vertex_average() - centroid).cross(_dFdc[_qp][j]);
     117        4228 :               _force_torque_c_jacobian_store[(6 * i + 0) * _total_dofs + _j_global] +=
     118        4228 :                   factor * _dFdc[_qp][j](0);
     119        4228 :               _force_torque_c_jacobian_store[(6 * i + 1) * _total_dofs + _j_global] +=
     120        4228 :                   factor * _dFdc[_qp][j](1);
     121        4228 :               _force_torque_c_jacobian_store[(6 * i + 2) * _total_dofs + _j_global] +=
     122        4228 :                   factor * _dFdc[_qp][j](2);
     123        4228 :               _force_torque_c_jacobian_store[(6 * i + 3) * _total_dofs + _j_global] +=
     124             :                   compute_torque_jacobian_c(0);
     125        4228 :               _force_torque_c_jacobian_store[(6 * i + 4) * _total_dofs + _j_global] +=
     126             :                   compute_torque_jacobian_c(1);
     127        4228 :               _force_torque_c_jacobian_store[(6 * i + 5) * _total_dofs + _j_global] +=
     128             :                   compute_torque_jacobian_c(2);
     129             :             }
     130             :         }
     131             : 
     132      863820 :   for (unsigned int i = 0; i < _op_num; ++i)
     133      575880 :     if (jvar == _vals_var[i])
     134      571200 :       for (unsigned int j = 0; j < _grain_num; ++j)
     135     1142400 :         for (unsigned int k = 0; k < _op_num; ++k)
     136      761600 :           if (j == op_to_grains[k])
     137             :           {
     138      111888 :             const auto centroid = _grain_tracker.getGrainCentroid(j);
     139      559440 :             for (_qp = 0; _qp < _qrule->n_points(); ++_qp)
     140      447552 :               if ((*_dFdgradeta[i])[_qp][j] != 0.0)
     141             :               {
     142        8456 :                 const Real factor = _JxW[_qp] * _coord[_qp] * (*_dFdgradeta[i])[_qp][k];
     143             :                 const RealGradient compute_torque_jacobian_eta =
     144        8456 :                     factor * (_current_elem->vertex_average() - centroid).cross(_grad_phi[_j][_qp]);
     145        8456 :                 _force_torque_eta_jacobian_store[i][(6 * j + 0) * _total_dofs + _j_global] +=
     146        8456 :                     factor * _grad_phi[_j][_qp](0);
     147        8456 :                 _force_torque_eta_jacobian_store[i][(6 * j + 1) * _total_dofs + _j_global] +=
     148        8456 :                     factor * _grad_phi[_j][_qp](1);
     149        8456 :                 _force_torque_eta_jacobian_store[i][(6 * j + 2) * _total_dofs + _j_global] +=
     150        8456 :                     factor * _grad_phi[_j][_qp](2);
     151        8456 :                 _force_torque_eta_jacobian_store[i][(6 * j + 3) * _total_dofs + _j_global] +=
     152             :                     compute_torque_jacobian_eta(0);
     153        8456 :                 _force_torque_eta_jacobian_store[i][(6 * j + 4) * _total_dofs + _j_global] +=
     154             :                     compute_torque_jacobian_eta(1);
     155        8456 :                 _force_torque_eta_jacobian_store[i][(6 * j + 5) * _total_dofs + _j_global] +=
     156             :                     compute_torque_jacobian_eta(2);
     157             :               }
     158             :           }
     159      287940 : }
     160             : 
     161             : void
     162         469 : ComputeGrainForceAndTorque::finalize()
     163             : {
     164         469 :   gatherSum(_force_torque_store);
     165        1391 :   for (unsigned int i = 0; i < _grain_num; ++i)
     166             :   {
     167         922 :     _force_values[i](0) = _force_torque_store[6 * i + 0];
     168         922 :     _force_values[i](1) = _force_torque_store[6 * i + 1];
     169         922 :     _force_values[i](2) = _force_torque_store[6 * i + 2];
     170         922 :     _torque_values[i](0) = _force_torque_store[6 * i + 3];
     171         922 :     _torque_values[i](1) = _force_torque_store[6 * i + 4];
     172         922 :     _torque_values[i](2) = _force_torque_store[6 * i + 5];
     173             :   }
     174             : 
     175         469 :   if (_fe_problem.currentlyComputingJacobian())
     176             :   {
     177         204 :     gatherSum(_force_torque_c_jacobian_store);
     178         612 :     for (unsigned int i = 0; i < _op_num; ++i)
     179         408 :       gatherSum(_force_torque_eta_jacobian_store[i]);
     180             :   }
     181         469 : }
     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         414 : ComputeGrainForceAndTorque::getForceValues() const
     201             : {
     202         414 :   return _force_values;
     203             : }
     204             : 
     205             : const std::vector<RealGradient> &
     206         414 : ComputeGrainForceAndTorque::getTorqueValues() const
     207             : {
     208         414 :   return _torque_values;
     209             : }
     210             : 
     211             : const std::vector<Real> &
     212      110268 : ComputeGrainForceAndTorque::getForceCJacobians() const
     213             : {
     214      110268 :   return _force_torque_c_jacobian_store;
     215             : }
     216             : const std::vector<std::vector<Real>> &
     217      220428 : ComputeGrainForceAndTorque::getForceEtaJacobians() const
     218             : {
     219      220428 :   return _force_torque_eta_jacobian_store;
     220             : }