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 "SingularShapeTensorEliminatorUserObjectPD.h"
      11             : #include "AuxiliarySystem.h"
      12             : 
      13             : registerMooseObject("PeridynamicsApp", SingularShapeTensorEliminatorUserObjectPD);
      14             : 
      15             : InputParameters
      16          10 : SingularShapeTensorEliminatorUserObjectPD::validParams()
      17             : {
      18          10 :   InputParameters params = GeneralUserObjectBasePD::validParams();
      19          10 :   params.addClassDescription("UserObject to eliminate the existance of singular shape tensor");
      20             : 
      21          10 :   return params;
      22           0 : }
      23             : 
      24           5 : SingularShapeTensorEliminatorUserObjectPD::SingularShapeTensorEliminatorUserObjectPD(
      25           5 :     const InputParameters & parameters)
      26           5 :   : GeneralUserObjectBasePD(parameters)
      27             : {
      28           5 : }
      29             : 
      30             : void
      31           6 : SingularShapeTensorEliminatorUserObjectPD::initialize()
      32             : {
      33           6 :   _aux.solution().close();
      34           6 : }
      35             : 
      36             : void
      37           6 : SingularShapeTensorEliminatorUserObjectPD::execute()
      38             : {
      39             :   bool converged = false;
      40             : 
      41             :   // Loop through the active local elements to check the shape tensor singularity
      42           6 :   auto first_elem = _mesh.getMesh().active_local_elements_begin();
      43           6 :   auto last_elem = _mesh.getMesh().active_local_elements_end();
      44             : 
      45           6 :   unsigned int singularity_count, loop_count = 0;
      46             :   bool should_print = true; // for printing purpose only
      47             : 
      48          18 :   while (!converged)
      49             :   {
      50          12 :     singularity_count = 0;
      51       10972 :     for (auto elem = first_elem; elem != last_elem; ++elem)
      52             :     {
      53             :       // shape tensor singularity check only applies to intact Edge2 elems
      54        5480 :       if ((*elem)->type() == 0 && _bond_status_var->getElementalValue(*elem) > 0.5)
      55             :       {
      56        2330 :         if (checkShapeTensorSingularity(*elem))
      57             :         {
      58          50 :           dof_id_type dof = (*elem)->dof_number(_aux.number(), _bond_status_var->number(), 0);
      59          50 :           _aux.solution().set(dof, 0); // treat bonds with singular shape tensor as broken
      60             : 
      61          50 :           singularity_count++;
      62             :         }
      63             :       }
      64             :     }
      65             : 
      66             :     gatherSum(singularity_count); // gather the value across processors
      67             : 
      68          12 :     if (singularity_count == 0)
      69             :       converged = true;
      70             :     else // sync aux across processors
      71             :     {
      72           6 :       _aux.solution().close();
      73           6 :       _aux.update();
      74             :     }
      75             : 
      76          12 :     if (singularity_count > 0 && should_print) // print once
      77             :     {
      78           7 :       _console << COLOR_MAGENTA << " Singular shape tensor detected! Elimination in process ... "
      79           7 :                << COLOR_DEFAULT << std::endl;
      80             :       should_print = false;
      81             :     }
      82             : 
      83          12 :     loop_count++;
      84          12 :     if ((loop_count == 1 && singularity_count != 0) || loop_count > 1)
      85          14 :       _console << COLOR_MAGENTA << "     Loop: " << loop_count
      86          14 :                << ", Singularities: " << singularity_count << COLOR_DEFAULT << std::endl;
      87             :   }
      88             : 
      89           6 :   if (loop_count > 1)
      90           8 :     _console << COLOR_MAGENTA << " Elimination done!" << COLOR_DEFAULT << std::endl;
      91           6 : }
      92             : 
      93             : void
      94           6 : SingularShapeTensorEliminatorUserObjectPD::finalize()
      95             : {
      96           6 :   _aux.solution().close();
      97           6 : }
      98             : 
      99             : bool
     100        2330 : SingularShapeTensorEliminatorUserObjectPD::checkShapeTensorSingularity(const Elem * elem)
     101             : {
     102             :   bool singular = false;
     103             : 
     104        2330 :   RankTwoTensor shape_tensor;
     105             :   Real vol_nb, weight_nb, horizon_nd;
     106             :   RealGradient origin_vec_nb;
     107        6990 :   for (unsigned int nd = 0; nd < _nnodes; ++nd)
     108             :   {
     109        4660 :     std::vector<dof_id_type> neighbors = _pdmesh.getNeighbors(elem->node_id(nd));
     110             :     std::vector<dof_id_type> bonds =
     111        4660 :         _pdmesh.getBonds(elem->node_id(nd)); // potentially includes ghosted elems
     112        4660 :     horizon_nd = _pdmesh.getHorizon(elem->node_id(nd));
     113             : 
     114             :     shape_tensor.zero();
     115        4660 :     if (_dim == 2)
     116        4660 :       shape_tensor(2, 2) = 1.0;
     117             : 
     118       93050 :     for (unsigned int nb = 0; nb < neighbors.size(); ++nb)
     119       88390 :       if (_bond_status_var->getElementalValue(_pdmesh.elemPtr(bonds[nb])) > 0.5)
     120             :       {
     121       60560 :         vol_nb = _pdmesh.getNodeVolume(neighbors[nb]);
     122       60560 :         origin_vec_nb =
     123       60560 :             _pdmesh.getNodeCoord(neighbors[nb]) - _pdmesh.getNodeCoord(elem->node_id(nd));
     124       60560 :         weight_nb = horizon_nd / origin_vec_nb.norm();
     125             : 
     126      181680 :         for (unsigned int k = 0; k < _dim; ++k)
     127      363360 :           for (unsigned int l = 0; l < _dim; ++l)
     128      242240 :             shape_tensor(k, l) += weight_nb * origin_vec_nb(k) * origin_vec_nb(l) * vol_nb;
     129             :       }
     130             : 
     131        4660 :     if (shape_tensor.det() == 0.)
     132             :       singular = true;
     133        4660 :   }
     134             : 
     135        2330 :   return singular;
     136             : }