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       15023 : ElementGroupCentroidPositions::validParams()
      18             : {
      19       15023 :   InputParameters params = Positions::validParams();
      20       15023 :   params += BlockRestrictable::validParams();
      21       30046 :   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       60092 :   params.addParam<MooseEnum>("grouping_type", groupTypeEnum(), "Type of group of elements");
      26       60092 :   params.addParam<std::vector<ExtraElementIDName>>("extra_id_name",
      27             :                                                    "Name(s) of the extra element ID(s) to use");
      28       60092 :   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       30046 :   params.set<bool>("auto_sort") = false;
      35             :   // Full replicated mesh is used on every rank to generate positions
      36       15023 :   params.set<bool>("auto_broadcast") = false;
      37             : 
      38       15023 :   return params;
      39           0 : }
      40             : 
      41          78 : ElementGroupCentroidPositions::ElementGroupCentroidPositions(const InputParameters & parameters)
      42             :   : Positions(parameters),
      43             :     BlockRestrictable(this),
      44         156 :     _mesh(_subproblem.mesh()),
      45         234 :     _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         104 :     _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         156 :     if (isParamValid("extra_id"))
      55         156 :       _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          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          78 :     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          78 :     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         130 :   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<dof_id_type>());
      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<dof_id_type> ids;
     117       71706 :       for (const auto & elem : _mesh.getMesh().active_local_element_ptr_range())
     118             :       {
     119       71680 :         auto eeid = id(*elem, _extra_id_indices[i], _blocks_in_use && i == 0);
     120       71680 :         if (eeid != DofObject::invalid_id)
     121       71680 :           ids.insert(eeid);
     122          26 :       }
     123          26 :       _mesh.comm().set_union(ids);
     124         130 :       for (const auto & id : ids)
     125         104 :         indices.push_back(id);
     126          26 :     }
     127             :   }
     128             : 
     129             :   // Allocate some vectors holding the volumes
     130          78 :   std::vector<Real> volumes;
     131          78 :   std::vector<std::vector<Real>> volumes_2d;
     132          78 :   std::vector<std::vector<std::vector<Real>>> volumes_3d;
     133          78 :   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          78 :   if (_extra_id_names.size() == 1)
     138             :   {
     139          26 :     _positions.resize(_extra_id_group_indices[0].size());
     140          26 :     volumes.resize(_positions.size());
     141             :   }
     142          52 :   else if (_extra_id_names.size() == 2)
     143             :   {
     144          13 :     _positions_2d.resize(_extra_id_group_indices[0].size());
     145          13 :     volumes_2d.resize(_positions_2d.size());
     146          39 :     for (auto & pos_group : _positions_2d)
     147          26 :       pos_group.resize(_extra_id_group_indices[1].size());
     148          39 :     for (auto & vol_group : volumes_2d)
     149          26 :       vol_group.resize(_extra_id_group_indices[1].size());
     150             :   }
     151          39 :   else if (_extra_id_names.size() == 3)
     152             :   {
     153          26 :     _positions_3d.resize(_extra_id_group_indices[0].size());
     154          26 :     volumes_3d.resize(_positions_3d.size());
     155          65 :     for (auto & vec_pos_group : _positions_3d)
     156             :     {
     157          39 :       vec_pos_group.resize(_extra_id_group_indices[1].size());
     158         117 :       for (auto & pos_group : vec_pos_group)
     159          78 :         pos_group.resize(_extra_id_group_indices[2].size());
     160             :     }
     161          65 :     for (auto & vec_vol_group : volumes_3d)
     162             :     {
     163          39 :       vec_vol_group.resize(_extra_id_group_indices[1].size());
     164         117 :       for (auto & vol_group : vec_vol_group)
     165          78 :         vol_group.resize(_extra_id_group_indices[2].size());
     166             :     }
     167             :   }
     168          13 :   else if (_extra_id_names.size() == 4)
     169             :   {
     170          13 :     _positions_4d.resize(_extra_id_group_indices[0].size());
     171          13 :     volumes_4d.resize(_positions_4d.size());
     172          39 :     for (auto & vec_vec_pos_group : _positions_4d)
     173             :     {
     174          26 :       vec_vec_pos_group.resize(_extra_id_group_indices[1].size());
     175          52 :       for (auto & vec_pos_group : vec_vec_pos_group)
     176             :       {
     177          26 :         vec_pos_group.resize(_extra_id_group_indices[2].size());
     178         130 :         for (auto & pos_group : vec_pos_group)
     179         104 :           pos_group.resize(_extra_id_group_indices[3].size());
     180             :       }
     181             :     }
     182          39 :     for (auto & vec_vec_vol_group : volumes_4d)
     183             :     {
     184          26 :       vec_vec_vol_group.resize(_extra_id_group_indices[1].size());
     185          52 :       for (auto & vec_vol_group : vec_vec_vol_group)
     186             :       {
     187          26 :         vec_vol_group.resize(_extra_id_group_indices[2].size());
     188         130 :         for (auto & vol_group : vec_vol_group)
     189         104 :           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       64000 :       [this](const std::vector<unsigned int> & indices) -> std::vector<Point> &
     199             :   {
     200             :     mooseAssert(indices[_extra_id_names.size() - 1] == 0, "Indexing issue");
     201       64000 :     if (_extra_id_names.size() == 1)
     202       51200 :       return _positions;
     203       12800 :     else if (_extra_id_names.size() == 2)
     204        3840 :       return _positions_2d[indices[0]];
     205        8960 :     else if (_extra_id_names.size() == 3)
     206        6400 :       return _positions_3d[indices[0]][indices[1]];
     207             :     else
     208        2560 :       return _positions_4d[indices[0]][indices[1]][indices[2]];
     209          78 :   };
     210       64000 :   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       64000 :     if (_extra_id_names.size() == 1)
     215       51200 :       return volumes;
     216       12800 :     else if (_extra_id_names.size() == 2)
     217        3840 :       return volumes_2d[indices[0]];
     218        8960 :     else if (_extra_id_names.size() == 3)
     219        6400 :       return volumes_3d[indices[0]][indices[1]];
     220             :     else
     221        2560 :       return volumes_4d[indices[0]][indices[1]][indices[2]];
     222          78 :   };
     223             : 
     224             :   std::vector<std::map<unsigned int, unsigned int>> positions_indexing(
     225          78 :       _extra_id_group_indices.size());
     226             : 
     227             :   // Make index maps to go from the extra element id to the positions array index
     228         260 :   for (const auto i : index_range(_extra_id_group_indices))
     229             :   {
     230         182 :     auto & indices = _extra_id_group_indices[i];
     231         182 :     auto j = 0;
     232         637 :     for (const auto extra_id : indices)
     233         455 :       positions_indexing[i][extra_id] = j++;
     234             :   }
     235             : 
     236      215118 :   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      215040 :     const auto centroid = elem->true_centroid();
     240      215040 :     const auto volume = elem->volume();
     241             : 
     242             :     // Keeps track of indices in multi-D arrays
     243      215040 :     std::vector<unsigned int> previous_indices(4);
     244             : 
     245      321920 :     for (const auto i : index_range(_extra_id_names))
     246             :     {
     247             :       auto iter =
     248      321920 :           positions_indexing[i].find(id(*elem, _extra_id_indices[i], _blocks_in_use && (i == 0)));
     249      321920 :       if (iter == positions_indexing[i].end())
     250      215040 :         break;
     251      170880 :       else if (_extra_id_names.size() == i + 1)
     252             :       {
     253       64000 :         getNestedPositions(previous_indices)[iter->second] += volume * centroid;
     254       64000 :         getNestedVolumes(previous_indices)[iter->second] += volume;
     255       64000 :         break;
     256             :       }
     257             :       else
     258      106880 :         previous_indices[i] = iter->second;
     259             :     }
     260      215118 :   }
     261             : 
     262             :   // Report the zero volumes
     263          78 :   unsigned int num_zeros = 0;
     264          78 :   if (_extra_id_names.size() == 1)
     265             :   {
     266          26 :     _mesh.comm().sum(volumes);
     267          26 :     _mesh.comm().sum(_positions);
     268          78 :     for (const auto & vol : volumes)
     269          52 :       if (MooseUtils::absoluteFuzzyEqual(vol, 0))
     270           0 :         num_zeros++;
     271             :   }
     272          78 :   if (_extra_id_names.size() == 2)
     273             :   {
     274          39 :     for (const auto i : index_range(volumes_2d))
     275             :     {
     276          26 :       _mesh.comm().sum(volumes_2d[i]);
     277          26 :       _mesh.comm().sum(_positions_2d[i]);
     278             :     }
     279          39 :     for (const auto & vol_vec : volumes_2d)
     280         104 :       for (const auto & vol : vol_vec)
     281          78 :         if (MooseUtils::absoluteFuzzyEqual(vol, 0))
     282           0 :           num_zeros++;
     283             :   }
     284          78 :   if (_extra_id_names.size() == 3)
     285             :   {
     286          65 :     for (const auto i : index_range(volumes_3d))
     287         117 :       for (const auto j : index_range(volumes_3d[i]))
     288             :       {
     289          78 :         _mesh.comm().sum(volumes_3d[i][j]);
     290          78 :         _mesh.comm().sum(_positions_3d[i][j]);
     291             :       }
     292          65 :     for (const auto & vol_vec_vec : volumes_3d)
     293         117 :       for (const auto & vol_vec : vol_vec_vec)
     294         364 :         for (const auto & vol : vol_vec)
     295         286 :           if (MooseUtils::absoluteFuzzyEqual(vol, 0))
     296          52 :             num_zeros++;
     297             :   }
     298          78 :   if (_extra_id_names.size() == 4)
     299             :   {
     300          39 :     for (const auto i : index_range(volumes_4d))
     301          52 :       for (const auto j : index_range(volumes_4d[i]))
     302         130 :         for (const auto k : index_range(volumes_4d[i][j]))
     303             :         {
     304         104 :           _mesh.comm().sum(volumes_4d[i][j][k]);
     305         104 :           _mesh.comm().sum(_positions_4d[i][j][k]);
     306             :         }
     307          39 :     for (const auto & vol_vec_vec_vec : volumes_4d)
     308          52 :       for (const auto & vol_vec_vec : vol_vec_vec_vec)
     309         130 :         for (const auto & vol_vec : vol_vec_vec)
     310         520 :           for (const auto & vol : vol_vec)
     311         416 :             if (MooseUtils::absoluteFuzzyEqual(vol, 0))
     312         312 :               num_zeros++;
     313             :   }
     314          78 :   if (num_zeros)
     315          26 :     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          78 :   if (_extra_id_names.size() == 1)
     321          78 :     for (MooseIndex(_positions) i = 0; i < _positions.size(); i++)
     322             :     {
     323          52 :       if (volumes[i] != 0)
     324          52 :         _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          52 :   else if (_extra_id_names.size() == 2)
     333          39 :     for (const auto i : index_range(_positions_2d))
     334         104 :       for (MooseIndex(_positions) j = 0; j < _positions_2d[i].size(); j++)
     335             :       {
     336          78 :         if (volumes_2d[i][j] != 0)
     337          78 :           _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          39 :   else if (_extra_id_names.size() == 3)
     346          65 :     for (const auto i : index_range(_positions_3d))
     347         117 :       for (const auto j : index_range(_positions_3d[i]))
     348         364 :         for (MooseIndex(_positions) k = 0; k < _positions_3d[i][j].size(); k++)
     349             :         {
     350         286 :           if (volumes_3d[i][j][k] != 0)
     351         234 :             _positions_3d[i][j][k] /= volumes_3d[i][j][k];
     352             :           else
     353             :           {
     354          52 :             _positions_3d[i][j].erase(_positions_3d[i][j].begin() + k);
     355          52 :             volumes_3d[i][j].erase(volumes_3d[i][j].begin() + k);
     356          52 :             k--;
     357             :           }
     358             :         }
     359          13 :   else if (_extra_id_names.size() == 4)
     360          39 :     for (const auto i : index_range(_positions_4d))
     361          52 :       for (const auto j : index_range(_positions_4d[i]))
     362         130 :         for (const auto k : index_range(_positions_4d[i][j]))
     363         520 :           for (MooseIndex(_positions) l = 0; l < _positions_4d[i][j][k].size(); l++)
     364             :           {
     365         416 :             if (volumes_4d[i][j][k][l] != 0)
     366         104 :               _positions_4d[i][j][k][l] /= volumes_4d[i][j][k][l];
     367             :             else
     368             :             {
     369         312 :               _positions_4d[i][j][k].erase(_positions_4d[i][j][k].begin() + l);
     370         312 :               volumes_4d[i][j][k].erase(volumes_4d[i][j][k].begin() + l);
     371         312 :               l--;
     372             :             }
     373             :           }
     374             : 
     375             :   // Fill the 1D position vector
     376          78 :   unrollMultiDPositions();
     377          78 :   _initialized = true;
     378          78 : }
     379             : 
     380             : dof_id_type
     381      393600 : 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      393600 :   if (use_subdomains)
     386      179200 :     return elem.subdomain_id();
     387             :   else
     388      214400 :     return elem.get_extra_integer(id_index);
     389             : }