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 "NodalPatchRecovery.h"
      11             : #include "SwapBackSentinel.h"
      12             : 
      13             : InputParameters
      14           0 : NodalPatchRecovery::validParams()
      15             : {
      16           0 :   InputParameters params = AuxKernel::validParams();
      17           0 :   MooseEnum orders("CONSTANT FIRST SECOND THIRD FOURTH");
      18           0 :   params.addParam<MooseEnum>("patch_polynomial_order",
      19             :                              orders,
      20             :                              "Polynomial order used in least squares fitting of material property "
      21             :                              "over the local patch of elements connected to a given node");
      22           0 :   params.addParamNamesToGroup("patch_polynomial_order", "Advanced");
      23             : 
      24             :   // TODO make this class work with relationship manager
      25             :   // params.registerRelationshipManagers("ElementSideNeighborLayers");
      26             :   // params.addPrivateParam<unsigned int short>("element_side_neighbor_layers", 2);
      27             : 
      28           0 :   return params;
      29           0 : }
      30             : 
      31           0 : NodalPatchRecovery::NodalPatchRecovery(const InputParameters & parameters)
      32             :   : AuxKernel(parameters),
      33           0 :     _patch_polynomial_order(
      34           0 :         isParamValid("patch_polynomial_order")
      35           0 :             ? static_cast<unsigned int>(getParam<MooseEnum>("patch_polynomial_order"))
      36           0 :             : static_cast<unsigned int>(_var.order())),
      37           0 :     _fe_problem(*parameters.get<FEProblemBase *>("_fe_problem_base")),
      38           0 :     _multi_index(multiIndex(_mesh.dimension(), _patch_polynomial_order))
      39             : {
      40             :   // if the patch polynomial order is lower than the variable order,
      41             :   // it is very likely that the patch recovery is not used at its max accuracy
      42           0 :   if (_patch_polynomial_order < static_cast<unsigned int>(_var.order()))
      43           0 :     mooseWarning("Specified 'patch_polynomial_order' is lower than the AuxVariable's order");
      44             : 
      45             :   // TODO remove the manual ghosting once relationship manager is working correctly
      46             :   // no need to ghost if this aux is elemental
      47           0 :   if (isNodal())
      48             :   {
      49           0 :     MeshBase & meshhelper = _mesh.getMesh();
      50           0 :     meshhelper.allow_renumbering(false);
      51           0 :     for (const auto & elem :
      52           0 :          as_range(meshhelper.semilocal_elements_begin(), meshhelper.semilocal_elements_end()))
      53           0 :       _fe_problem.addGhostedElem(elem->id());
      54             :   }
      55           0 : }
      56             : 
      57             : std::vector<std::vector<unsigned int>>
      58           0 : NodalPatchRecovery::nChooseK(unsigned int N, unsigned int K)
      59             : {
      60           0 :   std::vector<std::vector<unsigned int>> n_choose_k;
      61           0 :   std::vector<unsigned int> row;
      62           0 :   std::string bitmask(K, 1); // K leading 1's
      63           0 :   bitmask.resize(N, 0);      // N-K trailing 0's
      64             : 
      65             :   do
      66             :   {
      67           0 :     row.clear();
      68           0 :     row.push_back(0);
      69           0 :     for (unsigned int i = 0; i < N; ++i) // [0..N-1] integers
      70           0 :       if (bitmask[i])
      71           0 :         row.push_back(i + 1);
      72           0 :     row.push_back(N + 1);
      73           0 :     n_choose_k.push_back(row);
      74           0 :   } while (std::prev_permutation(bitmask.begin(), bitmask.end()));
      75             : 
      76           0 :   return n_choose_k;
      77           0 : }
      78             : 
      79             : std::vector<std::vector<unsigned int>>
      80           0 : NodalPatchRecovery::multiIndex(unsigned int dim, unsigned int order)
      81             : {
      82             :   // first row all zero
      83           0 :   std::vector<std::vector<unsigned int>> multi_index;
      84           0 :   std::vector<std::vector<unsigned int>> n_choose_k;
      85           0 :   std::vector<unsigned int> row(dim, 0);
      86           0 :   multi_index.push_back(row);
      87             : 
      88           0 :   if (dim == 1)
      89           0 :     for (unsigned int q = 1; q <= order; q++)
      90             :     {
      91           0 :       row[0] = q;
      92           0 :       multi_index.push_back(row);
      93             :     }
      94             :   else
      95           0 :     for (unsigned int q = 1; q <= order; q++)
      96             :     {
      97           0 :       n_choose_k = nChooseK(dim + q - 1, dim - 1);
      98           0 :       for (unsigned int r = 0; r < n_choose_k.size(); r++)
      99             :       {
     100           0 :         row.clear();
     101           0 :         for (unsigned int c = 1; c < n_choose_k[0].size(); c++)
     102           0 :           row.push_back(n_choose_k[r][c] - n_choose_k[r][c - 1] - 1);
     103           0 :         multi_index.push_back(row);
     104             :       }
     105             :     }
     106             : 
     107           0 :   return multi_index;
     108           0 : }
     109             : 
     110             : void
     111           0 : NodalPatchRecovery::computePVector(Point q_point)
     112             : {
     113           0 :   _P.resize(_multi_index.size());
     114             :   Real polynomial;
     115           0 :   for (unsigned int r = 0; r < _multi_index.size(); r++)
     116             :   {
     117           0 :     polynomial = 1.0;
     118           0 :     for (unsigned int c = 0; c < _multi_index[0].size(); c++)
     119           0 :       for (unsigned int p = 0; p < _multi_index[r][c]; p++)
     120           0 :         polynomial *= q_point(c);
     121           0 :     _P(r) = polynomial;
     122             :   }
     123           0 : }
     124             : 
     125             : void
     126           0 : NodalPatchRecovery::accumulateAMatrix()
     127             : {
     128           0 :   DenseMatrix<Number> PxP;
     129           0 :   PxP.resize(_multi_index.size(), _multi_index.size());
     130           0 :   PxP.outer_product(_P, _P);
     131           0 :   _A.add(1.0, PxP);
     132           0 : }
     133             : 
     134             : void
     135           0 : NodalPatchRecovery::accumulateBVector(Real val)
     136             : {
     137           0 :   _B.add(val, _P);
     138           0 : }
     139             : 
     140             : void
     141           0 : NodalPatchRecovery::reinitPatch()
     142             : {
     143             :   // clean and resize P, A, B, coef
     144           0 :   _P.resize(_multi_index.size());
     145           0 :   _A.resize(_multi_index.size(), _multi_index.size());
     146           0 :   _B.resize(_multi_index.size());
     147           0 :   _coef.resize(_multi_index.size());
     148             : 
     149             :   // activate dependent material properties
     150           0 :   std::unordered_set<unsigned int> needed_mat_props;
     151           0 :   const auto & mp_deps = getMatPropDependencies();
     152           0 :   needed_mat_props.insert(mp_deps.begin(), mp_deps.end());
     153           0 :   _fe_problem.setActiveMaterialProperties(needed_mat_props, _tid);
     154           0 : }
     155             : 
     156             : void
     157           0 : NodalPatchRecovery::compute()
     158             : {
     159             :   // Fall back on standard procedure for element variables
     160           0 :   if (!isNodal())
     161             :   {
     162           0 :     AuxKernel::compute();
     163           0 :     return;
     164             :   }
     165             : 
     166             :   // Limit current use of NodalPatchRecovery to a single processor
     167           0 :   if (_communicator.size() > 1)
     168           0 :     mooseError("The nodal patch recovery option, which calculates the Zienkiewicz-Zhu patch "
     169             :                "recovery for nodal variables (family = LAGRANGE), is not currently implemented for "
     170             :                "parallel runs. Run in serial if you must use the nodal patch capability");
     171             : 
     172             :   // Use Zienkiewicz-Zhu patch recovery for nodal variables
     173           0 :   reinitPatch();
     174             : 
     175             :   // get node-to-conneted-elem map
     176           0 :   const std::map<dof_id_type, std::vector<dof_id_type>> & node_to_elem_map = _mesh.nodeToElemMap();
     177           0 :   auto node_to_elem_pair = node_to_elem_map.find(_current_node->id());
     178             :   mooseAssert(node_to_elem_pair != node_to_elem_map.end(), "Missing entry in node to elem map");
     179           0 :   std::vector<dof_id_type> elem_ids = node_to_elem_pair->second;
     180             : 
     181             :   // consider the case for corner node
     182           0 :   if (elem_ids.size() == 1)
     183             :   {
     184           0 :     dof_id_type elem_id = elem_ids[0];
     185           0 :     for (auto & n : _mesh.elemPtr(elem_id)->node_ref_range())
     186             :     {
     187           0 :       node_to_elem_pair = node_to_elem_map.find(n.id());
     188           0 :       std::vector<dof_id_type> elem_ids_candidate = node_to_elem_pair->second;
     189           0 :       if (elem_ids_candidate.size() > elem_ids.size())
     190           0 :         elem_ids = elem_ids_candidate;
     191           0 :     }
     192             :   }
     193             : 
     194             :   // before we go, check if we have enough sample points for solving the least square fitting
     195           0 :   if (_q_point.size() * elem_ids.size() < _multi_index.size())
     196           0 :     mooseError("There are not enough sample points to recover the nodal value, try reducing the "
     197             :                "polynomial order or using a higher-order quadrature scheme.");
     198             : 
     199             :   // general treatment for side nodes and internal nodes
     200           0 :   for (auto elem_id : elem_ids)
     201             :   {
     202           0 :     const Elem * elem = _mesh.elemPtr(elem_id);
     203             : 
     204           0 :     if (!elem->is_semilocal(_mesh.processor_id()))
     205           0 :       mooseError("skipped non local elem!");
     206             : 
     207           0 :     _fe_problem.prepare(elem, _tid);
     208           0 :     _fe_problem.reinitElem(elem, _tid);
     209             : 
     210             :     // Set up Sentinel class so that, even if reinitMaterials() throws, we
     211             :     // still remember to swap back during stack unwinding.
     212           0 :     SwapBackSentinel sentinel(_fe_problem, &FEProblem::swapBackMaterials, _tid);
     213           0 :     _fe_problem.reinitMaterials(elem->subdomain_id(), _tid);
     214             : 
     215           0 :     for (_qp = 0; _qp < _q_point.size(); _qp++)
     216             :     {
     217           0 :       computePVector(_q_point[_qp]);
     218           0 :       accumulateAMatrix();
     219           0 :       accumulateBVector(computeValue());
     220             :     }
     221           0 :   }
     222             : 
     223           0 :   _fe_problem.clearActiveMaterialProperties(_tid);
     224             : 
     225             :   // solve for the coefficients of least square fitting
     226             :   // as we are on the physical coordinates, A is not always SPD
     227           0 :   _A.svd_solve(_B, _coef);
     228             : 
     229             :   // compute fitted nodal value
     230           0 :   computePVector(*_current_node);
     231           0 :   Real nodal_value = _P.dot(_coef);
     232             : 
     233             :   // set nodal value
     234           0 :   _fe_problem.reinitNode(_current_node, _tid);
     235           0 :   dof_id_type dof = _var.nodalDofIndex();
     236             :   {
     237           0 :     Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
     238           0 :     _solution.set(dof, nodal_value);
     239           0 :   }
     240           0 :   _var.setNodalValue(nodal_value);
     241           0 : }