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 "InterfaceMeshCut3DUserObject.h"
      11             : #include "XFEMMovingInterfaceVelocityBase.h"
      12             : #include "libmesh/enum_point_locator_type.h"
      13             : 
      14             : registerMooseObject("XFEMApp", InterfaceMeshCut3DUserObject);
      15             : 
      16             : InputParameters
      17          16 : InterfaceMeshCut3DUserObject::validParams()
      18             : {
      19          16 :   InputParameters params = InterfaceMeshCutUserObjectBase::validParams();
      20          16 :   params.addClassDescription("A userobject to cut a 3D mesh using a 2D cutter mesh.");
      21          16 :   return params;
      22           0 : }
      23             : 
      24           8 : InterfaceMeshCut3DUserObject::InterfaceMeshCut3DUserObject(const InputParameters & parameters)
      25           8 :   : InterfaceMeshCutUserObjectBase(parameters)
      26             : {
      27        4672 :   for (const auto & elem : _cutter_mesh->element_ptr_range())
      28        2328 :     if (elem->type() != TRI3)
      29           0 :       mooseError("InterfaceMeshCut3DUserObject currently only supports TRI3 elements in the "
      30           8 :                  "cutting mesh.");
      31           8 : }
      32             : 
      33             : void
      34          32 : InterfaceMeshCut3DUserObject::calculateNormals()
      35             : {
      36             :   _pseudo_normal.clear();
      37             : 
      38       18272 :   for (const auto & elem : _cutter_mesh->element_ptr_range())
      39             :   {
      40        9104 :     std::vector<Point> vertices{elem->node_ref(0), elem->node_ref(1), elem->node_ref(2)};
      41             :     std::array<Point, 7> normal;
      42        9104 :     Plane elem_plane(vertices[0], vertices[1], vertices[2]);
      43        9104 :     normal[0] = 2.0 * libMesh::pi * elem_plane.unit_normal(vertices[0]);
      44             : 
      45       36416 :     for (unsigned int i = 0; i < elem->n_nodes(); i++)
      46             :     {
      47             :       Point normal_at_node(0.0);
      48       27312 :       const Node & node = elem->node_ref(i);
      49             : 
      50      176832 :       for (const auto & node_neigh_elem_id : _node_to_elem_map[node.id()])
      51             :       {
      52      149520 :         const Elem & node_neigh_elem = _cutter_mesh->elem_ref(node_neigh_elem_id);
      53             :         std::vector<Point> vertices{
      54      149520 :             node_neigh_elem.node_ref(0), node_neigh_elem.node_ref(1), node_neigh_elem.node_ref(2)};
      55      149520 :         Plane elem_plane(vertices[0], vertices[1], vertices[2]);
      56      149520 :         unsigned int j = node_neigh_elem.local_node(node.id());
      57      149520 :         Point normal_at_node_j = elem_plane.unit_normal(vertices[0]);
      58      149520 :         unsigned int m = j + 1 < 3 ? j + 1 : j + 1 - 3;
      59      149520 :         unsigned int n = j + 2 < 3 ? j + 2 : j + 2 - 3;
      60             :         Point line_1 = node_neigh_elem.node_ref(j) - node_neigh_elem.node_ref(m);
      61             :         Point line_2 = node_neigh_elem.node_ref(j) - node_neigh_elem.node_ref(n);
      62             :         Real dot = line_1 * line_2;
      63             :         Real lenSq1 = line_1 * line_1;
      64             :         Real lenSq2 = line_2 * line_2;
      65      149520 :         Real angle = std::acos(dot / std::sqrt(lenSq1 * lenSq2));
      66             :         normal_at_node += normal_at_node_j * angle;
      67      149520 :       }
      68       27312 :       normal[1 + i] = normal_at_node;
      69             :     }
      70             : 
      71       36416 :     for (unsigned int i = 0; i < elem->n_sides(); i++)
      72             :     {
      73       27312 :       std::vector<Point> vertices{elem->node_ref(0), elem->node_ref(1), elem->node_ref(2)};
      74             : 
      75       27312 :       Plane elem_plane(vertices[0], vertices[1], vertices[2]);
      76       27312 :       Point normal_at_edge = libMesh::pi * elem_plane.unit_normal(vertices[0]);
      77             : 
      78       27312 :       const Elem * neighbor = elem->neighbor_ptr(i);
      79             : 
      80       27312 :       if (neighbor != nullptr)
      81             :       {
      82             :         std::vector<Point> vertices{
      83       25648 :             neighbor->node_ref(0), neighbor->node_ref(1), neighbor->node_ref(2)};
      84             : 
      85       25648 :         Plane elem_plane(vertices[0], vertices[1], vertices[2]);
      86       25648 :         normal_at_edge += libMesh::pi * elem_plane.unit_normal(vertices[0]);
      87       25648 :       }
      88       27312 :       normal[4 + i] = normal_at_edge;
      89       27312 :     }
      90        9104 :     _pseudo_normal.insert(std::make_pair(elem->id(), normal));
      91        9136 :   }
      92          32 : }
      93             : 
      94             : Point
      95        8240 : InterfaceMeshCut3DUserObject::nodeNormal(const unsigned int & node_id)
      96             : {
      97             :   Point normal(0.0);
      98             : 
      99       52496 :   for (const auto & node_neigh_elem_id : _node_to_elem_map[node_id])
     100             :   {
     101       44256 :     const auto & elem = _cutter_mesh->elem_ref(node_neigh_elem_id);
     102       44256 :     Plane elem_plane(elem.node_ref(0), elem.node_ref(1), elem.node_ref(2));
     103       44256 :     normal += elem_plane.unit_normal(elem.node_ref(0));
     104       44256 :   }
     105             : 
     106        8240 :   unsigned int num = _node_to_elem_map[node_id].size();
     107        8240 :   return normal / num;
     108             : }
     109             : 
     110             : bool
     111           0 : InterfaceMeshCut3DUserObject::cutElementByGeometry(const Elem * /*elem*/,
     112             :                                                    std::vector<Xfem::CutEdge> & /*cut_edges*/,
     113             :                                                    std::vector<Xfem::CutNode> & /*cut_nodes*/) const
     114             : {
     115           0 :   mooseError("invalid method for 3D mesh cutting");
     116             :   return false;
     117             : }
     118             : 
     119             : bool
     120         594 : InterfaceMeshCut3DUserObject::cutElementByGeometry(const Elem * elem,
     121             :                                                    std::vector<Xfem::CutFace> & cut_faces) const
     122             : {
     123             :   mooseAssert(elem->dim() == 3, "Dimension of element to be cut must be 3");
     124             : 
     125             :   bool elem_cut = false;
     126             : 
     127        4158 :   for (unsigned int i = 0; i < elem->n_sides(); ++i)
     128             :   {
     129             :     // This returns the lowest-order type of side.
     130        3564 :     std::unique_ptr<const Elem> curr_side = elem->side_ptr(i);
     131             : 
     132             :     mooseAssert(curr_side->dim() == 2, "Side dimension must be 2");
     133             : 
     134        3564 :     unsigned int n_edges = curr_side->n_sides();
     135             : 
     136             :     std::vector<unsigned int> cut_edges;
     137             :     std::vector<Real> cut_pos;
     138             : 
     139       17820 :     for (unsigned int j = 0; j < n_edges; j++)
     140             :     {
     141             :       // This returns the lowest-order type of side.
     142       14256 :       std::unique_ptr<const Elem> curr_edge = curr_side->side_ptr(j);
     143       14256 :       if (curr_edge->type() != EDGE2)
     144           0 :         mooseError("In cutElementByGeometry face edge must be EDGE2, but type is: ",
     145           0 :                    libMesh::Utility::enum_to_string(curr_edge->type()),
     146             :                    " base element type is: ",
     147           0 :                    libMesh::Utility::enum_to_string(elem->type()));
     148             :       const Node * node1 = curr_edge->node_ptr(0);
     149             :       const Node * node2 = curr_edge->node_ptr(1);
     150             : 
     151     8449056 :       for (const auto & cut_elem : _cutter_mesh->element_ptr_range())
     152             :       {
     153             :         std::vector<Point> vertices;
     154             : 
     155    16841088 :         for (auto & node : cut_elem->node_ref_range())
     156             :         {
     157    12630816 :           Point & this_point = node;
     158    12630816 :           vertices.push_back(this_point);
     159             :         }
     160             : 
     161             :         Point intersection;
     162     4210272 :         if (Xfem::intersectWithEdge(*node1, *node2, vertices, intersection) &&
     163     4210272 :             std::find(cut_edges.begin(), cut_edges.end(), j) == cut_edges.end())
     164             :         {
     165         672 :           cut_edges.push_back(j);
     166         672 :           cut_pos.push_back(Xfem::getRelativePosition(*node1, *node2, intersection));
     167             :         }
     168     4224528 :       }
     169       14256 :     }
     170             : 
     171             :     // if two edges of an element are cut, it is considered as an element being cut
     172        3564 :     if (cut_edges.size() == 2)
     173             :     {
     174             :       elem_cut = true;
     175             :       Xfem::CutFace mycut;
     176         336 :       mycut._face_id = i;
     177         336 :       mycut._face_edge.push_back(cut_edges[0]);
     178         336 :       mycut._face_edge.push_back(cut_edges[1]);
     179         336 :       mycut._position.push_back(cut_pos[0]);
     180         336 :       mycut._position.push_back(cut_pos[1]);
     181         336 :       cut_faces.push_back(mycut);
     182             :     }
     183        3564 :   }
     184             : 
     185         594 :   return elem_cut;
     186             : }
     187             : 
     188             : bool
     189           0 : InterfaceMeshCut3DUserObject::cutFragmentByGeometry(
     190             :     std::vector<std::vector<Point>> & /*frag_edges*/,
     191             :     std::vector<Xfem::CutEdge> & /*cut_edges*/) const
     192             : {
     193           0 :   mooseError("invalid method for 3D mesh cutting");
     194             :   return false;
     195             : }
     196             : 
     197             : bool
     198           0 : InterfaceMeshCut3DUserObject::cutFragmentByGeometry(
     199             :     std::vector<std::vector<Point>> & /*frag_faces*/,
     200             :     std::vector<Xfem::CutFace> & /*cut_faces*/) const
     201             : {
     202           0 :   mooseError("cutFragmentByGeometry not yet implemented for 3D mesh cutting");
     203             :   return false;
     204             : }
     205             : 
     206             : Real
     207      110712 : InterfaceMeshCut3DUserObject::calculateSignedDistance(Point p) const
     208             : {
     209             :   std::vector<Real> distance;
     210             :   Real min_dist = std::numeric_limits<Real>::max();
     211    67443528 :   for (const auto & cut_elem : _cutter_mesh->element_ptr_range())
     212             :   {
     213             :     std::vector<Point> vertices{
     214    33611052 :         cut_elem->node_ref(0), cut_elem->node_ref(1), cut_elem->node_ref(2)};
     215             :     unsigned int region;
     216             :     Point xp;
     217    33611052 :     Real dist = Xfem::pointTriangleDistance(
     218    33611052 :         p, cut_elem->node_ref(0), cut_elem->node_ref(1), cut_elem->node_ref(2), xp, region);
     219             : 
     220    33611052 :     distance.push_back(std::abs(dist));
     221             : 
     222    33611052 :     if (dist < std::abs(min_dist))
     223             :     {
     224             :       min_dist = dist;
     225     2558078 :       Point normal = (_pseudo_normal.find(cut_elem->id())->second)[region];
     226     1279039 :       if (normal * (p - xp) < 0.0)
     227      705150 :         min_dist *= -1.0;
     228             :     }
     229    33721764 :   }
     230      110712 :   std::sort(distance.begin(), distance.end());
     231             :   Real sum_dist = 0.0;
     232      221424 :   for (std::vector<Real>::iterator it = distance.begin(); it != distance.begin() + 1; ++it)
     233      110712 :     sum_dist += *it;
     234             : 
     235      110712 :   if (min_dist < 0.0)
     236       31122 :     return -sum_dist / 1.0;
     237             :   else
     238             :     return sum_dist / 1.0;
     239      110712 : }