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 "ComputeCrackTipEnrichmentIncrementalStrain.h"
      11             : #include "MooseMesh.h"
      12             : #include "libmesh/fe_interface.h"
      13             : #include "libmesh/string_to_enum.h"
      14             : 
      15             : #include "libmesh/quadrature.h"
      16             : 
      17             : registerMooseObject("XFEMApp", ComputeCrackTipEnrichmentIncrementalStrain);
      18             : 
      19             : InputParameters
      20          32 : ComputeCrackTipEnrichmentIncrementalStrain::validParams()
      21             : {
      22          32 :   InputParameters params = ComputeIncrementalStrainBase::validParams();
      23          32 :   params.addClassDescription(
      24             :       "Computes the crack tip enrichment strain for an incremental small-strain formulation.");
      25          64 :   params.addRequiredParam<std::vector<NonlinearVariableName>>("enrichment_displacements",
      26             :                                                               "The enrichment displacement");
      27             : 
      28          64 :   params.addRequiredParam<UserObjectName>("crack_front_definition",
      29             :                                           "The CrackFrontDefinition user object name");
      30             : 
      31          32 :   return params;
      32           0 : }
      33             : 
      34          24 : ComputeCrackTipEnrichmentIncrementalStrain::ComputeCrackTipEnrichmentIncrementalStrain(
      35          24 :     const InputParameters & parameters)
      36             :   : ComputeIncrementalStrainBase(parameters),
      37          24 :     EnrichmentFunctionCalculation(&getUserObject<CrackFrontDefinition>("crack_front_definition")),
      38          48 :     _enrich_disp(3),
      39          24 :     _grad_enrich_disp(3),
      40          24 :     _grad_enrich_disp_old(3),
      41          24 :     _enrich_variable(4),
      42          24 :     _phi(_assembly.phi()),
      43          24 :     _grad_phi(_assembly.gradPhi()),
      44          24 :     _grad_enrich_disp_tensor(
      45          24 :         declareProperty<RankTwoTensor>(_base_name + "grad_enrich_disp_tensor")),
      46          24 :     _grad_enrich_disp_tensor_old(
      47          24 :         getMaterialPropertyOld<RankTwoTensor>(_base_name + "grad_enrich_disp_tensor")),
      48          24 :     _B(4),
      49          24 :     _dBX(4),
      50          96 :     _dBx(4)
      51             : {
      52         120 :   for (unsigned int i = 0; i < _enrich_variable.size(); ++i)
      53          96 :     _enrich_variable[i].resize(_ndisp);
      54             : 
      55             :   const std::vector<NonlinearVariableName> & nl_vnames =
      56          24 :       getParam<std::vector<NonlinearVariableName>>("enrichment_displacements");
      57             : 
      58          24 :   if (_ndisp == 2 && nl_vnames.size() != 8)
      59           0 :     mooseError("The number of enrichment displacements should be total 8 for 2D.");
      60          24 :   else if (_ndisp == 3 && nl_vnames.size() != 12)
      61           0 :     mooseError("The number of enrichment displacements should be total 12 for 3D.");
      62             : 
      63          24 :   _nl = &(_fe_problem.getNonlinearSystem(/*nl_sys_num=*/0));
      64             : 
      65          72 :   for (unsigned int j = 0; j < _ndisp; ++j)
      66         240 :     for (unsigned int i = 0; i < 4; ++i)
      67         192 :       _enrich_variable[i][j] = &(_nl->getVariable(0, nl_vnames[j * 4 + i]));
      68             : 
      69          24 :   if (_ndisp == 2)
      70          24 :     _BI.resize(4); // QUAD4
      71           0 :   else if (_ndisp == 3)
      72           0 :     _BI.resize(8); // HEX8
      73             : 
      74         120 :   for (unsigned int i = 0; i < _BI.size(); ++i)
      75          96 :     _BI[i].resize(4);
      76          24 : }
      77             : 
      78             : void
      79        6276 : ComputeCrackTipEnrichmentIncrementalStrain::computeProperties()
      80             : {
      81        6276 :   FEType fe_type(Utility::string_to_enum<Order>("first"),
      82        6276 :                  Utility::string_to_enum<FEFamily>("lagrange"));
      83        6276 :   const unsigned int dim = _current_elem->dim();
      84        6276 :   std::unique_ptr<FEBase> fe(FEBase::build(dim, fe_type));
      85        6276 :   fe->attach_quadrature_rule(const_cast<QBase *>(_qrule));
      86        6276 :   _fe_phi = &(fe->get_phi());
      87        6276 :   _fe_dphi = &(fe->get_dphi());
      88             : 
      89        6276 :   if (isBoundaryMaterial())
      90           0 :     fe->reinit(_current_elem, _current_side);
      91             :   else
      92        6276 :     fe->reinit(_current_elem);
      93             : 
      94        6276 :   _sln = _nl->currentSolution();
      95             : 
      96       31380 :   for (unsigned int i = 0; i < _BI.size(); ++i)
      97       25104 :     crackTipEnrichementFunctionAtPoint(*(_current_elem->node_ptr(i)), _BI[i]);
      98             : 
      99       31380 :   for (_qp = 0; _qp < _qrule->n_points(); ++_qp)
     100             :   {
     101             :     // enrichment function
     102       25104 :     crackTipEnrichementFunctionAtPoint(_q_point[_qp], _B);
     103             :     unsigned int crack_front_point_index =
     104       25104 :         crackTipEnrichementFunctionDerivativeAtPoint(_q_point[_qp], _dBx);
     105             : 
     106      125520 :     for (unsigned int i = 0; i < 4; ++i)
     107      100416 :       rotateFromCrackFrontCoordsToGlobal(_dBx[i], _dBX[i], crack_front_point_index);
     108             : 
     109       75312 :     for (unsigned int m = 0; m < _ndisp; ++m)
     110             :     {
     111       50208 :       _enrich_disp[m] = 0.0;
     112             :       _grad_enrich_disp[m].zero();
     113             :       _grad_enrich_disp_old[m].zero();
     114      251040 :       for (unsigned int i = 0; i < _current_elem->n_nodes(); ++i)
     115             :       {
     116      200832 :         const Node * node_i = _current_elem->node_ptr(i);
     117     1004160 :         for (unsigned int j = 0; j < 4; ++j)
     118             :         {
     119      803328 :           dof_id_type dof = node_i->dof_number(_nl->number(), _enrich_variable[j][m]->number(), 0);
     120      803328 :           Real soln = (*_sln)(dof);
     121      803328 :           _enrich_disp[m] += (*_fe_phi)[i][_qp] * (_B[j] - _BI[i][j]) * soln;
     122      803328 :           RealVectorValue grad_B(_dBX[j]);
     123             : 
     124             :           _grad_enrich_disp[m] +=
     125      803328 :               ((*_fe_dphi)[i][_qp] * (_B[j] - _BI[i][j]) + (*_fe_phi)[i][_qp] * grad_B) * soln;
     126             :         }
     127             :       }
     128             :     }
     129             : 
     130       25104 :     _grad_enrich_disp_tensor[_qp] = RankTwoTensor::initializeFromRows(
     131             :         _grad_enrich_disp[0], _grad_enrich_disp[1], _grad_enrich_disp[2]);
     132             : 
     133             :     auto grad_disp_tensor = RankTwoTensor::initializeFromRows(
     134       25104 :         (*_grad_disp[0])[_qp], (*_grad_disp[1])[_qp], (*_grad_disp[2])[_qp]);
     135             : 
     136             :     auto grad_disp_tensor_old = RankTwoTensor::initializeFromRows(
     137       25104 :         (*_grad_disp_old[0])[_qp], (*_grad_disp_old[1])[_qp], (*_grad_disp_old[2])[_qp]);
     138             : 
     139       25104 :     _deformation_gradient[_qp] = grad_disp_tensor + _grad_enrich_disp_tensor[_qp];
     140       25104 :     _deformation_gradient[_qp].addIa(1.0);
     141             : 
     142       25104 :     _strain_increment[_qp] =
     143       25104 :         0.5 * ((grad_disp_tensor + _grad_enrich_disp_tensor[_qp]) +
     144       25104 :                (grad_disp_tensor + _grad_enrich_disp_tensor[_qp]).transpose()) -
     145       25104 :         0.5 * ((grad_disp_tensor_old + _grad_enrich_disp_tensor_old[_qp]) +
     146       25104 :                (grad_disp_tensor_old + _grad_enrich_disp_tensor_old[_qp]).transpose());
     147             : 
     148             :     // total strain
     149       25104 :     _total_strain[_qp] = _total_strain_old[_qp] + _strain_increment[_qp];
     150             : 
     151             :     // Remove the Eigen strain increment
     152       25104 :     subtractEigenstrainIncrementFromStrain(_strain_increment[_qp]);
     153             : 
     154             :     // strain rate
     155       25104 :     if (_dt > 0)
     156             :     {
     157       24672 :       _strain_rate[_qp] = _strain_increment[_qp] / _dt;
     158       24672 :       _grad_disp_rate[_qp] = (grad_disp_tensor - grad_disp_tensor_old) / _dt;
     159             :     }
     160             :     else
     161             :     {
     162         432 :       _strain_rate[_qp].zero();
     163         432 :       _grad_disp_rate[_qp].zero();
     164             :     }
     165             : 
     166             :     // Update strain in intermediate configuration: rotations are not needed
     167       25104 :     _mechanical_strain[_qp] = _mechanical_strain_old[_qp] + _strain_increment[_qp];
     168             : 
     169             :     // incremental small strain does not include rotation
     170       25104 :     _rotation_increment[_qp].setToIdentity();
     171             :   }
     172        6276 : }