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 "ComputeCrackTipEnrichmentSmallStrain.h"
      11             : #include "MooseMesh.h"
      12             : #include "libmesh/fe_interface.h"
      13             : #include "libmesh/string_to_enum.h"
      14             : 
      15             : #include "libmesh/quadrature_gauss.h"
      16             : 
      17             : registerMooseObject("XFEMApp", ComputeCrackTipEnrichmentSmallStrain);
      18             : 
      19             : InputParameters
      20          16 : ComputeCrackTipEnrichmentSmallStrain::validParams()
      21             : {
      22          16 :   InputParameters params = ComputeStrainBase::validParams();
      23          16 :   params.addClassDescription(
      24             :       "Computes the crack tip enrichment at a point within a small strain formulation.");
      25          32 :   params.addRequiredParam<std::vector<NonlinearVariableName>>("enrichment_displacements",
      26             :                                                               "The enrichment displacement");
      27          32 :   params.addRequiredParam<UserObjectName>("crack_front_definition",
      28             :                                           "The CrackFrontDefinition user object name");
      29          16 :   return params;
      30           0 : }
      31             : 
      32          12 : ComputeCrackTipEnrichmentSmallStrain::ComputeCrackTipEnrichmentSmallStrain(
      33          12 :     const InputParameters & parameters)
      34             :   : ComputeStrainBase(parameters),
      35          12 :     EnrichmentFunctionCalculation(&getUserObject<CrackFrontDefinition>("crack_front_definition")),
      36          24 :     _enrich_disp(3),
      37          12 :     _grad_enrich_disp(3),
      38          12 :     _enrich_variable(4),
      39          12 :     _phi(_assembly.phi()),
      40          12 :     _grad_phi(_assembly.gradPhi()),
      41          12 :     _B(4),
      42          12 :     _dBX(4),
      43          48 :     _dBx(4)
      44             : {
      45          60 :   for (unsigned int i = 0; i < _enrich_variable.size(); ++i)
      46          48 :     _enrich_variable[i].resize(_ndisp);
      47             : 
      48             :   const std::vector<NonlinearVariableName> & nl_vnames =
      49          12 :       getParam<std::vector<NonlinearVariableName>>("enrichment_displacements");
      50             : 
      51          12 :   if (_ndisp == 2 && nl_vnames.size() != 8)
      52           0 :     mooseError("The number of enrichment displacements should be total 8 for 2D.");
      53          12 :   else if (_ndisp == 3 && nl_vnames.size() != 12)
      54           0 :     mooseError("The number of enrichment displacements should be total 12 for 3D.");
      55             : 
      56          12 :   _nl = &(_fe_problem.getNonlinearSystem(/*nl_sys_num=*/0));
      57             : 
      58          36 :   for (unsigned int j = 0; j < _ndisp; ++j)
      59         120 :     for (unsigned int i = 0; i < 4; ++i)
      60          96 :       _enrich_variable[i][j] = &(_nl->getVariable(0, nl_vnames[j * 4 + i]));
      61             : 
      62          12 :   if (_ndisp == 2)
      63          12 :     _BI.resize(4); // QUAD4
      64           0 :   else if (_ndisp == 3)
      65           0 :     _BI.resize(8); // HEX8
      66             : 
      67          60 :   for (unsigned int i = 0; i < _BI.size(); ++i)
      68          48 :     _BI[i].resize(4);
      69          12 : }
      70             : 
      71             : void
      72       28588 : ComputeCrackTipEnrichmentSmallStrain::computeQpProperties()
      73             : {
      74       28588 :   crackTipEnrichementFunctionAtPoint(_q_point[_qp], _B);
      75             :   unsigned int crack_front_point_index =
      76       28588 :       crackTipEnrichementFunctionDerivativeAtPoint(_q_point[_qp], _dBx);
      77             : 
      78      142940 :   for (unsigned int i = 0; i < 4; ++i)
      79      114352 :     rotateFromCrackFrontCoordsToGlobal(_dBx[i], _dBX[i], crack_front_point_index);
      80             : 
      81       28588 :   _sln = _nl->currentSolution();
      82             : 
      83       85764 :   for (unsigned int m = 0; m < _ndisp; ++m)
      84             :   {
      85       57176 :     _enrich_disp[m] = 0.0;
      86             :     _grad_enrich_disp[m].zero();
      87      285880 :     for (unsigned int i = 0; i < _current_elem->n_nodes(); ++i)
      88             :     {
      89      228704 :       const Node * node_i = _current_elem->node_ptr(i);
      90     1143520 :       for (unsigned int j = 0; j < 4; ++j)
      91             :       {
      92      914816 :         dof_id_type dof = node_i->dof_number(_nl->number(), _enrich_variable[j][m]->number(), 0);
      93      914816 :         Real soln = (*_sln)(dof);
      94      914816 :         _enrich_disp[m] += (*_fe_phi)[i][_qp] * (_B[j] - _BI[i][j]) * soln;
      95      914816 :         RealVectorValue grad_B(_dBX[j]);
      96             :         _grad_enrich_disp[m] +=
      97      914816 :             ((*_fe_dphi)[i][_qp] * (_B[j] - _BI[i][j]) + (*_fe_phi)[i][_qp] * grad_B) * soln;
      98             :       }
      99             :     }
     100             :   }
     101             : 
     102             :   auto grad_tensor_enrich = RankTwoTensor::initializeFromRows(
     103       28588 :       _grad_enrich_disp[0], _grad_enrich_disp[1], _grad_enrich_disp[2]);
     104             : 
     105       28588 :   RankTwoTensor enrich_strain = (grad_tensor_enrich + grad_tensor_enrich.transpose()) / 2.0;
     106             : 
     107             :   auto grad_tensor = RankTwoTensor::initializeFromRows(
     108       28588 :       (*_grad_disp[0])[_qp], (*_grad_disp[1])[_qp], (*_grad_disp[2])[_qp]);
     109             : 
     110       28588 :   _total_strain[_qp] = (grad_tensor + grad_tensor.transpose()) / 2.0;
     111             : 
     112       28588 :   _total_strain[_qp] += enrich_strain;
     113             : 
     114       28588 :   _mechanical_strain[_qp] = _total_strain[_qp];
     115             : 
     116             :   // Remove the Eigen strain
     117       28588 :   for (auto es : _eigenstrains)
     118           0 :     _mechanical_strain[_qp] -= (*es)[_qp];
     119       28588 : }
     120             : 
     121             : void
     122        7767 : ComputeCrackTipEnrichmentSmallStrain::computeProperties()
     123             : {
     124        7767 :   FEType fe_type(Utility::string_to_enum<Order>("first"),
     125        7767 :                  Utility::string_to_enum<FEFamily>("lagrange"));
     126        7767 :   const unsigned int dim = _current_elem->dim();
     127        7767 :   std::unique_ptr<FEBase> fe(FEBase::build(dim, fe_type));
     128        7767 :   fe->attach_quadrature_rule(const_cast<QBase *>(_qrule));
     129        7767 :   _fe_phi = &(fe->get_phi());
     130        7767 :   _fe_dphi = &(fe->get_dphi());
     131             : 
     132        7767 :   if (isBoundaryMaterial())
     133        1240 :     fe->reinit(_current_elem, _current_side);
     134             :   else
     135        6527 :     fe->reinit(_current_elem);
     136             : 
     137       38835 :   for (unsigned int i = 0; i < _BI.size(); ++i)
     138       31068 :     crackTipEnrichementFunctionAtPoint(*(_current_elem->node_ptr(i)), _BI[i]);
     139             : 
     140        7767 :   ComputeStrainBase::computeProperties();
     141        7767 : }