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       26108 : ComputeCrackTipEnrichmentSmallStrain::computeQpProperties()
      73             : {
      74       26108 :   crackTipEnrichementFunctionAtPoint(_q_point[_qp], _B);
      75             :   unsigned int crack_front_point_index =
      76       26108 :       crackTipEnrichementFunctionDerivativeAtPoint(_q_point[_qp], _dBx);
      77             : 
      78      130540 :   for (unsigned int i = 0; i < 4; ++i)
      79      104432 :     rotateFromCrackFrontCoordsToGlobal(_dBx[i], _dBX[i], crack_front_point_index);
      80             : 
      81       26108 :   _sln = _nl->currentSolution();
      82             : 
      83       78324 :   for (unsigned int m = 0; m < _ndisp; ++m)
      84             :   {
      85       52216 :     _enrich_disp[m] = 0.0;
      86             :     _grad_enrich_disp[m].zero();
      87      261080 :     for (unsigned int i = 0; i < _current_elem->n_nodes(); ++i)
      88             :     {
      89      208864 :       const Node * node_i = _current_elem->node_ptr(i);
      90     1044320 :       for (unsigned int j = 0; j < 4; ++j)
      91             :       {
      92      835456 :         dof_id_type dof = node_i->dof_number(_nl->number(), _enrich_variable[j][m]->number(), 0);
      93      835456 :         Real soln = (*_sln)(dof);
      94      835456 :         _enrich_disp[m] += (*_fe_phi)[i][_qp] * (_B[j] - _BI[i][j]) * soln;
      95      835456 :         RealVectorValue grad_B(_dBX[j]);
      96             :         _grad_enrich_disp[m] +=
      97      835456 :             ((*_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       26108 :       _grad_enrich_disp[0], _grad_enrich_disp[1], _grad_enrich_disp[2]);
     104             : 
     105       26108 :   RankTwoTensor enrich_strain = (grad_tensor_enrich + grad_tensor_enrich.transpose()) / 2.0;
     106             : 
     107             :   auto grad_tensor = RankTwoTensor::initializeFromRows(
     108       26108 :       (*_grad_disp[0])[_qp], (*_grad_disp[1])[_qp], (*_grad_disp[2])[_qp]);
     109             : 
     110       26108 :   _total_strain[_qp] = (grad_tensor + grad_tensor.transpose()) / 2.0;
     111             : 
     112       26108 :   _total_strain[_qp] += enrich_strain;
     113             : 
     114       26108 :   _mechanical_strain[_qp] = _total_strain[_qp];
     115             : 
     116             :   // Remove the Eigen strain
     117       26108 :   for (auto es : _eigenstrains)
     118           0 :     _mechanical_strain[_qp] -= (*es)[_qp];
     119       26108 : }
     120             : 
     121             : void
     122        6527 : ComputeCrackTipEnrichmentSmallStrain::computeProperties()
     123             : {
     124        6527 :   FEType fe_type(Utility::string_to_enum<Order>("first"),
     125        6527 :                  Utility::string_to_enum<FEFamily>("lagrange"));
     126        6527 :   const unsigned int dim = _current_elem->dim();
     127        6527 :   std::unique_ptr<FEBase> fe(FEBase::build(dim, fe_type));
     128        6527 :   fe->attach_quadrature_rule(const_cast<QBase *>(_qrule));
     129        6527 :   _fe_phi = &(fe->get_phi());
     130        6527 :   _fe_dphi = &(fe->get_dphi());
     131             : 
     132        6527 :   if (isBoundaryMaterial())
     133           0 :     fe->reinit(_current_elem, _current_side);
     134             :   else
     135        6527 :     fe->reinit(_current_elem);
     136             : 
     137       32635 :   for (unsigned int i = 0; i < _BI.size(); ++i)
     138       26108 :     crackTipEnrichementFunctionAtPoint(*(_current_elem->node_ptr(i)), _BI[i]);
     139             : 
     140        6527 :   ComputeStrainBase::computeProperties();
     141        6527 : }