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 "InterfaceMeshCutUserObjectBase.h"
      11             : #include "XFEMMovingInterfaceVelocityBase.h"
      12             : 
      13             : InputParameters
      14          48 : InterfaceMeshCutUserObjectBase::validParams()
      15             : {
      16          48 :   InputParameters params = GeometricCutUserObject::validParams();
      17          96 :   params.addRequiredParam<MeshFileName>("mesh_file", "Mesh file for the XFEM geometric cut.");
      18          96 :   params.addParam<UserObjectName>("interface_velocity_uo",
      19             :                                   "The name of userobject that computes the velocity.");
      20          96 :   params.addParam<FunctionName>("interface_velocity_function",
      21             :                                 "The name of function that provides interface velocity.");
      22          96 :   params.addParam<bool>("output_exodus", true, "Whether to output exodus file of cutter mesh.");
      23          96 :   params.addParam<CutSubdomainID>(
      24          96 :       "negative_id", 0, "The CutSubdomainID corresponding to a non-positive signed distance");
      25          96 :   params.addParam<CutSubdomainID>(
      26          96 :       "positive_id", 1, "The CutSubdomainID corresponding to a positive signed distance");
      27          48 :   return params;
      28           0 : }
      29             : 
      30          24 : InterfaceMeshCutUserObjectBase::InterfaceMeshCutUserObjectBase(const InputParameters & parameters)
      31             :   : GeometricCutUserObject(parameters),
      32          24 :     _interface_velocity(nullptr),
      33          24 :     _func(nullptr),
      34          24 :     _mesh(_subproblem.mesh()),
      35             :     _exodus_io(nullptr),
      36          24 :     _explicit_system(nullptr),
      37             :     _equation_systems(nullptr),
      38          48 :     _output_exodus(getParam<bool>("output_exodus")),
      39          48 :     _negative_id(getParam<CutSubdomainID>("negative_id")),
      40          72 :     _positive_id(getParam<CutSubdomainID>("positive_id"))
      41             : {
      42          24 :   MeshFileName xfem_cutter_mesh_file = getParam<MeshFileName>("mesh_file");
      43          48 :   _cutter_mesh = std::make_shared<ReplicatedMesh>(_communicator);
      44          24 :   _cutter_mesh->read(xfem_cutter_mesh_file);
      45          48 :   _exodus_io = std::make_unique<ExodusII_IO>(*_cutter_mesh);
      46             : 
      47          60 :   if (!parameters.isParamSetByUser("interface_velocity_uo") &&
      48          48 :       !parameters.isParamSetByUser("interface_velocity_function"))
      49           0 :     paramError("Either `interface_velocity_uo` or `interface_velocity_function` must be provided.");
      50          48 :   else if (parameters.isParamSetByUser("interface_velocity_uo") &&
      51          36 :            parameters.isParamSetByUser("interface_velocity_function"))
      52           0 :     paramError(
      53             :         "'interface_velocity_uo` and `interface_velocity_function` cannot be both provided.");
      54             : 
      55          48 :   if (parameters.isParamSetByUser("interface_velocity_function"))
      56          12 :     _func = &getFunction("interface_velocity_function");
      57          24 : }
      58             : 
      59             : void
      60          24 : InterfaceMeshCutUserObjectBase::initialSetup()
      61             : {
      62          24 :   if (_func == nullptr)
      63             :   {
      64             :     const UserObject * uo =
      65          24 :         &(_fe_problem.getUserObjectBase(getParam<UserObjectName>("interface_velocity_uo")));
      66             : 
      67          12 :     if (dynamic_cast<const XFEMMovingInterfaceVelocityBase *>(uo) == nullptr)
      68           0 :       mooseError("UserObject casting to XFEMMovingInterfaceVelocityBase in "
      69             :                  "MovingLineSegmentCutSetUserObject");
      70             : 
      71          12 :     _interface_velocity = dynamic_cast<const XFEMMovingInterfaceVelocityBase *>(uo);
      72          12 :     const_cast<XFEMMovingInterfaceVelocityBase *>(_interface_velocity)->initialize();
      73             :   }
      74             : 
      75        1804 :   for (const auto & node : _cutter_mesh->node_ptr_range())
      76        1780 :     _initial_nodes_location[node->id()] = *node;
      77             : 
      78        5376 :   for (const auto & elem : _cutter_mesh->element_ptr_range())
      79       10320 :     for (unsigned int n = 0; n < elem->n_nodes(); n++)
      80        7680 :       _node_to_elem_map[elem->node_id(n)].push_back(elem->id());
      81             : 
      82          24 :   _cutter_mesh->prepare_for_use();
      83          24 :   _cutter_mesh->set_mesh_dimension(_mesh.dimension() - 1);
      84             : 
      85          24 :   if (_output_exodus)
      86             :   {
      87          48 :     _equation_systems = std::make_unique<EquationSystems>(*_cutter_mesh);
      88          24 :     _explicit_system = &(_equation_systems->add_system<ExplicitSystem>("InterfaceMeshSystem"));
      89             : 
      90          24 :     _explicit_system->add_variable("disp_x");
      91          24 :     _explicit_system->add_variable("disp_y");
      92             : 
      93          24 :     if (_mesh.dimension() == 3)
      94           8 :       _explicit_system->add_variable("disp_z");
      95             : 
      96          24 :     _equation_systems->init();
      97          72 :     _exodus_io->write_equation_systems(_app.getOutputFileBase() + "_" + name() + ".e",
      98             :                                        *_equation_systems);
      99             : 
     100          24 :     _var_num_disp_x = _explicit_system->variable_number("disp_x");
     101          24 :     _var_num_disp_y = _explicit_system->variable_number("disp_y");
     102          24 :     if (_mesh.dimension() == 3)
     103           8 :       _var_num_disp_z = _explicit_system->variable_number("disp_z");
     104             :   }
     105          24 : }
     106             : 
     107             : void
     108         100 : InterfaceMeshCutUserObjectBase::initialize()
     109             : {
     110         100 :   std::vector<Point> new_position(_cutter_mesh->n_nodes());
     111             : 
     112         200 :   _pl = _mesh.getPointLocator();
     113         100 :   _pl->enable_out_of_mesh_mode();
     114             : 
     115             :   std::map<unsigned int, Real> node_velocity;
     116             :   Real sum = 0.0;
     117             :   unsigned count = 0;
     118             :   // Loop all the nodes to calculate the velocity
     119       13376 :   for (const auto & node : _cutter_mesh->node_ptr_range())
     120             :   {
     121        6588 :     if ((*_pl)(*node) != nullptr)
     122             :     {
     123             :       Real velocity;
     124        4476 :       if (_func == nullptr)
     125             :         velocity =
     126        3312 :             _interface_velocity->computeMovingInterfaceVelocity(node->id(), nodeNormal(node->id()));
     127             :       else
     128        1164 :         velocity = _func->value(_t, *node);
     129             : 
     130             :       // only updates when t_step >0
     131        4476 :       if (_t_step <= 0)
     132             :         velocity = 0.0;
     133             : 
     134        4476 :       node_velocity[node->id()] = velocity;
     135        4476 :       sum += velocity;
     136        4476 :       count++;
     137             :     }
     138         100 :   }
     139             : 
     140         100 :   if (count == 0)
     141           0 :     mooseError("No node of the cutter mesh is found inside the computational domain.");
     142             : 
     143         100 :   Real average_velocity = sum / count;
     144             : 
     145       13376 :   for (const auto & node : _cutter_mesh->node_ptr_range())
     146             :   {
     147        6588 :     if ((*_pl)(*node) == nullptr)
     148        2112 :       node_velocity[node->id()] = average_velocity;
     149         100 :   }
     150             : 
     151       13376 :   for (const auto & node : _cutter_mesh->node_ptr_range())
     152             :   {
     153        6588 :     Point p = *node;
     154        6588 :     p += _dt * nodeNormal(node->id()) * node_velocity[node->id()];
     155        6588 :     new_position[node->id()] = p;
     156         100 :   }
     157       13376 :   for (const auto & node : _cutter_mesh->node_ptr_range())
     158        6688 :     _cutter_mesh->node_ref(node->id()) = new_position[node->id()];
     159             : 
     160         100 :   if (_output_exodus)
     161             :   {
     162             :     std::vector<dof_id_type> di;
     163       13376 :     for (const auto & node : _cutter_mesh->node_ptr_range())
     164             :     {
     165        6588 :       _explicit_system->get_dof_map().dof_indices(node, di, _var_num_disp_x);
     166        6588 :       _explicit_system->solution->set(
     167        6588 :           di[0], new_position[node->id()](0) - _initial_nodes_location[node->id()](0));
     168             : 
     169        6588 :       _explicit_system->get_dof_map().dof_indices(node, di, _var_num_disp_y);
     170        6588 :       _explicit_system->solution->set(
     171        6588 :           di[0], new_position[node->id()](1) - _initial_nodes_location[node->id()](1));
     172             : 
     173        6588 :       if (_mesh.dimension() == 3)
     174             :       {
     175        5408 :         _explicit_system->get_dof_map().dof_indices(node, di, _var_num_disp_z);
     176        5408 :         _explicit_system->solution->set(
     177        5408 :             di[0], new_position[node->id()](2) - _initial_nodes_location[node->id()](2));
     178             :       }
     179         100 :     }
     180             : 
     181         100 :     _explicit_system->solution->close();
     182             : 
     183         100 :     _exodus_io->append(true);
     184         100 :     _exodus_io->write_timestep(
     185         200 :         _app.getOutputFileBase() + "_" + name() + ".e", *_equation_systems, _t_step + 1, _t);
     186         100 :   }
     187             : 
     188         100 :   calculateNormals();
     189         100 : }
     190             : 
     191             : const std::vector<Point>
     192           0 : InterfaceMeshCutUserObjectBase::getCrackFrontPoints(unsigned int /*num_crack_front_points*/) const
     193             : {
     194           0 :   mooseError("getCrackFrontPoints() is not implemented for this object.");
     195             : }
     196             : 
     197             : const std::vector<RealVectorValue>
     198           0 : InterfaceMeshCutUserObjectBase::getCrackPlaneNormals(unsigned int /*num_crack_front_points*/) const
     199             : {
     200           0 :   mooseError("getCrackPlaneNormals() is not implemented for this object.");
     201             : }
     202             : 
     203             : CutSubdomainID
     204        1779 : InterfaceMeshCutUserObjectBase::getCutSubdomainID(const Node * node) const
     205             : {
     206        1779 :   return calculateSignedDistance(*node) > 0.0 ? _positive_id : _negative_id;
     207             : }