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 "InterfaceMeshCut2DUserObject.h"
      11             : #include "XFEMMovingInterfaceVelocityBase.h"
      12             : 
      13             : registerMooseObject("XFEMApp", InterfaceMeshCut2DUserObject);
      14             : 
      15             : InputParameters
      16          32 : InterfaceMeshCut2DUserObject::validParams()
      17             : {
      18          32 :   InputParameters params = InterfaceMeshCutUserObjectBase::validParams();
      19          32 :   params.addClassDescription("A userobject to cut a 2D mesh using a 1D cutter mesh.");
      20          32 :   return params;
      21           0 : }
      22             : 
      23          16 : InterfaceMeshCut2DUserObject::InterfaceMeshCut2DUserObject(const InputParameters & parameters)
      24          16 :   : InterfaceMeshCutUserObjectBase(parameters)
      25             : {
      26         704 :   for (const auto & elem : _cutter_mesh->element_ptr_range())
      27         336 :     if (elem->type() != EDGE2)
      28           0 :       mooseError(
      29          16 :           "InterfaceMeshCut2DUserObject currently only supports EDGE2 element in the cut mesh.");
      30          16 : }
      31             : 
      32             : void
      33          68 : InterfaceMeshCut2DUserObject::calculateNormals()
      34             : {
      35             :   _element_normals.clear();
      36             : 
      37        1256 :   for (const auto & elem : _cutter_mesh->element_ptr_range())
      38             :   {
      39        1120 :     Point a = elem->node_ref(1);
      40        1120 :     Point b = elem->node_ref(0);
      41             : 
      42        1120 :     Point normal_ab = Point(-(b - a)(1), (b - a)(0), 0);
      43        1120 :     normal_ab /= normal_ab.norm();
      44             : 
      45        1120 :     _element_normals.insert(std::make_pair(elem->id(), normal_ab));
      46          68 :   }
      47          68 : }
      48             : 
      49             : Point
      50        1660 : InterfaceMeshCut2DUserObject::nodeNormal(const unsigned int & node_id)
      51             : {
      52             :   Point normal(0.0);
      53             : 
      54        4860 :   for (const auto & node_neigh_elem_id : _node_to_elem_map[node_id])
      55             :   {
      56        3200 :     const auto & elem = _cutter_mesh->elem_ref(node_neigh_elem_id);
      57             : 
      58        3200 :     Point a = elem.node_ref(1);
      59        3200 :     Point b = elem.node_ref(0);
      60             : 
      61        3200 :     Point normal_ab = Point(-(b - a)(1), (b - a)(0), 0);
      62        3200 :     normal_ab /= normal_ab.norm();
      63             : 
      64             :     normal += normal_ab;
      65             :   }
      66             : 
      67        1660 :   unsigned int num = _node_to_elem_map[node_id].size();
      68             : 
      69        1660 :   if (num == 0)
      70           0 :     mooseError("InterfaceMeshCut2DUserObject, the node is not found in node_to_elem_map in "
      71             :                "calculting its normal.");
      72             : 
      73        1660 :   return normal / num;
      74             : }
      75             : 
      76             : bool
      77        1017 : InterfaceMeshCut2DUserObject::cutElementByGeometry(const Elem * elem,
      78             :                                                    std::vector<Xfem::CutEdge> & cut_edges,
      79             :                                                    std::vector<Xfem::CutNode> & cut_nodes) const
      80             : {
      81             :   mooseAssert(elem->dim() == 2, "Dimension of element to be cut must be 2");
      82             : 
      83             :   bool elem_cut = false;
      84             : 
      85       51642 :   for (const auto & cut_elem : _cutter_mesh->element_ptr_range())
      86             :   {
      87       24804 :     unsigned int n_sides = elem->n_sides();
      88             : 
      89      124020 :     for (unsigned int i = 0; i < n_sides; ++i)
      90             :     {
      91       99216 :       std::unique_ptr<const Elem> curr_side = elem->side_ptr(i);
      92             : 
      93             :       mooseAssert(curr_side->type() == EDGE2, "Element side type must be EDGE2.");
      94             : 
      95             :       const Node * node1 = curr_side->node_ptr(0);
      96             :       const Node * node2 = curr_side->node_ptr(1);
      97       99216 :       Real seg_int_frac = 0.0;
      98             : 
      99       99216 :       const std::pair<Point, Point> elem_endpoints(cut_elem->node_ref(0), cut_elem->node_ref(1));
     100             : 
     101       99216 :       if (Xfem::intersectSegmentWithCutLine(*node1, *node2, elem_endpoints, 1.0, seg_int_frac))
     102             :       {
     103         348 :         if (seg_int_frac > Xfem::tol && seg_int_frac < 1.0 - Xfem::tol)
     104             :         {
     105             :           elem_cut = true;
     106             :           Xfem::CutEdge mycut;
     107         348 :           mycut._id1 = node1->id();
     108         348 :           mycut._id2 = node2->id();
     109         348 :           mycut._distance = seg_int_frac;
     110         348 :           mycut._host_side_id = i;
     111         348 :           cut_edges.push_back(mycut);
     112         348 :         }
     113           0 :         else if (seg_int_frac < Xfem::tol)
     114             :         {
     115             :           elem_cut = true;
     116             :           Xfem::CutNode mycut;
     117           0 :           mycut._id = node1->id();
     118           0 :           mycut._host_id = i;
     119           0 :           cut_nodes.push_back(mycut);
     120             :         }
     121             :       }
     122       99216 :     }
     123        1017 :   }
     124        1017 :   return elem_cut;
     125             : }
     126             : 
     127             : bool
     128           0 : InterfaceMeshCut2DUserObject::cutElementByGeometry(const Elem * /* elem*/,
     129             :                                                    std::vector<Xfem::CutFace> & /*cut_faces*/) const
     130             : {
     131           0 :   mooseError("invalid method for InterfaceMeshCut2DUserObject");
     132             :   return false;
     133             : }
     134             : 
     135             : bool
     136           0 : InterfaceMeshCut2DUserObject::cutFragmentByGeometry(
     137             :     std::vector<std::vector<Point>> & /*frag_edges*/,
     138             :     std::vector<Xfem::CutEdge> & /*cut_edges*/) const
     139             : {
     140           0 :   mooseError("cutFragmentByGeometry not yet implemented for InterfaceMeshCut2DUserObject");
     141             :   return false;
     142             : }
     143             : 
     144             : bool
     145           0 : InterfaceMeshCut2DUserObject::cutFragmentByGeometry(
     146             :     std::vector<std::vector<Point>> & /*frag_faces*/,
     147             :     std::vector<Xfem::CutFace> & /*cut_faces*/) const
     148             : {
     149           0 :   mooseError("invalid method for InterfaceMeshCut2DUserObject");
     150             :   return false;
     151             : }
     152             : 
     153             : Real
     154       42795 : InterfaceMeshCut2DUserObject::calculateSignedDistance(Point p) const
     155             : {
     156             :   Real min_dist = std::numeric_limits<Real>::max();
     157             : 
     158     3206982 :   for (const auto & cut_elem : _cutter_mesh->element_ptr_range())
     159             :   {
     160     1560696 :     Point a = cut_elem->node_ref(0);
     161     1560696 :     Point b = cut_elem->node_ref(1);
     162             : 
     163             :     Point c = p - a;
     164     1560696 :     Point v = (b - a) / (b - a).norm();
     165     1560696 :     Real d = (b - a).norm();
     166             :     Real t = v * c;
     167             : 
     168             :     Real dist;
     169             :     Point nearest_point;
     170             : 
     171     1560696 :     if (t < 0)
     172             :     {
     173      742872 :       dist = (p - a).norm();
     174             :       nearest_point = a;
     175             :     }
     176      817824 :     else if (t > d)
     177             :     {
     178      733968 :       dist = (p - b).norm();
     179             :       nearest_point = b;
     180             :     }
     181             :     else
     182             :     {
     183             :       v *= t;
     184       83856 :       dist = (p - a - v).norm();
     185             :       nearest_point = (a + v);
     186             :     }
     187             : 
     188             :     Point p_nearest_point = nearest_point - p;
     189             : 
     190     1560696 :     Point normal_ab = Point(-(b - a)(1), (b - a)(0), 0);
     191             : 
     192     1560696 :     if (normal_ab * p_nearest_point < 0)
     193     1266774 :       dist = -dist;
     194             : 
     195     1560696 :     if (std::abs(dist) < std::abs(min_dist))
     196             :       min_dist = dist;
     197       42795 :   }
     198             : 
     199       42795 :   return min_dist;
     200             : }