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 "NodalPatchRecoveryBase.h"
      11             : #include "MathUtils.h"
      12             : 
      13             : #include <Eigen/Dense>
      14             : 
      15             : // TIMPI includes
      16             : #include "timpi/communicator.h"
      17             : #include "timpi/parallel_sync.h"
      18             : #include "libmesh/parallel_eigen.h"
      19             : 
      20             : InputParameters
      21       14290 : NodalPatchRecoveryBase::validParams()
      22             : {
      23       14290 :   InputParameters params = ElementUserObject::validParams();
      24             : 
      25       14290 :   MooseEnum orders("CONSTANT FIRST SECOND THIRD FOURTH");
      26       14290 :   params.addRequiredParam<MooseEnum>(
      27             :       "patch_polynomial_order",
      28             :       orders,
      29             :       "Polynomial order used in least squares fitting of material property "
      30             :       "over the local patch of elements connected to a given node");
      31             : 
      32       14290 :   params.addRelationshipManager("ElementSideNeighborLayers",
      33             :                                 Moose::RelationshipManagerType::ALGEBRAIC,
      34           0 :                                 [](const InputParameters &, InputParameters & rm_params)
      35          36 :                                 { rm_params.set<unsigned short>("layers") = 2; });
      36             : 
      37       14290 :   params.addParamNamesToGroup("patch_polynomial_order", "Advanced");
      38             : 
      39       28580 :   return params;
      40       14290 : }
      41             : 
      42          13 : NodalPatchRecoveryBase::NodalPatchRecoveryBase(const InputParameters & parameters)
      43             :   : ElementUserObject(parameters),
      44          13 :     _qp(0),
      45          13 :     _patch_polynomial_order(
      46          13 :         static_cast<unsigned int>(getParam<MooseEnum>("patch_polynomial_order"))),
      47          13 :     _multi_index(MathUtils::multiIndex(_mesh.dimension(), _patch_polynomial_order)),
      48          26 :     _q(_multi_index.size())
      49             : {
      50          13 : }
      51             : 
      52             : Real
      53        4840 : NodalPatchRecoveryBase::nodalPatchRecovery(const Point & x,
      54             :                                            const std::vector<dof_id_type> & elem_ids) const
      55             : {
      56             :   // Before we go, check if we have enough sample points for solving the least square fitting
      57        4840 :   if (_q_point.size() * elem_ids.size() < _q)
      58           0 :     mooseError("There are not enough sample points to recover the nodal value, try reducing the "
      59             :                "polynomial order or using a higher-order quadrature scheme.");
      60             : 
      61             :   // Assemble the least squares problem over the patch
      62        4840 :   RealEigenMatrix A = RealEigenMatrix::Zero(_q, _q);
      63        4840 :   RealEigenVector b = RealEigenVector::Zero(_q);
      64       21320 :   for (auto elem_id : elem_ids)
      65             :   {
      66       16480 :     A += libmesh_map_find(_Ae, elem_id);
      67       16480 :     b += libmesh_map_find(_be, elem_id);
      68             :   }
      69             : 
      70             :   // Solve the least squares fitting
      71        4840 :   RealEigenVector coef = A.completeOrthogonalDecomposition().solve(b);
      72             : 
      73             :   // Compute the fitted nodal value
      74        4840 :   RealEigenVector p = evaluateBasisFunctions(x);
      75        9680 :   return p.dot(coef);
      76        4840 : }
      77             : 
      78             : RealEigenVector
      79       20840 : NodalPatchRecoveryBase::evaluateBasisFunctions(const Point & q_point) const
      80             : {
      81       20840 :   RealEigenVector p(_q);
      82             :   Real polynomial;
      83       83360 :   for (unsigned int r = 0; r < _multi_index.size(); r++)
      84             :   {
      85       62520 :     polynomial = 1.0;
      86             :     mooseAssert(_multi_index[r].size() == _mesh.dimension(), "Wrong multi-index size.");
      87      187560 :     for (unsigned int c = 0; c < _multi_index[r].size(); c++)
      88      166720 :       for (unsigned int p = 0; p < _multi_index[r][c]; p++)
      89       41680 :         polynomial *= q_point(c);
      90       62520 :     p(r) = polynomial;
      91             :   }
      92       20840 :   return p;
      93           0 : }
      94             : 
      95             : void
      96          60 : NodalPatchRecoveryBase::initialize()
      97             : {
      98          60 :   _Ae.clear();
      99          60 :   _be.clear();
     100          60 : }
     101             : 
     102             : void
     103        4000 : NodalPatchRecoveryBase::execute()
     104             : {
     105        4000 :   RealEigenMatrix Ae = RealEigenMatrix::Zero(_q, _q);
     106        4000 :   RealEigenVector be = RealEigenVector::Zero(_q);
     107       20000 :   for (_qp = 0; _qp < _qrule->n_points(); _qp++)
     108             :   {
     109       16000 :     RealEigenVector p = evaluateBasisFunctions(_q_point[_qp]);
     110       16000 :     Ae += p * p.transpose();
     111       16000 :     be += computeValue() * p;
     112       16000 :   }
     113             : 
     114        4000 :   dof_id_type elem_id = _current_elem->id();
     115        4000 :   _Ae[elem_id] = Ae;
     116        4000 :   _be[elem_id] = be;
     117        4000 : }
     118             : 
     119             : void
     120           5 : NodalPatchRecoveryBase::threadJoin(const UserObject & uo)
     121             : {
     122           5 :   const auto & npr = static_cast<const NodalPatchRecoveryBase &>(uo);
     123           5 :   _Ae.insert(npr._Ae.begin(), npr._Ae.end());
     124           5 :   _be.insert(npr._be.begin(), npr._be.end());
     125           5 : }
     126             : 
     127             : void
     128          55 : NodalPatchRecoveryBase::finalize()
     129             : {
     130             :   // When calling nodalPatchRecovery, we may need to know _Ae and _be on algebraically ghosted
     131             :   // elements. However, this userobject is only run on local elements, so we need to query those
     132             :   // information from other processors in this finalize() method.
     133             : 
     134             :   // Populate algebraically ghosted elements to query
     135          55 :   std::unordered_map<processor_id_type, std::vector<dof_id_type>> query_ids;
     136          55 :   const ConstElemRange evaluable_elem_range = _fe_problem.getEvaluableElementRange();
     137        4555 :   for (auto elem : evaluable_elem_range)
     138        4500 :     if (elem->processor_id() != processor_id())
     139         500 :       query_ids[elem->processor_id()].push_back(elem->id());
     140             : 
     141             :   typedef std::pair<RealEigenMatrix, RealEigenVector> AbPair;
     142             : 
     143             :   // Answer queries received from other processors
     144          30 :   auto gather_data = [this](const processor_id_type /*pid*/,
     145             :                             const std::vector<dof_id_type> & elem_ids,
     146        1000 :                             std::vector<AbPair> & ab_pairs)
     147             :   {
     148         530 :     for (const auto i : index_range(elem_ids))
     149             :     {
     150         500 :       const auto elem_id = elem_ids[i];
     151         500 :       ab_pairs.emplace_back(libmesh_map_find(_Ae, elem_id), libmesh_map_find(_be, elem_id));
     152             :     }
     153          30 :   };
     154             : 
     155             :   // Gather answers received from other processors
     156          30 :   auto act_on_data = [this](const processor_id_type /*pid*/,
     157             :                             const std::vector<dof_id_type> & elem_ids,
     158        1000 :                             const std::vector<AbPair> & ab_pairs)
     159             :   {
     160         530 :     for (const auto i : index_range(elem_ids))
     161             :     {
     162         500 :       const auto elem_id = elem_ids[i];
     163         500 :       const auto & [Ae, be] = ab_pairs[i];
     164         500 :       _Ae[elem_id] = Ae;
     165         500 :       _be[elem_id] = be;
     166             :     }
     167          30 :   };
     168             : 
     169          55 :   libMesh::Parallel::pull_parallel_vector_data<AbPair>(
     170          55 :       _communicator, query_ids, gather_data, act_on_data, 0);
     171          55 : }