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 "MeshCoarseningUtils.h"
      11             : #include "Conversion.h"
      12             : #include "MooseError.h"
      13             : 
      14             : #include "libmesh/enum_elem_type.h"
      15             : #include "libmesh/remote_elem.h"
      16             : 
      17             : using namespace libMesh;
      18             : 
      19             : namespace MeshCoarseningUtils
      20             : {
      21             : bool
      22       12829 : getFineElementsFromInteriorNode(const libMesh::Node & interior_node,
      23             :                                 const libMesh::Node & reference_node,
      24             :                                 const libMesh::Elem & fine_elem,
      25             :                                 std::vector<const libMesh::Node *> & tentative_coarse_nodes,
      26             :                                 std::set<const libMesh::Elem *> & fine_elements)
      27             : {
      28       12829 :   const auto elem_type = fine_elem.type();
      29             : 
      30             :   // Add point neighbors of interior node to list of potentially refined elements
      31             :   // NOTE: we could potentially replace this with a simple call to point_neighbors
      32             :   // on a fine element with the interior node. It's not clear which approach is more
      33             :   // resilient to meshes with slits from discarded adaptivity information
      34       12829 :   fine_elements.insert(&fine_elem);
      35       86447 :   for (const auto neigh : fine_elem.neighbor_ptr_range())
      36             :   {
      37       73618 :     if (!neigh || neigh == libMesh::remote_elem)
      38       27082 :       continue;
      39       46536 :     const auto node_index = neigh->get_node_index(&interior_node);
      40       46536 :     if (node_index != libMesh::invalid_uint && neigh->is_vertex(node_index))
      41             :     {
      42             :       // Get the neighbor's neighbors, to catch point non-side neighbors
      43             :       // This is needed in 2D to get all quad neighbors
      44       32174 :       fine_elements.insert(neigh);
      45      219586 :       for (const auto neigh_two : neigh->neighbor_ptr_range())
      46             :       {
      47      187412 :         if (!neigh_two || neigh_two == libMesh::remote_elem)
      48       59972 :           continue;
      49      127440 :         const auto node_index_2 = neigh_two->get_node_index(&interior_node);
      50      127440 :         if (node_index_2 != libMesh::invalid_uint && neigh_two->is_vertex(node_index_2))
      51             :         {
      52             :           // Get the neighbor's neighbors
      53       85192 :           fine_elements.insert(neigh_two);
      54             : 
      55             :           // Get the neighbor's neighbors' neighbors, to catch point non-side neighbors
      56             :           // This is needed for 3D to get all hex neighbors
      57      586106 :           for (const auto neigh_three : neigh_two->neighbor_ptr_range())
      58             :           {
      59      500914 :             if (!neigh_three || neigh_three == libMesh::remote_elem)
      60      149158 :               continue;
      61      351756 :             const auto node_index_3 = neigh_three->get_node_index(&interior_node);
      62      351756 :             if (node_index_3 != libMesh::invalid_uint && neigh_three->is_vertex(node_index_3))
      63      233969 :               fine_elements.insert(neigh_three);
      64             :           }
      65             :         }
      66             :       }
      67             :     }
      68             :   }
      69             : 
      70             :   // If the fine elements are not all of the same type, we do not know how to get the opposite node
      71             :   // of the interior node in the fine elements
      72       91532 :   for (auto elem : fine_elements)
      73       78703 :     if (elem && fine_elem.type() != elem_type)
      74           0 :       return false;
      75             : 
      76       12829 :   if (elem_type == libMesh::QUAD4 || elem_type == libMesh::QUAD8 || elem_type == libMesh::QUAD9)
      77             :   {
      78             :     // We need 4 elements around the interior node
      79        1678 :     if (fine_elements.size() != 4)
      80         531 :       return false;
      81             : 
      82             :     // We need to order the fine elements so when we get the coarse element nodes they form
      83             :     // a non-twisted element
      84        1147 :     tentative_coarse_nodes.resize(4);
      85             : 
      86             :     // The exterior nodes are the opposite nodes of the interior_node!
      87        1147 :     unsigned int neighbor_i = 0;
      88        5735 :     for (auto neighbor : fine_elements)
      89             :     {
      90        4588 :       const auto interior_node_number = neighbor->get_node_index(&interior_node);
      91        4588 :       unsigned int opposite_node_index = (interior_node_number + 2) % 4;
      92             : 
      93        4588 :       tentative_coarse_nodes[neighbor_i++] = neighbor->node_ptr(opposite_node_index);
      94             :     }
      95             : 
      96             :     // Re-order nodes so that they will form a decent quad
      97             :     libMesh::Point axis =
      98        1147 :         (fine_elem.vertex_average() - interior_node).cross(interior_node - reference_node);
      99        1147 :     reorderNodes(tentative_coarse_nodes, interior_node, reference_node, axis);
     100        1147 :     return true;
     101             :   }
     102             :   // For hexes, similar strategy but we need to pick 4 nodes to form a side, then 4 other nodes
     103             :   // facing those initial nodes
     104       11151 :   else if (elem_type == libMesh::HEX8)
     105             :   {
     106             :     // We need 8 elements around the interior node
     107       11151 :     if (fine_elements.size() != 8)
     108        4032 :       return false;
     109             : 
     110        7119 :     tentative_coarse_nodes.resize(4);
     111             : 
     112             :     // Pick a node (mid-face for the coarse element) to form the base
     113             :     // We must use the same element reproducibly, despite the pointers being ordered in the set
     114             :     // We use the element id to choose the same element consistently
     115        7119 :     const Elem * one_fine_elem = nullptr;
     116        7119 :     unsigned int max_id = 0;
     117       64071 :     for (const auto elem_ptr : fine_elements)
     118       56952 :       if (elem_ptr->id() > max_id)
     119             :       {
     120       50003 :         max_id = elem_ptr->id();
     121       50003 :         one_fine_elem = elem_ptr;
     122             :       }
     123        7119 :     const auto interior_node_index = one_fine_elem->get_node_index(&interior_node);
     124             : 
     125             :     // Find any side which contains the interior node
     126        7119 :     unsigned int an_interior_node_side = 0;
     127        7119 :     for (const auto s : make_range(one_fine_elem->n_sides()))
     128        7119 :       if (one_fine_elem->is_node_on_side(interior_node_index, s))
     129             :       {
     130        7119 :         an_interior_node_side = s;
     131        7119 :         break;
     132             :       }
     133             :     // A node near a face of the coarse element seek is on the same side, but opposite from the
     134             :     // interior node
     135             :     const auto center_face_node_index =
     136        7119 :         one_fine_elem->opposite_node(interior_node_index, an_interior_node_side);
     137        7119 :     const auto center_face_node = one_fine_elem->node_ptr(center_face_node_index);
     138             : 
     139             :     // We gather the coarse element nodes from the fine elements that share the center face node we
     140             :     // just selected
     141        7119 :     unsigned int neighbor_i = 0;
     142        7119 :     std::vector<const libMesh::Elem *> other_fine_elems;
     143       64071 :     for (auto neighbor : fine_elements)
     144             :     {
     145       56952 :       if (neighbor->get_node_index(center_face_node) == libMesh::invalid_uint)
     146             :       {
     147       28476 :         other_fine_elems.push_back(neighbor);
     148       28476 :         continue;
     149             :       }
     150             :       // The coarse element node is the opposite nodes of the interior_node in a fine element
     151       28476 :       const auto interior_node_number = neighbor->get_node_index(&interior_node);
     152             :       unsigned int opposite_node_index =
     153       28476 :           getOppositeNodeIndex(neighbor->type(), interior_node_number);
     154             : 
     155       28476 :       tentative_coarse_nodes[neighbor_i++] = neighbor->node_ptr(opposite_node_index);
     156             :     }
     157             : 
     158             :     // Center face node was not shared with 4 elements
     159             :     // We could try again on any of the 5 other coarse element faces but we don't insist for now
     160        7119 :     if (neighbor_i != 4 || other_fine_elems.size() != 4)
     161           0 :       return false;
     162             : 
     163             :     // Sort the coarse nodes so we are reproducibly picking the same ordering of nodes
     164       43038 :     auto cmp_node = [](const Node * a, const Node * b) { return a->id() < b->id(); };
     165        7119 :     std::sort(tentative_coarse_nodes.begin(), tentative_coarse_nodes.end(), cmp_node);
     166             : 
     167             :     // Re-order nodes so that they will form a decent quad
     168             :     // Pick the reference node for the rotation frame as the face center
     169        7119 :     const libMesh::Point clock_start = *tentative_coarse_nodes[0];
     170        7119 :     libMesh::Point axis = interior_node - *center_face_node;
     171        7119 :     reorderNodes(tentative_coarse_nodes, *center_face_node, clock_start, axis);
     172             : 
     173             :     // Look through the 4 other fine elements to finish the coarse hex element nodes
     174       35595 :     for (const auto coarse_node_index : make_range(4))
     175             :     {
     176             :       // Find the fine element containing each coarse node already found & ordered
     177       28476 :       const Elem * fine_elem = nullptr;
     178      161941 :       for (auto elem : fine_elements)
     179      161941 :         if (elem->get_node_index(tentative_coarse_nodes[coarse_node_index]) !=
     180             :             libMesh::invalid_uint)
     181             :         {
     182       28476 :           fine_elem = elem;
     183       28476 :           break;
     184             :         }
     185             :       mooseAssert(fine_elem, "Search for fine element should have worked");
     186             : 
     187             :       // Find the other fine element opposite the element containing the node (they share a side)
     188       28476 :       const Elem * fine_neighbor = nullptr;
     189      142380 :       for (auto neighbor : other_fine_elems)
     190             :         // Side neighbor. This requires the mesh to have correct neighbors
     191             :         // For meshes with lost AMR information, this wont work
     192      113904 :         if (neighbor->which_neighbor_am_i(fine_elem) != libMesh::invalid_uint)
     193             :         {
     194       28476 :           if (fine_neighbor)
     195           0 :             mooseError("Found two neighbors");
     196       28476 :           fine_neighbor = neighbor;
     197             :         }
     198             :       // the fine element in the base of the coarse hex is not a neighbor to any element
     199             :       // in the top part. The mesh is probably slit in the middle of the potential coarse hex
     200             :       // element. We wont support this for now.
     201       28476 :       if (!fine_neighbor)
     202           0 :         return false;
     203             : 
     204             :       // Get the coarse node, opposite the interior node in that fine element
     205       28476 :       const auto interior_node_index_neighbor = fine_neighbor->get_node_index(&interior_node);
     206       28476 :       tentative_coarse_nodes.push_back(fine_neighbor->node_ptr(
     207       28476 :           getOppositeNodeIndex(fine_neighbor->type(), interior_node_index_neighbor)));
     208             :     }
     209             : 
     210             :     // Found 8 fine elements and 8 coarse element nodes as expected
     211        7119 :     if (tentative_coarse_nodes.size() == 8)
     212        7119 :       return true;
     213             :     else
     214           0 :       return false;
     215        7119 :   }
     216             :   else
     217           0 :     mooseError("Not implemented for element type " + Moose::stringify(elem_type));
     218             : }
     219             : 
     220             : void
     221        8354 : reorderNodes(std::vector<const libMesh::Node *> & nodes,
     222             :              const libMesh::Point & origin,
     223             :              const libMesh::Point & clock_start,
     224             :              libMesh::Point & axis)
     225             : {
     226             :   mooseAssert(axis.norm() != 0, "Invalid rotation axis when ordering nodes");
     227             :   mooseAssert(origin != clock_start, "Invalid starting direction when ordering nodes");
     228             : 
     229             :   // We'll need to order the coarse nodes based on the clock-wise order of the elements
     230             :   // Define a frame in which to compute the angles of the fine elements centers
     231             :   // angle 0 is the [interior node, non-conformal node] vertex
     232        8354 :   auto start_clock = origin - clock_start;
     233        8354 :   start_clock /= start_clock.norm();
     234        8354 :   axis /= axis.norm();
     235             : 
     236        8354 :   std::vector<std::pair<unsigned int, libMesh::Real>> nodes_angles(nodes.size());
     237       41770 :   for (const auto angle_i : index_range(nodes))
     238             :   {
     239             :     mooseAssert(nodes[angle_i], "Nodes cant be nullptr");
     240       33416 :     auto vec = *nodes[angle_i] - origin;
     241       33416 :     vec /= vec.norm();
     242       33416 :     const auto angle = atan2(vec.cross(start_clock) * axis, vec * start_clock);
     243       33416 :     nodes_angles[angle_i] = std::make_pair(angle_i, angle);
     244             :   }
     245             : 
     246             :   // sort by angle, so it goes around the interior node
     247        8354 :   std::sort(nodes_angles.begin(),
     248             :             nodes_angles.end(),
     249       62382 :             [](auto & left, auto & right) { return left.second < right.second; });
     250             : 
     251             :   // Re-sort the nodes based on their angle
     252        8354 :   std::vector<const libMesh::Node *> new_nodes(nodes.size());
     253       41770 :   for (const auto & old_index : index_range(nodes))
     254       33416 :     new_nodes[old_index] = nodes[nodes_angles[old_index].first];
     255       41770 :   for (const auto & index : index_range(nodes))
     256       33416 :     nodes[index] = new_nodes[index];
     257        8354 : }
     258             : 
     259             : unsigned int
     260       92835 : getOppositeNodeIndex(libMesh::ElemType elem_type, unsigned int node_index)
     261             : {
     262       92835 :   switch (elem_type)
     263             :   {
     264        3159 :     case QUAD4:
     265        3159 :       return (node_index + 2) % 4;
     266       89676 :     case HEX8:
     267             :     {
     268             :       mooseAssert(node_index < 8, "Node index too high: " + std::to_string(node_index));
     269       89676 :       return std::vector<unsigned int>({6, 7, 4, 5, 2, 3, 0, 1})[node_index];
     270             :     }
     271           0 :     default:
     272           0 :       mooseError("Unsupported element type for retrieving the opposite node");
     273             :   }
     274             : }
     275             : }