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 "SurfaceSubdomainsFromAllNormalsGenerator.h"
      11             : #include "InputParameters.h"
      12             : #include "MooseMeshUtils.h"
      13             : #include "CastUniquePointer.h"
      14             : 
      15             : #include "libmesh/mesh_generation.h"
      16             : #include "libmesh/mesh.h"
      17             : #include "libmesh/string_to_enum.h"
      18             : #include "libmesh/elem.h"
      19             : 
      20             : #include <typeinfo>
      21             : 
      22             : registerMooseObject("MooseApp", SurfaceSubdomainsFromAllNormalsGenerator);
      23             : 
      24             : InputParameters
      25        3237 : SurfaceSubdomainsFromAllNormalsGenerator::validParams()
      26             : {
      27        3237 :   InputParameters params = SurfaceMeshGeneratorBase::validParams();
      28             : 
      29       19422 :   params.renameParam("max_paint_size_centroids", "max_subdomain_size_centroids", "");
      30        9711 :   params.addParam<bool>(
      31             :       "contiguous_assignments_only",
      32        6474 :       false,
      33             :       "Whether to only group elements in a subdomain using the 'flooding' algorithm. "
      34             :       "Note: this also sets 'flood_elements_once' to true if the user did not explicitly pass "
      35             :       "false.");
      36             : 
      37             :   // Post flooding operations/cleanup
      38        9711 :   params.addParam<bool>("select_max_neighbor_element_subdomains",
      39        6474 :                         false,
      40             :                         "Whether to perform a final subdomain assignment phase where each element "
      41             :                         "is assigned to the subdomain "
      42             :                         "that holds the most neighbors of the element with a similar normal.");
      43        9711 :   params.addParam<bool>(
      44             :       "separate_elements_connected_by_a_single_node",
      45        6474 :       false,
      46             :       "Whether to perform a final subdomain assignment phase where any element that is not "
      47             :       "connected to other elements from the same subdomain is re-assigned to either the subdomain "
      48             :       "that all side neighbor elements belong to, or a new subdomain if there is no such subdomain "
      49             :       "including all its side neighbors. Note that normals are not checked in this phase, and the "
      50             :       "subdomain size criterion is assumed to be already met.");
      51             : 
      52       12948 :   params.addParamNamesToGroup("contiguous_assignments_only", "Other flooding");
      53             :   // NOTE: these post-flooding / post-treatment operations could be moved to a base class.
      54             :   // They could also become their own mesh generator!
      55       12948 :   params.addParamNamesToGroup(
      56             :       "separate_elements_connected_by_a_single_node select_max_neighbor_element_subdomains",
      57             :       "Post flooding subdomain re-assignments");
      58             : 
      59             :   // There can be many of them, and we don't control the number in this generator
      60        3237 :   params.suppressParameter<std::vector<SubdomainName>>("new_subdomain");
      61             :   // Using normals is the base principle of this mesh generator
      62        9711 :   params.addPrivateParam<bool>("_using_normal", true);
      63             : 
      64        3237 :   params.addClassDescription(
      65             :       "Adds subdomains to surface (2D) elements in the (3D) mesh based on unique normals.");
      66        3237 :   return params;
      67           0 : }
      68             : 
      69          88 : SurfaceSubdomainsFromAllNormalsGenerator::SurfaceSubdomainsFromAllNormalsGenerator(
      70          88 :     const InputParameters & parameters)
      71             :   : SurfaceMeshGeneratorBase(parameters),
      72         176 :     _contiguous_assignments_only(getParam<bool>("contiguous_assignments_only"))
      73             : {
      74         104 :   if (_contiguous_assignments_only && !isParamSetByUser("flood_elements_once"))
      75           0 :     _flood_only_once = true;
      76          88 : }
      77             : 
      78             : std::unique_ptr<MeshBase>
      79          88 : SurfaceSubdomainsFromAllNormalsGenerator::generate()
      80             : {
      81          88 :   std::unique_ptr<MeshBase> mesh = std::move(_input);
      82          88 :   if (!mesh->is_serial())
      83           0 :     mooseError(
      84             :         "SurfaceSubdomainsFromAllNormalsGenerator is not implemented for non serialized meshes");
      85          88 :   setup(*mesh);
      86             : 
      87          88 :   unsigned int num_neighborless = 0;
      88         264 :   const bool user_specified_normal = isParamSetByUser("normal");
      89             : 
      90             :   // Make sure neighbors are known
      91          88 :   if (!mesh->preparation().has_neighbor_ptrs)
      92          80 :     mesh->find_neighbors();
      93             : 
      94             :   // We'll need to loop over all of the elements to find ones that match this normal.
      95             :   // We can't rely on flood catching them all in one go. We have to flood from multiple elements
      96             :   // in case there are disconnected groups of elements.
      97             :   // Depending on user parameters, we may be "flooding" multiple times the same elements
      98        2392 :   for (auto & elem : mesh->element_ptr_range())
      99             :   {
     100             :     // Nothing to do with edges or 3D elements
     101        2304 :     if (elem->dim() != 2)
     102         320 :       continue;
     103             :     // Likely an issue with the mesh
     104        2304 :     if (_contiguous_assignments_only && elem->n_neighbors() == 0)
     105             :     {
     106           0 :       num_neighborless++;
     107           0 :       mooseWarning("Element '" + std::to_string(elem->id()) + "' has no neighbors");
     108             :     }
     109             : 
     110             :     // Check if element should be used to paint from
     111        2304 :     if (_included_subdomain_ids.size() &&
     112           0 :         std::find(_included_subdomain_ids.begin(),
     113             :                   _included_subdomain_ids.end(),
     114        2304 :                   elem->subdomain_id()) == _included_subdomain_ids.end())
     115           0 :       continue;
     116             : 
     117             :     // Compute the normal
     118        2304 :     const auto elem_normal = get2DElemNormal(elem);
     119             : 
     120             :     // Four options for which subdomain to paint with:
     121             :     // 1) See if we've seen this normal before (linear search), within some tolerance
     122             :     // NOTE: these elements might not be neighbors of the current element,
     123             :     // so the patch may not be contiguous
     124        2304 :     const std::map<SubdomainID, RealVectorValue>::value_type * item = nullptr;
     125        2304 :     bool sub_id_found = false;
     126        2304 :     if (!_contiguous_assignments_only)
     127        7264 :       for (const auto & id_pair : _subdomain_to_normal_map)
     128        6368 :         if (elementSatisfiesRequirements(
     129        6368 :                 elem, id_pair.second, *_subdomain_to_starting_elem[id_pair.first], elem_normal))
     130             :         {
     131        1024 :           sub_id_found = true;
     132        1024 :           item = &id_pair;
     133        1024 :           break;
     134             :         }
     135             :     // 2) Look at the neighbors, if a majority of them have the same subdomain and a similar normal,
     136             :     // use that subdomain to paint elements with
     137             :     // NOTE: compatible with 'contiguous_assignments_only' option
     138             :     // NOTE: this overrides the result of the previous search
     139             :     SubdomainID neighbor_majority_sub_id;
     140        2304 :     if (_check_painted_neighor_normals)
     141             :     {
     142           0 :       std::map<SubdomainID, unsigned int> sub_id_neighbors;
     143             :       // Try to flood from each side with the same subdomain
     144             :       // Look for the neighbor subdomain id with the most neighbors
     145           0 :       for (const auto neighbor : make_range(elem->n_sides()))
     146           0 :         if (elem->neighbor_ptr(neighbor) &&
     147           0 :             _acted_upon_once.find(elem->neighbor_ptr(neighbor)) != _acted_upon_once.end() &&
     148           0 :             elementSatisfiesRequirements(elem,
     149           0 :                                          get2DElemNormal(elem->neighbor_ptr(neighbor)),
     150           0 :                                          *elem->neighbor_ptr(neighbor),
     151             :                                          elem_normal))
     152             :         {
     153           0 :           sub_id_found = true;
     154           0 :           sub_id_neighbors[elem->neighbor_ptr(neighbor)->subdomain_id()]++;
     155             :         }
     156             : 
     157           0 :       unsigned int max_of_subid = 0;
     158           0 :       for (const auto & [key, item] : sub_id_neighbors)
     159           0 :         if (item >= max_of_subid)
     160             :         {
     161           0 :           max_of_subid = item;
     162           0 :           neighbor_majority_sub_id = key;
     163             :         }
     164             : 
     165             :       // Note: the max distance from starting element of the subdomain flooding to this
     166             :       // element should be checked in the call to flood()
     167           0 :     }
     168             :     // 3) Flood only with the selected fixed normal
     169             :     // Note: this steers away from an "FromAllNormals" generator, it's more like a "FromNormal"
     170             :     // generator capability.
     171        2304 :     if (user_specified_normal && !elementSatisfiesRequirements(elem, _normal, *elem, elem_normal))
     172         320 :       continue;
     173             :     // 4) Flood with a new subdomain every time ('contiguous_assignments_only' option)
     174             : 
     175             :     // Finalize flooding parameter selection
     176        1984 :     Elem * starting_element = elem;
     177             :     subdomain_id_type flooding_sub_id;
     178        1984 :     Point flooding_normal = elem_normal;
     179        1984 :     if (item)
     180             :     {
     181        1024 :       starting_element = _subdomain_to_starting_elem[item->first];
     182        1024 :       flooding_sub_id = item->first;
     183        1024 :       flooding_normal = item->second;
     184             :     }
     185         960 :     else if (sub_id_found)
     186             :     {
     187           0 :       starting_element = _subdomain_to_starting_elem[neighbor_majority_sub_id];
     188           0 :       flooding_sub_id = neighbor_majority_sub_id;
     189             :     }
     190             :     else
     191             :     {
     192         960 :       flooding_sub_id = MooseMeshUtils::getNextFreeSubdomainID(*mesh);
     193         960 :       _subdomain_to_normal_map[flooding_sub_id] = elem_normal;
     194         960 :       _subdomain_to_starting_elem[flooding_sub_id] = elem;
     195             :     }
     196             :     // User input of a normal takes precedence
     197        1984 :     if (user_specified_normal)
     198             :     {
     199          64 :       flooding_normal = _normal;
     200          64 :       _subdomain_to_normal_map[flooding_sub_id] = flooding_normal;
     201             :     }
     202             : 
     203             :     // Flood with the previously created subdomains and normals
     204        1984 :     flood(elem, flooding_normal, *starting_element, flooding_sub_id, *mesh);
     205          88 :   }
     206             : 
     207             :   // Check all elements once with only neighbor averaging
     208         264 :   if (getParam<bool>("select_max_neighbor_element_subdomains"))
     209         200 :     for (auto & elem : mesh->element_ptr_range())
     210             :     {
     211             :       // Nothing to do with edges or 3D elements
     212         192 :       if (elem->dim() != 2)
     213           0 :         continue;
     214             : 
     215             :       // Compute the normal
     216         192 :       const auto normal = get2DElemNormal(elem);
     217             : 
     218         192 :       bool sub_id_found = false;
     219         192 :       SubdomainID sub_id = 0;
     220         192 :       std::map<SubdomainID, unsigned int> sub_id_neighbors;
     221             : 
     222             :       // NOTE: we use side neighbors instead of point neighbors as side neighbors
     223             :       // are always connected by 2 nodes to the rest of the subdomains, and elements
     224             :       // connected by a single node to a subdomain cause issues when triangulating
     225             :       // that subdomain.
     226             : 
     227             :       // Try to flood from each side with the same subdomain
     228             :       // Look for the neighbor subdomain id with the most neighbors
     229         960 :       for (const auto neighbor_i : make_range(elem->n_sides()))
     230             :       {
     231         768 :         const auto neighbor = elem->neighbor_ptr(neighbor_i);
     232             :         // NOTE: to avoid exceeding the patch size, we pass the element that was used to start
     233             :         // painting the patch to check the 'distance'/'patch size' criterion
     234        1536 :         if (neighbor &&
     235         768 :             elementSatisfiesRequirements(elem,
     236         768 :                                          get2DElemNormal(neighbor),
     237         768 :                                          *_subdomain_to_starting_elem[neighbor->subdomain_id()],
     238             :                                          normal))
     239             :         {
     240         656 :           sub_id_found = true;
     241         656 :           sub_id_neighbors[neighbor->subdomain_id()]++;
     242             :         }
     243             :       }
     244             : 
     245         192 :       unsigned int max_of_subid = 0;
     246         464 :       for (const auto & [key, item] : sub_id_neighbors)
     247         272 :         if (item >= max_of_subid)
     248             :         {
     249         256 :           max_of_subid = item;
     250         256 :           sub_id = key;
     251             :         }
     252         192 :       if (sub_id_found)
     253         192 :         elem->subdomain_id() = sub_id;
     254         200 :     }
     255             : 
     256         264 :   if (getParam<bool>("separate_elements_connected_by_a_single_node"))
     257         200 :     for (auto & elem : mesh->element_ptr_range())
     258             :     {
     259             :       // Nothing to do with edges or 3D elements
     260         192 :       if (elem->dim() != 2)
     261           0 :         continue;
     262             : 
     263         192 :       bool connected_to_a_neighbor = false;
     264         960 :       for (const auto neighbor : make_range(elem->n_sides()))
     265        1536 :         if (elem->neighbor_ptr(neighbor) &&
     266         768 :             elem->subdomain_id() == elem->neighbor_ptr(neighbor)->subdomain_id())
     267         400 :           connected_to_a_neighbor = true;
     268             : 
     269         192 :       if (!connected_to_a_neighbor)
     270             :       {
     271          16 :         bool same_subdomain = true;
     272          16 :         subdomain_id_type common_sub = std::numeric_limits<subdomain_id_type>::max();
     273             : 
     274             :         // If all side-neighbors have the same subdomain, use that subdomain instead
     275             :         // NOTE: we don't check the normal criteria here
     276          80 :         for (const auto neighbor : make_range(elem->n_sides()))
     277          64 :           if (elem->neighbor_ptr(neighbor))
     278             :           {
     279          64 :             if (common_sub == std::numeric_limits<subdomain_id_type>::max())
     280          16 :               common_sub = elem->neighbor_ptr(neighbor)->subdomain_id();
     281          48 :             else if (common_sub != elem->neighbor_ptr(neighbor)->subdomain_id())
     282          32 :               same_subdomain = false;
     283             :           }
     284             : 
     285          16 :         if (same_subdomain && (common_sub != std::numeric_limits<subdomain_id_type>::max()))
     286           0 :           elem->subdomain_id() = common_sub;
     287             :         else
     288          16 :           elem->subdomain_id() = MooseMeshUtils::getNextFreeSubdomainID(*mesh);
     289             :       }
     290           8 :     }
     291             : 
     292          88 :   if (_contiguous_assignments_only && num_neighborless)
     293           0 :     mooseWarning("Several subdomains were created for neighborless elements: " +
     294           0 :                  std::to_string(num_neighborless));
     295             : 
     296             :   // subdomain assignments have changed
     297          88 :   mesh->unset_is_prepared();
     298         176 :   return dynamic_pointer_cast<MeshBase>(mesh);
     299          88 : }
     300             : 
     301             : void
     302        3656 : SurfaceSubdomainsFromAllNormalsGenerator::actOnElem(Elem * const elem,
     303             :                                                     const Point &,
     304             :                                                     const subdomain_id_type & sub_id,
     305             :                                                     MeshBase &)
     306             : {
     307        3656 :   elem->subdomain_id() = sub_id;
     308        3656 : }