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 "ElementGroupCentroidPositions.h"
      11             : 
      12             : #include "libmesh/parallel_algebra.h"
      13             : 
      14             : registerMooseObject("MooseApp", ElementGroupCentroidPositions);
      15             : 
      16             : InputParameters
      17       14409 : ElementGroupCentroidPositions::validParams()
      18             : {
      19       14409 :   InputParameters params = Positions::validParams();
      20       14409 :   params += BlockRestrictable::validParams();
      21       14409 :   params.addClassDescription("Gets the Positions of the centroid of groups of elements. "
      22             :                              "Groups may be defined using subdomains or element extra ids.");
      23             : 
      24             :   // To enable extra element ID groups
      25       14409 :   params.addParam<MooseEnum>("grouping_type", groupTypeEnum(), "Type of group of elements");
      26       14409 :   params.addParam<std::vector<ExtraElementIDName>>("extra_id_name",
      27             :                                                    "Name(s) of the extra element ID(s) to use");
      28       14409 :   params.addParam<std::vector<std::vector<unsigned int>>>(
      29             :       "extra_id",
      30             :       "Specific ID(s), for each extra id name, for grouping elements. "
      31             :       "If empty, all *valid* ids will be used to bin");
      32             : 
      33             :   // Order already arises from the block/ID bins
      34       14409 :   params.set<bool>("auto_sort") = false;
      35             :   // Full replicated mesh is used on every rank to generate positions
      36       14409 :   params.set<bool>("auto_broadcast") = false;
      37             : 
      38       14409 :   return params;
      39           0 : }
      40             : 
      41          72 : ElementGroupCentroidPositions::ElementGroupCentroidPositions(const InputParameters & parameters)
      42             :   : Positions(parameters),
      43             :     BlockRestrictable(this),
      44         144 :     _mesh(_fe_problem.mesh()),
      45          72 :     _group_type(getParam<MooseEnum>("grouping_type"))
      46             : {
      47             :   // We are not excluding using both block restriction and extra element ids
      48          72 :   if (_group_type == "extra_id" || _group_type == "block_and_extra_id")
      49             :   {
      50          48 :     _extra_id_names = getParam<std::vector<ExtraElementIDName>>("extra_id_name");
      51         156 :     for (const auto & name : _extra_id_names)
      52         108 :       _extra_id_indices.push_back(_mesh.getMesh().get_elem_integer_index(name));
      53             : 
      54          48 :     if (isParamValid("extra_id"))
      55          48 :       _extra_id_group_indices = getParam<std::vector<std::vector<unsigned int>>>("extra_id");
      56             :     else
      57           0 :       _extra_id_group_indices.resize(_extra_id_names.size());
      58             : 
      59          48 :     if (_extra_id_group_indices.size() != _extra_id_names.size())
      60           0 :       paramError("extra_id",
      61             :                  "Number of extra id names and the indices to select must match. "
      62             :                  "If you want all indices for an extra id, use an empty vector entry");
      63             : 
      64             :     // Can only have so many groups though, considering 4D max capability
      65          48 :     if (_extra_id_indices.size() > unsigned(3 + blockRestricted()))
      66           0 :       mooseError("Positions currently supports only up to 4D storage");
      67             :   }
      68             :   else
      69             :   {
      70          24 :     if (isParamValid("extra_id_name"))
      71           0 :       paramError("extra_id_name",
      72             :                  "An extra id name was specified but elements are not grouped by extra ids");
      73          24 :     if (isParamValid("extra_ids"))
      74           0 :       paramError("extra_ids",
      75             :                  "An extra id was specified but elements are not grouped by extra ids");
      76             :   }
      77             : 
      78             :   // Insert subdomain as an extra element id to simplify computation logic
      79          72 :   if (blockRestricted() || !isParamValid("extra_id_name"))
      80             :   {
      81          60 :     _blocks_in_use = true;
      82          60 :     _extra_id_names.insert(_extra_id_names.begin(), "block");
      83          60 :     _extra_id_indices.insert(_extra_id_indices.begin(), std::numeric_limits<unsigned short>::max());
      84          60 :     _extra_id_group_indices.insert(_extra_id_group_indices.begin(), std::vector<unsigned int>());
      85             :     // Add real block restriction
      86          60 :     if (blockRestricted())
      87         132 :       for (const auto & block : blockIDs())
      88          84 :         _extra_id_group_indices[0].push_back(block);
      89             :     // Just add all blocks
      90             :     else
      91          48 :       for (const auto & block : meshBlockIDs())
      92          36 :         _extra_id_group_indices[0].push_back(block);
      93             :   }
      94             :   else
      95          12 :     _blocks_in_use = false;
      96             : 
      97             :   // Mesh is ready at construction
      98          72 :   initialize();
      99             :   // Sort the output (by XYZ) if requested
     100          72 :   finalize();
     101          72 : }
     102             : 
     103             : void
     104          72 : ElementGroupCentroidPositions::initialize()
     105             : {
     106          72 :   clearPositions();
     107             :   // By default, initialize should be called on meshChanged()
     108             : 
     109             :   // If users did not specify a value for an extra element integer, they want all the bins
     110             :   // We'll need to loop through the mesh to find the element extra ids
     111         240 :   for (const auto i : index_range(_extra_id_group_indices))
     112             :   {
     113         168 :     auto & indices = _extra_id_group_indices[i];
     114         168 :     if (indices.empty())
     115             :     {
     116          24 :       std::set<unsigned int> ids;
     117       64536 :       for (const auto & elem : _mesh.getMesh().active_local_element_ptr_range())
     118       64536 :         ids.insert(id(*elem, _extra_id_indices[i], _blocks_in_use && i == 0));
     119          24 :       _mesh.comm().set_union(ids);
     120         120 :       for (const auto & id : ids)
     121          96 :         indices.push_back(id);
     122          24 :     }
     123             :   }
     124             : 
     125             :   // Allocate some vectors holding the volumes
     126          72 :   std::vector<Real> volumes;
     127          72 :   std::vector<std::vector<Real>> volumes_2d;
     128          72 :   std::vector<std::vector<std::vector<Real>>> volumes_3d;
     129          72 :   std::vector<std::vector<std::vector<std::vector<Real>>>> volumes_4d;
     130             : 
     131             :   // Default indexing: blocks, extra_id_name_1, extra_id_name_2 ...
     132             :   // Allocate vectors holding the positions
     133          72 :   if (_extra_id_names.size() == 1)
     134             :   {
     135          24 :     _positions.resize(_extra_id_group_indices[0].size());
     136          24 :     volumes.resize(_positions.size());
     137             :   }
     138          48 :   else if (_extra_id_names.size() == 2)
     139             :   {
     140          12 :     _positions_2d.resize(_extra_id_group_indices[0].size());
     141          12 :     volumes_2d.resize(_positions_2d.size());
     142          36 :     for (auto & pos_group : _positions_2d)
     143          24 :       pos_group.resize(_extra_id_group_indices[1].size());
     144          36 :     for (auto & vol_group : volumes_2d)
     145          24 :       vol_group.resize(_extra_id_group_indices[1].size());
     146             :   }
     147          36 :   else if (_extra_id_names.size() == 3)
     148             :   {
     149          24 :     _positions_3d.resize(_extra_id_group_indices[0].size());
     150          24 :     volumes_3d.resize(_positions_3d.size());
     151          60 :     for (auto & vec_pos_group : _positions_3d)
     152             :     {
     153          36 :       vec_pos_group.resize(_extra_id_group_indices[1].size());
     154         108 :       for (auto & pos_group : vec_pos_group)
     155          72 :         pos_group.resize(_extra_id_group_indices[2].size());
     156             :     }
     157          60 :     for (auto & vec_vol_group : volumes_3d)
     158             :     {
     159          36 :       vec_vol_group.resize(_extra_id_group_indices[1].size());
     160         108 :       for (auto & vol_group : vec_vol_group)
     161          72 :         vol_group.resize(_extra_id_group_indices[2].size());
     162             :     }
     163             :   }
     164          12 :   else if (_extra_id_names.size() == 4)
     165             :   {
     166          12 :     _positions_4d.resize(_extra_id_group_indices[0].size());
     167          12 :     volumes_4d.resize(_positions_4d.size());
     168          36 :     for (auto & vec_vec_pos_group : _positions_4d)
     169             :     {
     170          24 :       vec_vec_pos_group.resize(_extra_id_group_indices[1].size());
     171          48 :       for (auto & vec_pos_group : vec_vec_pos_group)
     172             :       {
     173          24 :         vec_pos_group.resize(_extra_id_group_indices[2].size());
     174         120 :         for (auto & pos_group : vec_pos_group)
     175          96 :           pos_group.resize(_extra_id_group_indices[3].size());
     176             :       }
     177             :     }
     178          36 :     for (auto & vec_vec_vol_group : volumes_4d)
     179             :     {
     180          24 :       vec_vec_vol_group.resize(_extra_id_group_indices[1].size());
     181          48 :       for (auto & vec_vol_group : vec_vec_vol_group)
     182             :       {
     183          24 :         vec_vol_group.resize(_extra_id_group_indices[2].size());
     184         120 :         for (auto & vol_group : vec_vol_group)
     185          96 :           vol_group.resize(_extra_id_group_indices[3].size());
     186             :       }
     187             :     }
     188             :   }
     189             :   else
     190           0 :     mooseError("Too much dimensionality for positions");
     191             : 
     192             :   // To simplify retrieving the final positions vector
     193             :   auto getNestedPositions =
     194      134784 :       [this](const std::vector<unsigned int> & indices) -> std::vector<Point> &
     195             :   {
     196             :     mooseAssert(indices[_extra_id_names.size() - 1] == 0, "Indexing issue");
     197       57600 :     if (_extra_id_names.size() == 1)
     198       46080 :       return _positions;
     199       11520 :     else if (_extra_id_names.size() == 2)
     200        3456 :       return _positions_2d[indices[0]];
     201        8064 :     else if (_extra_id_names.size() == 3)
     202        5760 :       return _positions_3d[indices[0]][indices[1]];
     203             :     else
     204        2304 :       return _positions_4d[indices[0]][indices[1]][indices[2]];
     205          72 :   };
     206       57600 :   auto getNestedVolumes = [this, &volumes, &volumes_2d, &volumes_3d, &volumes_4d](
     207       88704 :                               const std::vector<unsigned int> & indices) -> std::vector<Real> &
     208             :   {
     209             :     mooseAssert(indices[_extra_id_names.size() - 1] == 0, "Indexing issue");
     210       57600 :     if (_extra_id_names.size() == 1)
     211       46080 :       return volumes;
     212       11520 :     else if (_extra_id_names.size() == 2)
     213        3456 :       return volumes_2d[indices[0]];
     214        8064 :     else if (_extra_id_names.size() == 3)
     215        5760 :       return volumes_3d[indices[0]][indices[1]];
     216             :     else
     217        2304 :       return volumes_4d[indices[0]][indices[1]][indices[2]];
     218          72 :   };
     219             : 
     220             :   std::vector<std::map<unsigned int, unsigned int>> positions_indexing(
     221          72 :       _extra_id_group_indices.size());
     222             : 
     223             :   // Make index maps to go from the extra element id to the positions array index
     224         240 :   for (const auto i : index_range(_extra_id_group_indices))
     225             :   {
     226         168 :     auto & indices = _extra_id_group_indices[i];
     227         168 :     auto j = 0;
     228         588 :     for (const auto extra_id : indices)
     229         420 :       positions_indexing[i][extra_id] = j++;
     230             :   }
     231             : 
     232      387144 :   for (const auto & elem : _mesh.getMesh().active_local_element_ptr_range())
     233             :   {
     234             :     // Pre-compute the centroid, this is expensive but may be used for multiple element ids
     235      193536 :     const auto centroid = elem->true_centroid();
     236      193536 :     const auto volume = elem->volume();
     237             : 
     238             :     // Keeps track of indices in multi-D arrays
     239      193536 :     std::vector<unsigned int> previous_indices(4);
     240             : 
     241      289728 :     for (const auto i : index_range(_extra_id_names))
     242             :     {
     243             :       auto iter =
     244      289728 :           positions_indexing[i].find(id(*elem, _extra_id_indices[i], _blocks_in_use && (i == 0)));
     245      289728 :       if (iter == positions_indexing[i].end())
     246      193536 :         break;
     247      153792 :       else if (_extra_id_names.size() == i + 1)
     248             :       {
     249       57600 :         getNestedPositions(previous_indices)[iter->second] += volume * centroid;
     250       57600 :         getNestedVolumes(previous_indices)[iter->second] += volume;
     251       57600 :         break;
     252             :       }
     253             :       else
     254       96192 :         previous_indices[i] = iter->second;
     255             :     }
     256      193608 :   }
     257             : 
     258             :   // Report the zero volumes
     259          72 :   unsigned int num_zeros = 0;
     260          72 :   if (_extra_id_names.size() == 1)
     261             :   {
     262          24 :     _mesh.comm().sum(volumes);
     263          24 :     _mesh.comm().sum(_positions);
     264          72 :     for (const auto & vol : volumes)
     265          48 :       if (MooseUtils::absoluteFuzzyEqual(vol, 0))
     266           0 :         num_zeros++;
     267             :   }
     268          72 :   if (_extra_id_names.size() == 2)
     269             :   {
     270          36 :     for (const auto i : index_range(volumes_2d))
     271             :     {
     272          24 :       _mesh.comm().sum(volumes_2d[i]);
     273          24 :       _mesh.comm().sum(_positions_2d[i]);
     274             :     }
     275          36 :     for (const auto & vol_vec : volumes_2d)
     276          96 :       for (const auto & vol : vol_vec)
     277          72 :         if (MooseUtils::absoluteFuzzyEqual(vol, 0))
     278           0 :           num_zeros++;
     279             :   }
     280          72 :   if (_extra_id_names.size() == 3)
     281             :   {
     282          60 :     for (const auto i : index_range(volumes_3d))
     283         108 :       for (const auto j : index_range(volumes_3d[i]))
     284             :       {
     285          72 :         _mesh.comm().sum(volumes_3d[i][j]);
     286          72 :         _mesh.comm().sum(_positions_3d[i][j]);
     287             :       }
     288          60 :     for (const auto & vol_vec_vec : volumes_3d)
     289         108 :       for (const auto & vol_vec : vol_vec_vec)
     290         336 :         for (const auto & vol : vol_vec)
     291         264 :           if (MooseUtils::absoluteFuzzyEqual(vol, 0))
     292          48 :             num_zeros++;
     293             :   }
     294          72 :   if (_extra_id_names.size() == 4)
     295             :   {
     296          36 :     for (const auto i : index_range(volumes_4d))
     297          48 :       for (const auto j : index_range(volumes_4d[i]))
     298         120 :         for (const auto k : index_range(volumes_4d[i][j]))
     299             :         {
     300          96 :           _mesh.comm().sum(volumes_4d[i][j][k]);
     301          96 :           _mesh.comm().sum(_positions_4d[i][j][k]);
     302             :         }
     303          36 :     for (const auto & vol_vec_vec_vec : volumes_4d)
     304          48 :       for (const auto & vol_vec_vec : vol_vec_vec_vec)
     305         120 :         for (const auto & vol_vec : vol_vec_vec)
     306         480 :           for (const auto & vol : vol_vec)
     307         384 :             if (MooseUtils::absoluteFuzzyEqual(vol, 0))
     308         288 :               num_zeros++;
     309             :   }
     310          72 :   if (num_zeros)
     311          24 :     mooseWarning(std::to_string(num_zeros) +
     312             :                  " zero volume bins detected during group centroid position calculation. "
     313             :                  "The corresponding positions will be removed from consideration.");
     314             : 
     315             :   // Renormalize by the total bin volumes
     316          72 :   if (_extra_id_names.size() == 1)
     317          72 :     for (MooseIndex(_positions) i = 0; i < _positions.size(); i++)
     318             :     {
     319          48 :       if (volumes[i] != 0)
     320          48 :         _positions[i] /= volumes[i];
     321             :       else
     322             :       {
     323           0 :         _positions.erase(_positions.begin() + i);
     324           0 :         volumes.erase(volumes.begin() + i);
     325           0 :         i--;
     326             :       }
     327             :     }
     328          48 :   else if (_extra_id_names.size() == 2)
     329          36 :     for (const auto i : index_range(_positions_2d))
     330          96 :       for (MooseIndex(_positions) j = 0; j < _positions_2d[i].size(); j++)
     331             :       {
     332          72 :         if (volumes_2d[i][j] != 0)
     333          72 :           _positions_2d[i][j] /= volumes_2d[i][j];
     334             :         else
     335             :         {
     336           0 :           _positions_2d[i].erase(_positions_2d[i].begin() + j);
     337           0 :           volumes_2d[i].erase(volumes_2d[i].begin() + j);
     338           0 :           j--;
     339             :         }
     340             :       }
     341          36 :   else if (_extra_id_names.size() == 3)
     342          60 :     for (const auto i : index_range(_positions_3d))
     343         108 :       for (const auto j : index_range(_positions_3d[i]))
     344         336 :         for (MooseIndex(_positions) k = 0; k < _positions_3d[i][j].size(); k++)
     345             :         {
     346         264 :           if (volumes_3d[i][j][k] != 0)
     347         216 :             _positions_3d[i][j][k] /= volumes_3d[i][j][k];
     348             :           else
     349             :           {
     350          48 :             _positions_3d[i][j].erase(_positions_3d[i][j].begin() + k);
     351          48 :             volumes_3d[i][j].erase(volumes_3d[i][j].begin() + k);
     352          48 :             k--;
     353             :           }
     354             :         }
     355          12 :   else if (_extra_id_names.size() == 4)
     356          36 :     for (const auto i : index_range(_positions_4d))
     357          48 :       for (const auto j : index_range(_positions_4d[i]))
     358         120 :         for (const auto k : index_range(_positions_4d[i][j]))
     359         480 :           for (MooseIndex(_positions) l = 0; l < _positions_4d[i][j][k].size(); l++)
     360             :           {
     361         384 :             if (volumes_4d[i][j][k][l] != 0)
     362          96 :               _positions_4d[i][j][k][l] /= volumes_4d[i][j][k][l];
     363             :             else
     364             :             {
     365         288 :               _positions_4d[i][j][k].erase(_positions_4d[i][j][k].begin() + l);
     366         288 :               volumes_4d[i][j][k].erase(volumes_4d[i][j][k].begin() + l);
     367         288 :               l--;
     368             :             }
     369             :           }
     370             : 
     371             :   // Fill the 1D position vector
     372          72 :   unrollMultiDPositions();
     373          72 :   _initialized = true;
     374          72 : }
     375             : 
     376             : unsigned int
     377      354240 : ElementGroupCentroidPositions::id(const Elem & elem, unsigned int id_index, bool use_subdomains)
     378             : {
     379             :   mooseAssert(!use_subdomains || (id_index == std::numeric_limits<unsigned short>::max()),
     380             :               "We do not expect a valid element extra integer index for subdomains");
     381      354240 :   if (use_subdomains)
     382      161280 :     return elem.subdomain_id();
     383             :   else
     384      192960 :     return elem.get_extra_integer(id_index);
     385             : }