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       14421 : ElementGroupCentroidPositions::validParams()
      18             : {
      19       14421 :   InputParameters params = Positions::validParams();
      20       14421 :   params += BlockRestrictable::validParams();
      21       14421 :   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       14421 :   params.addParam<MooseEnum>("grouping_type", groupTypeEnum(), "Type of group of elements");
      26       14421 :   params.addParam<std::vector<ExtraElementIDName>>("extra_id_name",
      27             :                                                    "Name(s) of the extra element ID(s) to use");
      28       14421 :   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       14421 :   params.set<bool>("auto_sort") = false;
      35             :   // Full replicated mesh is used on every rank to generate positions
      36       14421 :   params.set<bool>("auto_broadcast") = false;
      37             : 
      38       14421 :   return params;
      39           0 : }
      40             : 
      41          78 : ElementGroupCentroidPositions::ElementGroupCentroidPositions(const InputParameters & parameters)
      42             :   : Positions(parameters),
      43             :     BlockRestrictable(this),
      44         156 :     _mesh(_fe_problem.mesh()),
      45          78 :     _group_type(getParam<MooseEnum>("grouping_type"))
      46             : {
      47             :   // We are not excluding using both block restriction and extra element ids
      48          78 :   if (_group_type == "extra_id" || _group_type == "block_and_extra_id")
      49             :   {
      50          52 :     _extra_id_names = getParam<std::vector<ExtraElementIDName>>("extra_id_name");
      51         169 :     for (const auto & name : _extra_id_names)
      52         117 :       _extra_id_indices.push_back(_mesh.getMesh().get_elem_integer_index(name));
      53             : 
      54          52 :     if (isParamValid("extra_id"))
      55          52 :       _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          52 :     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          52 :     if (_extra_id_indices.size() > unsigned(3 + blockRestricted()))
      66           0 :       mooseError("Positions currently supports only up to 4D storage");
      67             :   }
      68             :   else
      69             :   {
      70          26 :     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          26 :     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          78 :   if (blockRestricted() || !isParamValid("extra_id_name"))
      80             :   {
      81          65 :     _blocks_in_use = true;
      82          65 :     _extra_id_names.insert(_extra_id_names.begin(), "block");
      83          65 :     _extra_id_indices.insert(_extra_id_indices.begin(), std::numeric_limits<unsigned short>::max());
      84          65 :     _extra_id_group_indices.insert(_extra_id_group_indices.begin(), std::vector<unsigned int>());
      85             :     // Add real block restriction
      86          65 :     if (blockRestricted())
      87         143 :       for (const auto & block : blockIDs())
      88          91 :         _extra_id_group_indices[0].push_back(block);
      89             :     // Just add all blocks
      90             :     else
      91          52 :       for (const auto & block : meshBlockIDs())
      92          39 :         _extra_id_group_indices[0].push_back(block);
      93             :   }
      94             :   else
      95          13 :     _blocks_in_use = false;
      96             : 
      97             :   // Mesh is ready at construction
      98          78 :   initialize();
      99             :   // Sort the output (by XYZ) if requested
     100          78 :   finalize();
     101          78 : }
     102             : 
     103             : void
     104          78 : ElementGroupCentroidPositions::initialize()
     105             : {
     106          78 :   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         260 :   for (const auto i : index_range(_extra_id_group_indices))
     112             :   {
     113         182 :     auto & indices = _extra_id_group_indices[i];
     114         182 :     if (indices.empty())
     115             :     {
     116          26 :       std::set<unsigned int> ids;
     117       71706 :       for (const auto & elem : _mesh.getMesh().active_local_element_ptr_range())
     118       71706 :         ids.insert(id(*elem, _extra_id_indices[i], _blocks_in_use && i == 0));
     119          26 :       _mesh.comm().set_union(ids);
     120         130 :       for (const auto & id : ids)
     121         104 :         indices.push_back(id);
     122          26 :     }
     123             :   }
     124             : 
     125             :   // Allocate some vectors holding the volumes
     126          78 :   std::vector<Real> volumes;
     127          78 :   std::vector<std::vector<Real>> volumes_2d;
     128          78 :   std::vector<std::vector<std::vector<Real>>> volumes_3d;
     129          78 :   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          78 :   if (_extra_id_names.size() == 1)
     134             :   {
     135          26 :     _positions.resize(_extra_id_group_indices[0].size());
     136          26 :     volumes.resize(_positions.size());
     137             :   }
     138          52 :   else if (_extra_id_names.size() == 2)
     139             :   {
     140          13 :     _positions_2d.resize(_extra_id_group_indices[0].size());
     141          13 :     volumes_2d.resize(_positions_2d.size());
     142          39 :     for (auto & pos_group : _positions_2d)
     143          26 :       pos_group.resize(_extra_id_group_indices[1].size());
     144          39 :     for (auto & vol_group : volumes_2d)
     145          26 :       vol_group.resize(_extra_id_group_indices[1].size());
     146             :   }
     147          39 :   else if (_extra_id_names.size() == 3)
     148             :   {
     149          26 :     _positions_3d.resize(_extra_id_group_indices[0].size());
     150          26 :     volumes_3d.resize(_positions_3d.size());
     151          65 :     for (auto & vec_pos_group : _positions_3d)
     152             :     {
     153          39 :       vec_pos_group.resize(_extra_id_group_indices[1].size());
     154         117 :       for (auto & pos_group : vec_pos_group)
     155          78 :         pos_group.resize(_extra_id_group_indices[2].size());
     156             :     }
     157          65 :     for (auto & vec_vol_group : volumes_3d)
     158             :     {
     159          39 :       vec_vol_group.resize(_extra_id_group_indices[1].size());
     160         117 :       for (auto & vol_group : vec_vol_group)
     161          78 :         vol_group.resize(_extra_id_group_indices[2].size());
     162             :     }
     163             :   }
     164          13 :   else if (_extra_id_names.size() == 4)
     165             :   {
     166          13 :     _positions_4d.resize(_extra_id_group_indices[0].size());
     167          13 :     volumes_4d.resize(_positions_4d.size());
     168          39 :     for (auto & vec_vec_pos_group : _positions_4d)
     169             :     {
     170          26 :       vec_vec_pos_group.resize(_extra_id_group_indices[1].size());
     171          52 :       for (auto & vec_pos_group : vec_vec_pos_group)
     172             :       {
     173          26 :         vec_pos_group.resize(_extra_id_group_indices[2].size());
     174         130 :         for (auto & pos_group : vec_pos_group)
     175         104 :           pos_group.resize(_extra_id_group_indices[3].size());
     176             :       }
     177             :     }
     178          39 :     for (auto & vec_vec_vol_group : volumes_4d)
     179             :     {
     180          26 :       vec_vec_vol_group.resize(_extra_id_group_indices[1].size());
     181          52 :       for (auto & vec_vol_group : vec_vec_vol_group)
     182             :       {
     183          26 :         vec_vol_group.resize(_extra_id_group_indices[2].size());
     184         130 :         for (auto & vol_group : vec_vol_group)
     185         104 :           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      149760 :       [this](const std::vector<unsigned int> & indices) -> std::vector<Point> &
     195             :   {
     196             :     mooseAssert(indices[_extra_id_names.size() - 1] == 0, "Indexing issue");
     197       64000 :     if (_extra_id_names.size() == 1)
     198       51200 :       return _positions;
     199       12800 :     else if (_extra_id_names.size() == 2)
     200        3840 :       return _positions_2d[indices[0]];
     201        8960 :     else if (_extra_id_names.size() == 3)
     202        6400 :       return _positions_3d[indices[0]][indices[1]];
     203             :     else
     204        2560 :       return _positions_4d[indices[0]][indices[1]][indices[2]];
     205          78 :   };
     206       64000 :   auto getNestedVolumes = [this, &volumes, &volumes_2d, &volumes_3d, &volumes_4d](
     207       98560 :                               const std::vector<unsigned int> & indices) -> std::vector<Real> &
     208             :   {
     209             :     mooseAssert(indices[_extra_id_names.size() - 1] == 0, "Indexing issue");
     210       64000 :     if (_extra_id_names.size() == 1)
     211       51200 :       return volumes;
     212       12800 :     else if (_extra_id_names.size() == 2)
     213        3840 :       return volumes_2d[indices[0]];
     214        8960 :     else if (_extra_id_names.size() == 3)
     215        6400 :       return volumes_3d[indices[0]][indices[1]];
     216             :     else
     217        2560 :       return volumes_4d[indices[0]][indices[1]][indices[2]];
     218          78 :   };
     219             : 
     220             :   std::vector<std::map<unsigned int, unsigned int>> positions_indexing(
     221          78 :       _extra_id_group_indices.size());
     222             : 
     223             :   // Make index maps to go from the extra element id to the positions array index
     224         260 :   for (const auto i : index_range(_extra_id_group_indices))
     225             :   {
     226         182 :     auto & indices = _extra_id_group_indices[i];
     227         182 :     auto j = 0;
     228         637 :     for (const auto extra_id : indices)
     229         455 :       positions_indexing[i][extra_id] = j++;
     230             :   }
     231             : 
     232      430158 :   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      215040 :     const auto centroid = elem->true_centroid();
     236      215040 :     const auto volume = elem->volume();
     237             : 
     238             :     // Keeps track of indices in multi-D arrays
     239      215040 :     std::vector<unsigned int> previous_indices(4);
     240             : 
     241      321920 :     for (const auto i : index_range(_extra_id_names))
     242             :     {
     243             :       auto iter =
     244      321920 :           positions_indexing[i].find(id(*elem, _extra_id_indices[i], _blocks_in_use && (i == 0)));
     245      321920 :       if (iter == positions_indexing[i].end())
     246      215040 :         break;
     247      170880 :       else if (_extra_id_names.size() == i + 1)
     248             :       {
     249       64000 :         getNestedPositions(previous_indices)[iter->second] += volume * centroid;
     250       64000 :         getNestedVolumes(previous_indices)[iter->second] += volume;
     251       64000 :         break;
     252             :       }
     253             :       else
     254      106880 :         previous_indices[i] = iter->second;
     255             :     }
     256      215118 :   }
     257             : 
     258             :   // Report the zero volumes
     259          78 :   unsigned int num_zeros = 0;
     260          78 :   if (_extra_id_names.size() == 1)
     261             :   {
     262          26 :     _mesh.comm().sum(volumes);
     263          26 :     _mesh.comm().sum(_positions);
     264          78 :     for (const auto & vol : volumes)
     265          52 :       if (MooseUtils::absoluteFuzzyEqual(vol, 0))
     266           0 :         num_zeros++;
     267             :   }
     268          78 :   if (_extra_id_names.size() == 2)
     269             :   {
     270          39 :     for (const auto i : index_range(volumes_2d))
     271             :     {
     272          26 :       _mesh.comm().sum(volumes_2d[i]);
     273          26 :       _mesh.comm().sum(_positions_2d[i]);
     274             :     }
     275          39 :     for (const auto & vol_vec : volumes_2d)
     276         104 :       for (const auto & vol : vol_vec)
     277          78 :         if (MooseUtils::absoluteFuzzyEqual(vol, 0))
     278           0 :           num_zeros++;
     279             :   }
     280          78 :   if (_extra_id_names.size() == 3)
     281             :   {
     282          65 :     for (const auto i : index_range(volumes_3d))
     283         117 :       for (const auto j : index_range(volumes_3d[i]))
     284             :       {
     285          78 :         _mesh.comm().sum(volumes_3d[i][j]);
     286          78 :         _mesh.comm().sum(_positions_3d[i][j]);
     287             :       }
     288          65 :     for (const auto & vol_vec_vec : volumes_3d)
     289         117 :       for (const auto & vol_vec : vol_vec_vec)
     290         364 :         for (const auto & vol : vol_vec)
     291         286 :           if (MooseUtils::absoluteFuzzyEqual(vol, 0))
     292          52 :             num_zeros++;
     293             :   }
     294          78 :   if (_extra_id_names.size() == 4)
     295             :   {
     296          39 :     for (const auto i : index_range(volumes_4d))
     297          52 :       for (const auto j : index_range(volumes_4d[i]))
     298         130 :         for (const auto k : index_range(volumes_4d[i][j]))
     299             :         {
     300         104 :           _mesh.comm().sum(volumes_4d[i][j][k]);
     301         104 :           _mesh.comm().sum(_positions_4d[i][j][k]);
     302             :         }
     303          39 :     for (const auto & vol_vec_vec_vec : volumes_4d)
     304          52 :       for (const auto & vol_vec_vec : vol_vec_vec_vec)
     305         130 :         for (const auto & vol_vec : vol_vec_vec)
     306         520 :           for (const auto & vol : vol_vec)
     307         416 :             if (MooseUtils::absoluteFuzzyEqual(vol, 0))
     308         312 :               num_zeros++;
     309             :   }
     310          78 :   if (num_zeros)
     311          26 :     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          78 :   if (_extra_id_names.size() == 1)
     317          78 :     for (MooseIndex(_positions) i = 0; i < _positions.size(); i++)
     318             :     {
     319          52 :       if (volumes[i] != 0)
     320          52 :         _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          52 :   else if (_extra_id_names.size() == 2)
     329          39 :     for (const auto i : index_range(_positions_2d))
     330         104 :       for (MooseIndex(_positions) j = 0; j < _positions_2d[i].size(); j++)
     331             :       {
     332          78 :         if (volumes_2d[i][j] != 0)
     333          78 :           _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          39 :   else if (_extra_id_names.size() == 3)
     342          65 :     for (const auto i : index_range(_positions_3d))
     343         117 :       for (const auto j : index_range(_positions_3d[i]))
     344         364 :         for (MooseIndex(_positions) k = 0; k < _positions_3d[i][j].size(); k++)
     345             :         {
     346         286 :           if (volumes_3d[i][j][k] != 0)
     347         234 :             _positions_3d[i][j][k] /= volumes_3d[i][j][k];
     348             :           else
     349             :           {
     350          52 :             _positions_3d[i][j].erase(_positions_3d[i][j].begin() + k);
     351          52 :             volumes_3d[i][j].erase(volumes_3d[i][j].begin() + k);
     352          52 :             k--;
     353             :           }
     354             :         }
     355          13 :   else if (_extra_id_names.size() == 4)
     356          39 :     for (const auto i : index_range(_positions_4d))
     357          52 :       for (const auto j : index_range(_positions_4d[i]))
     358         130 :         for (const auto k : index_range(_positions_4d[i][j]))
     359         520 :           for (MooseIndex(_positions) l = 0; l < _positions_4d[i][j][k].size(); l++)
     360             :           {
     361         416 :             if (volumes_4d[i][j][k][l] != 0)
     362         104 :               _positions_4d[i][j][k][l] /= volumes_4d[i][j][k][l];
     363             :             else
     364             :             {
     365         312 :               _positions_4d[i][j][k].erase(_positions_4d[i][j][k].begin() + l);
     366         312 :               volumes_4d[i][j][k].erase(volumes_4d[i][j][k].begin() + l);
     367         312 :               l--;
     368             :             }
     369             :           }
     370             : 
     371             :   // Fill the 1D position vector
     372          78 :   unrollMultiDPositions();
     373          78 :   _initialized = true;
     374          78 : }
     375             : 
     376             : unsigned int
     377      393600 : 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      393600 :   if (use_subdomains)
     382      179200 :     return elem.subdomain_id();
     383             :   else
     384      214400 :     return elem.get_extra_integer(id_index);
     385             : }