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        3205 : ElementGroupCentroidPositions::validParams()
      18             : {
      19        3205 :   InputParameters params = Positions::validParams();
      20        3205 :   params += BlockRestrictable::validParams();
      21        6410 :   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       12820 :   params.addParam<MooseEnum>("grouping_type", groupTypeEnum(), "Type of group of elements");
      26       12820 :   params.addParam<std::vector<ExtraElementIDName>>("extra_id_name",
      27             :                                                    "Name(s) of the extra element ID(s) to use");
      28       12820 :   params.addParam<std::vector<std::vector<dof_id_type>>>(
      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        6410 :   params.set<bool>("auto_sort") = false;
      35             :   // Full replicated mesh is used on every rank to generate positions
      36        3205 :   params.set<bool>("auto_broadcast") = false;
      37             : 
      38        3205 :   return params;
      39           0 : }
      40             : 
      41          72 : ElementGroupCentroidPositions::ElementGroupCentroidPositions(const InputParameters & parameters)
      42             :   : Positions(parameters),
      43             :     BlockRestrictable(this),
      44         144 :     _mesh(_subproblem.mesh()),
      45         216 :     _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          96 :     _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         144 :     if (isParamValid("extra_id"))
      55         144 :       _extra_id_group_indices = getParam<std::vector<std::vector<dof_id_type>>>("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          72 :     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          72 :     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         120 :   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<dof_id_type>());
      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<dof_id_type> ids;
     117       64536 :       for (const auto & elem : _mesh.getMesh().active_local_element_ptr_range())
     118             :       {
     119       64512 :         auto eeid = id(*elem, _extra_id_indices[i], _blocks_in_use && i == 0);
     120       64512 :         if (eeid != DofObject::invalid_id)
     121       64512 :           ids.insert(eeid);
     122          24 :       }
     123          24 :       _mesh.comm().set_union(ids);
     124         120 :       for (const auto & id : ids)
     125          96 :         indices.push_back(id);
     126          24 :     }
     127             :   }
     128             : 
     129             :   // Allocate some vectors holding the volumes
     130          72 :   std::vector<Real> volumes;
     131          72 :   std::vector<std::vector<Real>> volumes_2d;
     132          72 :   std::vector<std::vector<std::vector<Real>>> volumes_3d;
     133          72 :   std::vector<std::vector<std::vector<std::vector<Real>>>> volumes_4d;
     134             : 
     135             :   // Default indexing: blocks, extra_id_name_1, extra_id_name_2 ...
     136             :   // Allocate vectors holding the positions
     137          72 :   if (_extra_id_names.size() == 1)
     138             :   {
     139          24 :     _positions.resize(_extra_id_group_indices[0].size());
     140          24 :     volumes.resize(_positions.size());
     141             :   }
     142          48 :   else if (_extra_id_names.size() == 2)
     143             :   {
     144          12 :     _positions_2d.resize(_extra_id_group_indices[0].size());
     145          12 :     volumes_2d.resize(_positions_2d.size());
     146          36 :     for (auto & pos_group : _positions_2d)
     147          24 :       pos_group.resize(_extra_id_group_indices[1].size());
     148          36 :     for (auto & vol_group : volumes_2d)
     149          24 :       vol_group.resize(_extra_id_group_indices[1].size());
     150             :   }
     151          36 :   else if (_extra_id_names.size() == 3)
     152             :   {
     153          24 :     _positions_3d.resize(_extra_id_group_indices[0].size());
     154          24 :     volumes_3d.resize(_positions_3d.size());
     155          60 :     for (auto & vec_pos_group : _positions_3d)
     156             :     {
     157          36 :       vec_pos_group.resize(_extra_id_group_indices[1].size());
     158         108 :       for (auto & pos_group : vec_pos_group)
     159          72 :         pos_group.resize(_extra_id_group_indices[2].size());
     160             :     }
     161          60 :     for (auto & vec_vol_group : volumes_3d)
     162             :     {
     163          36 :       vec_vol_group.resize(_extra_id_group_indices[1].size());
     164         108 :       for (auto & vol_group : vec_vol_group)
     165          72 :         vol_group.resize(_extra_id_group_indices[2].size());
     166             :     }
     167             :   }
     168          12 :   else if (_extra_id_names.size() == 4)
     169             :   {
     170          12 :     _positions_4d.resize(_extra_id_group_indices[0].size());
     171          12 :     volumes_4d.resize(_positions_4d.size());
     172          36 :     for (auto & vec_vec_pos_group : _positions_4d)
     173             :     {
     174          24 :       vec_vec_pos_group.resize(_extra_id_group_indices[1].size());
     175          48 :       for (auto & vec_pos_group : vec_vec_pos_group)
     176             :       {
     177          24 :         vec_pos_group.resize(_extra_id_group_indices[2].size());
     178         120 :         for (auto & pos_group : vec_pos_group)
     179          96 :           pos_group.resize(_extra_id_group_indices[3].size());
     180             :       }
     181             :     }
     182          36 :     for (auto & vec_vec_vol_group : volumes_4d)
     183             :     {
     184          24 :       vec_vec_vol_group.resize(_extra_id_group_indices[1].size());
     185          48 :       for (auto & vec_vol_group : vec_vec_vol_group)
     186             :       {
     187          24 :         vec_vol_group.resize(_extra_id_group_indices[2].size());
     188         120 :         for (auto & vol_group : vec_vol_group)
     189          96 :           vol_group.resize(_extra_id_group_indices[3].size());
     190             :       }
     191             :     }
     192             :   }
     193             :   else
     194           0 :     mooseError("Too much dimensionality for positions");
     195             : 
     196             :   // To simplify retrieving the final positions vector
     197             :   auto getNestedPositions =
     198       57600 :       [this](const std::vector<unsigned int> & indices) -> std::vector<Point> &
     199             :   {
     200             :     mooseAssert(indices[_extra_id_names.size() - 1] == 0, "Indexing issue");
     201       57600 :     if (_extra_id_names.size() == 1)
     202       46080 :       return _positions;
     203       11520 :     else if (_extra_id_names.size() == 2)
     204        3456 :       return _positions_2d[indices[0]];
     205        8064 :     else if (_extra_id_names.size() == 3)
     206        5760 :       return _positions_3d[indices[0]][indices[1]];
     207             :     else
     208        2304 :       return _positions_4d[indices[0]][indices[1]][indices[2]];
     209          72 :   };
     210       57600 :   auto getNestedVolumes = [this, &volumes, &volumes_2d, &volumes_3d, &volumes_4d](
     211             :                               const std::vector<unsigned int> & indices) -> std::vector<Real> &
     212             :   {
     213             :     mooseAssert(indices[_extra_id_names.size() - 1] == 0, "Indexing issue");
     214       57600 :     if (_extra_id_names.size() == 1)
     215       46080 :       return volumes;
     216       11520 :     else if (_extra_id_names.size() == 2)
     217        3456 :       return volumes_2d[indices[0]];
     218        8064 :     else if (_extra_id_names.size() == 3)
     219        5760 :       return volumes_3d[indices[0]][indices[1]];
     220             :     else
     221        2304 :       return volumes_4d[indices[0]][indices[1]][indices[2]];
     222          72 :   };
     223             : 
     224             :   std::vector<std::map<unsigned int, unsigned int>> positions_indexing(
     225          72 :       _extra_id_group_indices.size());
     226             : 
     227             :   // Make index maps to go from the extra element id to the positions array index
     228         240 :   for (const auto i : index_range(_extra_id_group_indices))
     229             :   {
     230         168 :     auto & indices = _extra_id_group_indices[i];
     231         168 :     auto j = 0;
     232         588 :     for (const auto extra_id : indices)
     233         420 :       positions_indexing[i][extra_id] = j++;
     234             :   }
     235             : 
     236      193608 :   for (const auto & elem : _mesh.getMesh().active_local_element_ptr_range())
     237             :   {
     238             :     // Pre-compute the centroid, this is expensive but may be used for multiple element ids
     239      193536 :     const auto centroid = elem->true_centroid();
     240      193536 :     const auto volume = elem->volume();
     241             : 
     242             :     // Keeps track of indices in multi-D arrays
     243      193536 :     std::vector<unsigned int> previous_indices(4);
     244             : 
     245      289728 :     for (const auto i : index_range(_extra_id_names))
     246             :     {
     247             :       auto iter =
     248      289728 :           positions_indexing[i].find(id(*elem, _extra_id_indices[i], _blocks_in_use && (i == 0)));
     249      289728 :       if (iter == positions_indexing[i].end())
     250      193536 :         break;
     251      153792 :       else if (_extra_id_names.size() == i + 1)
     252             :       {
     253       57600 :         getNestedPositions(previous_indices)[iter->second] += volume * centroid;
     254       57600 :         getNestedVolumes(previous_indices)[iter->second] += volume;
     255       57600 :         break;
     256             :       }
     257             :       else
     258       96192 :         previous_indices[i] = iter->second;
     259             :     }
     260      193608 :   }
     261             : 
     262             :   // Report the zero volumes
     263          72 :   unsigned int num_zeros = 0;
     264          72 :   if (_extra_id_names.size() == 1)
     265             :   {
     266          24 :     _mesh.comm().sum(volumes);
     267          24 :     _mesh.comm().sum(_positions);
     268          72 :     for (const auto & vol : volumes)
     269          48 :       if (MooseUtils::absoluteFuzzyEqual(vol, 0))
     270           0 :         num_zeros++;
     271             :   }
     272          72 :   if (_extra_id_names.size() == 2)
     273             :   {
     274          36 :     for (const auto i : index_range(volumes_2d))
     275             :     {
     276          24 :       _mesh.comm().sum(volumes_2d[i]);
     277          24 :       _mesh.comm().sum(_positions_2d[i]);
     278             :     }
     279          36 :     for (const auto & vol_vec : volumes_2d)
     280          96 :       for (const auto & vol : vol_vec)
     281          72 :         if (MooseUtils::absoluteFuzzyEqual(vol, 0))
     282           0 :           num_zeros++;
     283             :   }
     284          72 :   if (_extra_id_names.size() == 3)
     285             :   {
     286          60 :     for (const auto i : index_range(volumes_3d))
     287         108 :       for (const auto j : index_range(volumes_3d[i]))
     288             :       {
     289          72 :         _mesh.comm().sum(volumes_3d[i][j]);
     290          72 :         _mesh.comm().sum(_positions_3d[i][j]);
     291             :       }
     292          60 :     for (const auto & vol_vec_vec : volumes_3d)
     293         108 :       for (const auto & vol_vec : vol_vec_vec)
     294         336 :         for (const auto & vol : vol_vec)
     295         264 :           if (MooseUtils::absoluteFuzzyEqual(vol, 0))
     296          48 :             num_zeros++;
     297             :   }
     298          72 :   if (_extra_id_names.size() == 4)
     299             :   {
     300          36 :     for (const auto i : index_range(volumes_4d))
     301          48 :       for (const auto j : index_range(volumes_4d[i]))
     302         120 :         for (const auto k : index_range(volumes_4d[i][j]))
     303             :         {
     304          96 :           _mesh.comm().sum(volumes_4d[i][j][k]);
     305          96 :           _mesh.comm().sum(_positions_4d[i][j][k]);
     306             :         }
     307          36 :     for (const auto & vol_vec_vec_vec : volumes_4d)
     308          48 :       for (const auto & vol_vec_vec : vol_vec_vec_vec)
     309         120 :         for (const auto & vol_vec : vol_vec_vec)
     310         480 :           for (const auto & vol : vol_vec)
     311         384 :             if (MooseUtils::absoluteFuzzyEqual(vol, 0))
     312         288 :               num_zeros++;
     313             :   }
     314          72 :   if (num_zeros)
     315          24 :     mooseWarning(std::to_string(num_zeros) +
     316             :                  " zero volume bins detected during group centroid position calculation. "
     317             :                  "The corresponding positions will be removed from consideration.");
     318             : 
     319             :   // Renormalize by the total bin volumes
     320          72 :   if (_extra_id_names.size() == 1)
     321          72 :     for (MooseIndex(_positions) i = 0; i < _positions.size(); i++)
     322             :     {
     323          48 :       if (volumes[i] != 0)
     324          48 :         _positions[i] /= volumes[i];
     325             :       else
     326             :       {
     327           0 :         _positions.erase(_positions.begin() + i);
     328           0 :         volumes.erase(volumes.begin() + i);
     329           0 :         i--;
     330             :       }
     331             :     }
     332          48 :   else if (_extra_id_names.size() == 2)
     333          36 :     for (const auto i : index_range(_positions_2d))
     334          96 :       for (MooseIndex(_positions) j = 0; j < _positions_2d[i].size(); j++)
     335             :       {
     336          72 :         if (volumes_2d[i][j] != 0)
     337          72 :           _positions_2d[i][j] /= volumes_2d[i][j];
     338             :         else
     339             :         {
     340           0 :           _positions_2d[i].erase(_positions_2d[i].begin() + j);
     341           0 :           volumes_2d[i].erase(volumes_2d[i].begin() + j);
     342           0 :           j--;
     343             :         }
     344             :       }
     345          36 :   else if (_extra_id_names.size() == 3)
     346          60 :     for (const auto i : index_range(_positions_3d))
     347         108 :       for (const auto j : index_range(_positions_3d[i]))
     348         336 :         for (MooseIndex(_positions) k = 0; k < _positions_3d[i][j].size(); k++)
     349             :         {
     350         264 :           if (volumes_3d[i][j][k] != 0)
     351         216 :             _positions_3d[i][j][k] /= volumes_3d[i][j][k];
     352             :           else
     353             :           {
     354          48 :             _positions_3d[i][j].erase(_positions_3d[i][j].begin() + k);
     355          48 :             volumes_3d[i][j].erase(volumes_3d[i][j].begin() + k);
     356          48 :             k--;
     357             :           }
     358             :         }
     359          12 :   else if (_extra_id_names.size() == 4)
     360          36 :     for (const auto i : index_range(_positions_4d))
     361          48 :       for (const auto j : index_range(_positions_4d[i]))
     362         120 :         for (const auto k : index_range(_positions_4d[i][j]))
     363         480 :           for (MooseIndex(_positions) l = 0; l < _positions_4d[i][j][k].size(); l++)
     364             :           {
     365         384 :             if (volumes_4d[i][j][k][l] != 0)
     366          96 :               _positions_4d[i][j][k][l] /= volumes_4d[i][j][k][l];
     367             :             else
     368             :             {
     369         288 :               _positions_4d[i][j][k].erase(_positions_4d[i][j][k].begin() + l);
     370         288 :               volumes_4d[i][j][k].erase(volumes_4d[i][j][k].begin() + l);
     371         288 :               l--;
     372             :             }
     373             :           }
     374             : 
     375             :   // Fill the 1D position vector
     376          72 :   unrollMultiDPositions();
     377          72 :   _initialized = true;
     378          72 : }
     379             : 
     380             : dof_id_type
     381      354240 : ElementGroupCentroidPositions::id(const Elem & elem, unsigned int id_index, bool use_subdomains)
     382             : {
     383             :   mooseAssert(!use_subdomains || (id_index == std::numeric_limits<unsigned short>::max()),
     384             :               "We do not expect a valid element extra integer index for subdomains");
     385      354240 :   if (use_subdomains)
     386      161280 :     return elem.subdomain_id();
     387             :   else
     388      192960 :     return elem.get_extra_integer(id_index);
     389             : }