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 "GeometricCut3DUserObject.h"
      11             : 
      12             : // MOOSE includes
      13             : #include "MooseError.h"
      14             : 
      15             : #include "libmesh/string_to_enum.h"
      16             : 
      17             : // XFEM includes
      18             : #include "XFEMFuncs.h"
      19             : 
      20             : InputParameters
      21         112 : GeometricCut3DUserObject::validParams()
      22             : {
      23         112 :   InputParameters params = GeometricCutUserObject::validParams();
      24         112 :   params.addClassDescription("Base class for 3D XFEM Geometric Cut UserObjects");
      25         112 :   return params;
      26           0 : }
      27             : 
      28          56 : GeometricCut3DUserObject::GeometricCut3DUserObject(const InputParameters & parameters)
      29          56 :   : GeometricCutUserObject(parameters), _center(), _normal()
      30             : {
      31          56 : }
      32             : 
      33             : bool
      34           0 : GeometricCut3DUserObject::cutElementByGeometry(const Elem * /*elem*/,
      35             :                                                std::vector<Xfem::CutEdge> & /*cut_edges*/,
      36             :                                                std::vector<Xfem::CutNode> & /*cut_nodes*/) const
      37             : {
      38           0 :   mooseError("Invalid method: must use vector of element faces for 3D mesh cutting");
      39             :   return false;
      40             : }
      41             : 
      42             : bool
      43       16956 : GeometricCut3DUserObject::cutElementByGeometry(const Elem * elem,
      44             :                                                std::vector<Xfem::CutFace> & cut_faces) const
      45             : // TODO: Time evolving cuts not yet supported in 3D (hence the lack of use of the time variable)
      46             : {
      47             :   bool cut_elem = false;
      48             : 
      49       97316 :   for (unsigned int i = 0; i < elem->n_sides(); ++i)
      50             :   {
      51             :     // This returns the lowest-order type of side.
      52       80360 :     std::unique_ptr<const Elem> curr_side = elem->side_ptr(i);
      53       80360 :     if (curr_side->dim() != 2)
      54           0 :       mooseError("In cutElementByGeometry dimension of side must be 2, but it is ",
      55           0 :                  curr_side->dim());
      56       80360 :     unsigned int n_edges = curr_side->n_sides();
      57             : 
      58             :     std::vector<unsigned int> cut_edges;
      59             :     std::vector<Real> cut_pos;
      60             : 
      61      359048 :     for (unsigned int j = 0; j < n_edges; j++)
      62             :     {
      63             :       // This returns the lowest-order type of side.
      64      278688 :       std::unique_ptr<const Elem> curr_edge = curr_side->side_ptr(j);
      65      278688 :       if (curr_edge->type() != EDGE2)
      66           0 :         mooseError("In cutElementByGeometry face edge must be EDGE2, but type is: ",
      67           0 :                    libMesh::Utility::enum_to_string(curr_edge->type()),
      68             :                    " base element type is: ",
      69           0 :                    libMesh::Utility::enum_to_string(elem->type()));
      70             :       const Node * node1 = curr_edge->node_ptr(0);
      71             :       const Node * node2 = curr_edge->node_ptr(1);
      72             : 
      73             :       Point intersection;
      74      278688 :       if (intersectWithEdge(*node1, *node2, intersection))
      75             :       {
      76       20804 :         cut_edges.push_back(j);
      77       20804 :         cut_pos.push_back(getRelativePosition(*node1, *node2, intersection));
      78             :       }
      79      278688 :     }
      80             : 
      81       80360 :     if (cut_edges.size() == 2)
      82             :     {
      83             :       cut_elem = true;
      84             :       Xfem::CutFace mycut;
      85       10130 :       mycut._face_id = i;
      86       10130 :       mycut._face_edge.push_back(cut_edges[0]);
      87       10130 :       mycut._face_edge.push_back(cut_edges[1]);
      88       10130 :       mycut._position.push_back(cut_pos[0]);
      89       10130 :       mycut._position.push_back(cut_pos[1]);
      90       10130 :       cut_faces.push_back(mycut);
      91             :     }
      92       80360 :   }
      93             : 
      94       16956 :   return cut_elem;
      95             : }
      96             : 
      97             : bool
      98           0 : GeometricCut3DUserObject::cutFragmentByGeometry(std::vector<std::vector<Point>> & /*frag_edges*/,
      99             :                                                 std::vector<Xfem::CutEdge> & /*cut_edges*/) const
     100             : {
     101           0 :   mooseError("Invalid method: must use vector of element faces for 3D mesh cutting");
     102             :   return false;
     103             : }
     104             : 
     105             : bool
     106           0 : GeometricCut3DUserObject::cutFragmentByGeometry(std::vector<std::vector<Point>> & /*frag_faces*/,
     107             :                                                 std::vector<Xfem::CutFace> & /*cut_faces*/) const
     108             : {
     109             :   // TODO: Need this for branching in 3D
     110           0 :   mooseError("cutFragmentByGeometry not yet implemented for 3D mesh cutting");
     111             :   return false;
     112             : }
     113             : 
     114             : bool
     115      278688 : GeometricCut3DUserObject::intersectWithEdge(const Point & p1, const Point & p2, Point & pint) const
     116             : {
     117             :   bool has_intersection = false;
     118      278688 :   double plane_point[3] = {_center(0), _center(1), _center(2)};
     119      278688 :   double plane_normal[3] = {_normal(0), _normal(1), _normal(2)};
     120      278688 :   double edge_point1[3] = {p1(0), p1(1), p1(2)};
     121      278688 :   double edge_point2[3] = {p2(0), p2(1), p2(2)};
     122      278688 :   double cut_point[3] = {0.0, 0.0, 0.0};
     123             : 
     124      278688 :   if (Xfem::plane_normal_line_exp_int_3d(
     125             :           plane_point, plane_normal, edge_point1, edge_point2, cut_point) == 1)
     126             :   {
     127      141032 :     Point temp_p(cut_point[0], cut_point[1], cut_point[2]);
     128      141032 :     if (isInsideCutPlane(temp_p) && isInsideEdge(p1, p2, temp_p))
     129             :     {
     130       20804 :       pint = temp_p;
     131             :       has_intersection = true;
     132             :     }
     133             :   }
     134      278688 :   return has_intersection;
     135             : }
     136             : 
     137             : bool
     138       69876 : GeometricCut3DUserObject::isInsideEdge(const Point & p1, const Point & p2, const Point & p) const
     139             : {
     140             :   Real dotp1 = (p1 - p) * (p2 - p1);
     141             :   Real dotp2 = (p2 - p) * (p2 - p1);
     142       69876 :   return (dotp1 * dotp2 <= 0.0);
     143             : }
     144             : 
     145             : Real
     146       20804 : GeometricCut3DUserObject::getRelativePosition(const Point & p1,
     147             :                                               const Point & p2,
     148             :                                               const Point & p) const
     149             : {
     150             :   // get the relative position of p from p1
     151       20804 :   Real full_len = (p2 - p1).norm();
     152       20804 :   Real len_p1_p = (p - p1).norm();
     153       20804 :   return len_p1_p / full_len;
     154             : }