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             : // MOOSE includes
      11             : #include "MeshBaseDiagnosticsUtils.h"
      12             : #include "ConsoleStream.h"
      13             : #include "MooseError.h"
      14             : 
      15             : #include "libmesh/elem.h"
      16             : #include "libmesh/node.h"
      17             : #include "libmesh/mesh_base.h"
      18             : #include "libmesh/point_locator_base.h"
      19             : 
      20             : namespace MeshBaseDiagnosticsUtils
      21             : {
      22             : void
      23         192 : checkNonConformalMesh(const std::unique_ptr<MeshBase> & mesh,
      24             :                       const ConsoleStream & console,
      25             :                       const unsigned int num_outputs,
      26             :                       const Real conformality_tol,
      27             :                       unsigned int & num_nonconformal_nodes)
      28             : {
      29         192 :   num_nonconformal_nodes = 0;
      30         192 :   auto pl = mesh->sub_point_locator();
      31         192 :   pl->set_close_to_point_tol(conformality_tol);
      32             : 
      33         192 :   if (!mesh->is_serial())
      34           0 :     mooseError("Only serialized/replicated meshes are supported");
      35             : 
      36             :   // loop on nodes, assumes a replicated mesh
      37       34560 :   for (auto & node : mesh->node_ptr_range())
      38             :   {
      39             :     // find all the elements around this node
      40       17184 :     std::set<const Elem *> elements;
      41       17184 :     (*pl)(*node, elements);
      42             : 
      43             :     // loop through the set of elements near this node
      44       88136 :     for (auto & elem : elements)
      45             :     {
      46             :       // If the node is not part of this element's nodes, it is a
      47             :       // case of non-conformality
      48       70952 :       bool found_conformal = false;
      49             : 
      50      306344 :       for (auto & elem_node : elem->node_ref_range())
      51             :       {
      52      306272 :         if (*node == elem_node)
      53             :         {
      54       70880 :           found_conformal = true;
      55       70880 :           break;
      56             :         }
      57             :       }
      58       70952 :       if (!found_conformal)
      59             :       {
      60          72 :         num_nonconformal_nodes++;
      61          72 :         if (num_nonconformal_nodes < num_outputs)
      62          72 :           console << "Non-conformality detected at  : " << *node << std::endl;
      63           0 :         else if (num_nonconformal_nodes == num_outputs)
      64           0 :           console << "Maximum output reached, log is silenced" << std::endl;
      65             :       }
      66             :     }
      67       17376 :   }
      68         192 :   pl->unset_close_to_point_tol();
      69         192 : }
      70             : 
      71             : bool
      72      463909 : checkFirstOrderEdgeOverlap(const Elem & edge1,
      73             :                            const Elem & edge2,
      74             :                            Point & intersection_point,
      75             :                            const Real intersection_tol)
      76             : {
      77             :   // check that the two elements are of type EDGE2
      78             :   mooseAssert(edge1.type() == EDGE2, "Elements must be of type EDGE2");
      79             :   mooseAssert(edge2.type() == EDGE2, "Elements must be of type EDGE2");
      80             :   // Get nodes from the two edges
      81      463909 :   const Point & p1 = edge1.point(0);
      82      463909 :   const Point & p2 = edge1.point(1);
      83      463909 :   const Point & p3 = edge2.point(0);
      84      463909 :   const Point & p4 = edge2.point(1);
      85             : 
      86             :   // Check that the two edges are not sharing a node
      87      463909 :   if (p1 == p3 || p1 == p4 || p2 == p3 || p2 == p4)
      88      159872 :     return false;
      89             : 
      90             :   /*
      91             :     There's a chance that they overlap. Find shortest line that connects two edges and if its length
      92             :     is close enough to 0 return true The shortest line between the two edges will be perpendicular
      93             :     to both.
      94             :   */
      95      304037 :   const auto d1343 = (p1 - p3) * (p4 - p3);
      96      304037 :   const auto d4321 = (p4 - p3) * (p2 - p1);
      97      304037 :   const auto d1321 = (p1 - p3) * (p2 - p1);
      98      304037 :   const auto d4343 = (p4 - p3) * (p4 - p3);
      99      304037 :   const auto d2121 = (p2 - p1) * (p2 - p1);
     100             : 
     101      304037 :   const auto denominator = d2121 * d4343 - d4321 * d4321;
     102      304037 :   const auto numerator = d1343 * d4321 - d1321 * d4343;
     103             : 
     104      304037 :   if (std::fabs(denominator) < intersection_tol)
     105             :     // This indicates that the two lines are parallel so they don't intersect
     106       42709 :     return false;
     107             : 
     108      261328 :   const auto mua = numerator / denominator;
     109      261328 :   const auto mub = (d1343 + (mua * d4321)) / d4343;
     110             : 
     111             :   // Use these values to solve for the two points that define the shortest line segment
     112      261328 :   const auto pa = p1 + mua * (p2 - p1);
     113      261328 :   const auto pb = p3 + mub * (p4 - p3);
     114             : 
     115             :   // This method assume the two lines are infinite. This check to make sure na and nb are part of
     116             :   // their respective line segments
     117      261328 :   if (mua < 0 || mua > 1)
     118       94696 :     return false;
     119      166632 :   if (mub < 0 || mub > 1)
     120       57688 :     return false;
     121             : 
     122             :   // Calculate distance between these two nodes
     123      108944 :   const auto distance = (pa - pb).norm();
     124      108944 :   if (distance < intersection_tol)
     125             :   {
     126          16 :     intersection_point = pa;
     127          16 :     return true;
     128             :   }
     129             :   else
     130      108928 :     return false;
     131             : }
     132             : }